STEPS: Slip Time Earthquake Path Simulations Applied to the San Andreas and Toe Jam Hill Faults to Redefine Geologic Slip Rate Uncertainty

Geologic slip rates are a time‐averaged measurement of fault displacement calculated over hundreds to million‐year time scales and are a primary input for probabilistic seismic hazard analyses, which forecast expected ground shaking in future earthquakes. Despite their utility for seismic hazard calculations, longer‐term geologic slip rates represent a time‐averaged measure of the tempo of strain release and do not measure variability across earthquake cycles. We have developed a numerical approach called STEPS (Slip Time Earthquake Path Simulations), which is built upon field‐based observations and explicitly incorporates realistic variations in displacement per event and variability in the recurrence interval between earthquakes. The STEPS approach, which simulates strain release through time, relies on representing earthquake cycles as stairsteps, rather than straight‐line paths, connecting per earthquake time‐displacement coordinates. We simulate earthquake histories based on these input constraints using two examples: the Carrizo section of the San Andreas fault and the Toe Jam Hill fault of the Seattle fault zone. We find that modeled slip rate distributions agree with slip rates reported for the sites of interest by the original investigators, while providing a slip rate distribution that reflects the variability of earthquake frequency and displacement. The STEPS approach provides an estimate of fault slip rate uncertainty based on a synthetic suite of plausible time‐displacement paths resulting from individual earthquakes, rather than measurement uncertainties associated with offset features. When considering this simulated earthquake behavior between measurements, the uncertainty associated with earthquake paths is greater than that calculated by the long‐term rate.

, and can be measured as either horizontal or vertical separations. Offset features are assigned ages using either relative or absolute methods such as radiocarbon (e.g., Cowan, 1990), luminescence (e.g., Delano et al., 2017), and cosmogenic radionuclides (e.g., Le et al., 2007). In the absence of absolute geochronometers, relative methodologies such as scarp morphology (e.g., Hanks & Schwartz, 1987), soil profile analysis (e.g., Rockwell et al., 1984), regional Quaternary stratigraphic/geomorphic correlations (e.g., Koehler, 2019), bedrock relationships (e.g., Wells et al., 2020), and weathering rind thicknesses (e.g., Knuepfer, 1992), among other methods, may be applied. Most ideally, the age associated with a given offset feature is close in time to either the genesis or abandonment of that feature, indicating the time after which the feature accrued offset. Ideally, the genesis or abandonment age of an offset feature would be bracketed by multiple ages that pre-and postdate the time interval of interest (e.g., Blisniuk et al., 2012). Geologic slip rates are typically made at a discrete geographic location, given that offset landforms are observed at a specific location and the age constraining that offset comes from that feature or a nearby deposit. Although geologic slip rates are constrained at points, the rates are often generalized to describe the behavior of a fault segment or even an entire fault. A general assumption is that geologic slip rates represent primarily coseismic offset, but also possibly an unknown and minor contribution from postseismic slip or interseismic creep.
In this study, we recalculate slip rate uncertainties by incorporating singe event variations in recurrence interval and displacements, moving toward the ability to capture the full range of complex fault behaviors for both active tectonics and seismic hazard studies. To this end, we develop a modeling approach called STEPS (Slip Time Earthquake Path Simulations). STEPS is a Monte Carlo simulation that incorporates the variation among possible time-displacement paths (i.e., earthquake paths) when considering a distribution of possible recurrence intervals and single event displacements. STEPS leverages known information such as the timing and displacement of the most recent earthquake, multiple dated offset markers if available, and the fact that current displacement is 0 m (e.g., the origin of a displacement-time plot in a Cartesian coordinate system). With these input data, STEPS then builds a simulated set of possible earthquake paths that expand the mean rate provided by geologic slip rates. This process allows the attribution of slip rate uncertainty using variability in recurrence interval and single event displacement in addition to measurement uncertainties. Once earthquake paths are constructed, they can be extrapolated further into the past, or in the future; the future extrapolation allows for construction of a probability distribution for the timing of the next event, a key parameter for seismic hazard analyses.

Previous Work in Understanding the Representation of Geologic Slip Rates
Geologic slip rate uncertainty distributions, temporal variability, and implementation within seismic hazard models have been explored in previous studies. Bird (2007) explored uncertainty bounds of geologic slip rates in the western United States utilizing geophysical and geologic inferences and fault-specific convolutions of age and offset distributions and concluded that very few faults have "well-constrained" slip rates. In that work, Bird (2007) defines well-constrained slip rates as those where the width of the 95% confidence interval is less than the median slip rate. Thompson et al. (2002) and Zechar and Frankel (2009) advocate for utilizing a distribution for a given offset, as well as the associated age, to arrive at a site-specific distribution that provides a more process-informed view of geologic slip rate. Following this work, Gold and Cowgill (2011) developed a Monte Carlo process by which they describe fault slip histories by utilizing distributions in both age and offset space for multiple offset features, and then sampling locations within 10.1029/2021GC009848 3 of 29 these distributions to generate numerous simulations of slip histories between the field data. Other studies have examined how per-event variability can influence slip rate estimates. For example, Styron (2019) completed numerical simulations of faults with variable recurrence intervals, testing various distributions of earthquake return time, and found that the distribution shape did not impact slip rate variability through time, but the coefficient of variation (CV) of those distributions and number of earthquakes over which slip rates are averaged had large impacts on the likelihood that a given sample of observed earthquakes would yield an accurate long-term slip rate. Fitzenz et al. (2010) explored probabilistically how single-event distributions combine to match geologic slip rates, though their focus was on the trimmed per event posteriors rather than the resulting slip rate variability.
Geologic slip rate is a key input for seismic hazard assessments and, aside from ground motion modeling, places a first-order constraint on seismic hazard at a given location (e.g., Field et al., 2013;Petersen et al., 2020;Stirling et al., 2012;Wesnousky et al., 1984). In an analysis of the relationship between background seismicity rates and Quaternary geologic slip rates on intraplate faults in Japan, Wesnousky et al. (1984) established that geologic slip rate is a useful and reliable proxy for earthquake rate in seismic hazard calculations when compared to seismicity. This comparison demonstrates that, with all else held equal, an increase in geologic slip rate leads to an increase in calculated seismic hazard. However, if a geologic slip rate is known to vary by a factor of 4, as was observed on the Wellington fault (New Zealand) over 100-to 100,000-year time scales (Ninis et al., 2013), the resultant hazard calculation will depend on which slip rate is used. Van Dissen et al. (2020) showed that implementing the faster-than-average slip rate along the Wellington fault (19 mm/yr) increased peak ground acceleration (PGA) (2% exceedance in 50 years) in the city of Wellington by 50% compared to PGA when calculated using the slower-than-average slip rate (1.3 mm/yr). In this case, by holding all other factors equal (i.e., fault geometry, ground motion parameters, basin thickness) and only changing slip rate along the Wellington fault, slip rate has a large influence on resultant PGA. Temporally variable slip rates can be accommodated in hazard computations through different logic tree branches. For example, in the U.S. National Seismic Hazard Map (NSHM), the New Madrid seismic zone is treated with the possibility of being in sequence or out of sequence (i.e., seismic cluster or seismic lull) or not having earthquake clustering at all (Petersen et al., 2015(Petersen et al., , 2020. However, despite the known limitations of using a single integrated slip rate to infer seismic hazard, more work is needed on how to best measure, report, and utilize slip rate uncertainties.

Challenges With Geologic Slip Rates
We identify and discuss four challenges with the typical methodologies of measuring, calculating, and reporting geologic slip rates that we explore in the STEPS approach: (a) linear approximation of earthquake cycle behavior; (b) averaging slip rate through the origin; (c) number of seismic cycles included in a slip rate; and (d) uncertainty determination and reporting ( Figure 1).

Linear Fits of Earthquake Behavior
Geologic slip rates are intended to describe time-averaged earthquake behavior measured at a specific location. Despite earthquakes occurring as punctuated events in time with intervening interseismic cycles of variable lengths, a straight line with constant slope connecting two observation points is commonly used to calculate slip rates (i.e., mean slip rate). These two points are, in many cases, one observation of a dated displaced feature in the landscape and the origin of a displacement-time plot (Figure 1a). Often, two straight lines illustrate the maximum average slip rate (calculated from the minimum time, maximum displacement coordinate of the observation box to the origin) and the minimum average slip rate (calculated from the maximum time, minimum displacement coordinate of the observation box to the origin) permitted from the data. A benefit to this approach is it establishes which faults move slower or faster on average and requires little to no interpretation, and typically only the limiting uncertainties in age and offset are used. This approach does not capture variation in recurrence and displacement at the level of individual earthquake cycles. Over time, earthquake paths may deviate from the swath of slip rates bracketed by the minimum and maximum. Slip rates calculated along earthquake paths over shorter and longer time intervals can be greater than and less than the apparent maximum and minimum slip rates from the swath outline, respectively. Earthquake paths indicate that variability in slip rate through time is not captured with the simple approximation of earthquake occurrence through time as a set of single lines with constant slope (e.g., Wallace, 1987).

Averaging Slip Rate Through the Origin
Slip rates are typically averaged through the present day, which is the origin of a displacement-time plot (Figure 1). The origin represents the most certain observation we have regarding fault behavior, showing the lack of a surface-rupturing earthquake occurring on the day a given displacement-time plot is generated. However, a drawback of this convenience is that it can impose a bias on slip rates, particularly for slow strain-rate faults (e.g., Weldon & Sieh, 1985). As the open interval continues, the origin is mobile along the time (x) axis, such that the observation of 0 m displacement remains constant while time advances to the left (Figure 1b). All slip rates-both minimum and maximum bounds-only get slower as time goes on if the singular point of the origin is used in slip rate calculations. Given the direct correlation between slip rate and seismic hazard (e.g., Van Dissen et al., 2020;Wesnousky et al., 1984), this observation implies hazard also decreases as time passes, which may be an unlikely assumption (e.g., Matthews et al., 2002;Milner et al., 2020;Salditch et al., 2020). Such a bias is particularly pronounced on low strain rate faults (i.e., <1 mm/yr), because the chance of the open interval lasting many hundreds to thousands of years is far greater than that same possibility on a high strain rate fault (i.e., >5 mm/yr). Although some studies consider the strain accumulated during the formation of the displaced geologic marker and the most recent in upper right-hand side of plot) with swath (light gray diagonal polygonal) connecting observation to the origin. A hypothetical earthquake path is shown as a blue stairway with individual earthquakes marked by stars. (b) As in A, with additional illustration of "the origin problem." Yellow circles paired with gray swaths illustrate a decreasing slip rate as the length of the open interval continues. (c) Schematic figure highlighting the difference in number of events averaged over a hypothetical high strain rate fault (purple earthquake path/left y-axis) compared to a hypothetical slow strain rate fault (green earthquake path/right y-axis). (d) Three representations of geologic slip rate uncertainty. i: Single value of slip rate with no reported uncertainty. ii: Uniform distribution of slip rate, with minimum, preferred (mid-point) and maximum rates given equal probability. iii: Unique distribution of slip rate that reports a median, 68% confidence interval (CI), and 95% CI.
The disparity in number of events averaged across a geologic slip rate may have large implications for estimating slip rate uncertainty along a given fault (e.g., Styron, 2019). In the case of the high strain rate fault (e.g., San Andreas fault), many events are integrated into the displacement of a single marker from the end of the last glacial maximum. In contrast, a few-event record from a low strain rate fault (e.g., Toe Jam Hill fault) may overrepresent one mode of strain release (whether a cluster of earthquakes or a seismic lull). In either case, the utility of these slip rate estimates may be limited when trying to understand the spectrum of possible fault behaviors. However, geologists must draw the best possible conclusions based on the available data, and rates are reported with variable numbers of events given permissible field observations.

Uncertainty Determination and Reporting
Geologic slip rates may be reported as a single value, as a range with a minimum, maximum, and mid-point preferred value (i.e., boxcar distribution), or a continuous distribution derived from measurement uncertainty distributions (Figure 1d). When geologic slip rate is reported as a single value, there is typically no uncertainty in the estimates of displacement or age (Figure 1d). Given the lack of measurement uncertainty in this case, reporting slip rate as a single value is not satisfactory. A uniform distribution of slip rates is an improvement on reporting only a single value. An appropriate use of this distribution is when there is no preferred restoration nor preferred age of a given displaced feature, and only a minimum or maximum can be determined for either time or displacement or both.
Offset and age measurements used to calculate slip rates are ideally reported with measurement errors that capture the range of permissible offset feature reconstructions and the analytical uncertainties of offset feature ages. Increasingly, slip rate studies incorporate complex distributions of both displacement and time measurements and produce a unique probability distribution for slip rate (Figure 1d) (e.g., Burgette et al., 2020;Nicol et al., 2016;Zechar & Frankel, 2009;Zinke et al., 2019). When deriving this unique distribution, the median, 68% confidence interval (CI) and 95% CI may also be reported. It is worth noting that reduction of uncertainties in displacement and time (i.e., more precise measurements of offsets and ages) tend to lead to a narrowing of mean slip rate uncertainty. While precision in displacement-time measurements is valued and warranted, the potential overinterpretation of constant fault behavior between, for instance, a very precise displacement-time measurement and the origin may not be supported. Overly precise displacement-time measurements, though, may underestimate the process uncertainty of how a given feature formed or consider fault zone evolution through time (e.g., Reitman et al., 2019). For instance, at the Wallace Creek site on the San Andreas fault, Sieh and Jahns (1984) report uncertainty of ±3 m on an offset measurement of 475 m (<1% measurement uncertainty), a very small margin of error on this offset channel. A related issue is that the age derived for a landform, even if very precise, may significantly pre-or postdate the time over which the marker accumulated displacement (e.g., DuRoss et al., 2018DuRoss et al., , 2020McGill et al., 2009). Overly precise measurements may result in slip rate uncertainties smaller than justifiable by the resolution of geomorphic recorders or geochronologic methods, given that relating absolute ages of a landform to the actual age of a landform remains an outstanding challenge in active tectonics.
It is also unclear if slip rate is sufficiently described with a mean and standard deviation through time given a potential continuum of fault behavior from strictly periodic to bursty strain release (e.g., Chen et al., 2020;Cowie et al., 2012;Mouslopoulou et al., 2009;Salditch et al., 2020). On a mechanistic level, measurement uncertainties represent our uncertainty in the earthquake cycle itself (i.e., process uncertainty). These complications aside, reported geologic slip rate uncertainty only captures uncertainty in displacement and age measurements.

Potential STEPS Forward
Previous studies focus on displaced features but generally give little attention to the displacement-time space between dated offset observations, that is, the potential paths between observations. The work we describe here is designed to interpolate between observations points by randomly sampling a synthetically derived, geologically based catalog of surface-rupturing earthquakes, without any earthquake physics or renewal models imposed on the sampling process. To address the four challenges described above, we have developed a Monte Carlo simulation that provides a technique for addressing these distinct issues, called STEPS (Slip Time Earthquake Path Simulations). Our approach models displacement-time paths (i.e., earthquake paths) visualized as stairs between published field measurements of dated offset features. This model addresses the above four challenges in the following ways: 1. The space between dated offset observations is filled with earthquake paths plotted as stairs, not as lines with a constant slope, that allow for the interpolation of simulated earthquakes between observations. 2. Slip rate calculations are averaged through the displacement-time observation of the most recent event (MRE) and not the origin. 3. Utilizing constant sampling intervals along earthquake paths provides slip rates that have less bias due to an inconsistent number of events comprising a given slip rate. 4. By conducting a large number of realizations of this simulation, we explore geologic slip rate uncertainty that encompasses variability of the earthquake-level behavior, as opposed to a long-term mean slip rate derived from only measurement uncertainty.
We present the STEPS approach in detail and apply the method to field observations from a high strain rate fault (Carrizo section of the San Andreas fault) and a low strain rate fault (Toe Jam Hill fault of the Seattle fault zone). We show that STEPS-derived slip rate distributions agree with previously reported preferred slip-rate values. We derive and report the 68 and 95% confidence intervals (CI) of these STEPS-derived rates. We test the results against field-observed paleoseismic chronologies. Finally, we present one potential application of this numerical approach in deriving a probability distribution for the timing of the next future event.

Methods
The methodology we use is an analytical numerical solution that generates a distribution of earthquake paths utilizing predefined distributions of single event displacement and recurrence interval. The numerical model is written in MATLAB and is available on GitLab at https://doi.org/10.5066/P9NZKRXH. In the following sections, we outline the STEPS method, discuss results of the method imposed on a hypothetical fault to briefly explore the parameter space, and describe the field-based input data from the San Andreas and Toe Jam Hill faults used for the results section of this manuscript.

The STEPS Approach
The STEPS approach relies on the following input data for a given site along a fault (Figures 2a and 2b): 1. Displacement-time measurements (i.e., dated offsets typically used for long-term mean slip rate estimates) 2. Timing and displacement of the most recent event 3. Distributions of recurrence interval (RI) and single event displacement (SED) Given these input data, the STEPS approach provides the basic output as a set of earthquake paths that satisfy the displacement-time measurements (Figure 2c). From the results, many derivative products can be calculated, including a distribution of slip rates along these earthquake paths (Figures 2d-2f). The goal of this approach is to provide slip rates recalculated from field measurements over specific time intervals meaningful for seismic hazard analysis.
In the STEPS approach, input data in the form of displacement-time boxes are plotted, similar to most existing slip rate estimation techniques ( Figure 2a). The displacement-time observations, made in the field as a measurement of an offset geomorphic feature (see Introduction for more details), are represented as uniform distributions of uncertainty in the measurements of a given dated offset feature. Data pertaining to the most recent event are represented in the same manner as older displacement-time measurements. Additionally, applications of the STEPS approach could include unique uncertainty distributions to the input displacement-time observations. The minimum number of displacement-time observations for the STEPS approach without utilizing the origin as an observation point is two, in order to initiate the earthquake paths in one box and terminate them in another.
Next, earthquake paths that connect displacement-time observations are generated using predefined sampling distributions (Figure 2b). In Figure 2, we present possible distributions of SED and RI as normal and exponential distributions, respectively. The normal distribution of SED is selected to represent the variability of slip along rupture length in a given earthquake (e.g., Zielke et al., 2010). We use the normal distribution as a more general unimodal form of the lognormal distribution derived by Hecker et al. (2013) to allow for more symmetric displacements at the tails. The choice of the normal distribution for SED broadly incorporates the lognormal SED distribution preferred by Hecker et al. (2013) without assuming that the displacements in the most recent few events are inclusive of the entire long-term behavior of a fault. The exponential distribution of RI is selected to represent the assumption that earthquake occurrence through time is random (CV = 1) (e.g., Field et al., 2013;Kagan & Jackson, 1991;Salditch et al., 2020;Styron, 2019). The STEPS approach is interchangeable such that virtually any distribution can be ascribed to SED and RI if more information would indicate a preferred distribution choice over another. In our discussion of Figures 3 and 4, we demonstrate that the choice of distribution in SED and RI is not as influential on the earthquake path wander as other parameters, such as size of displacement-time observation box. As a test Swath (gray polygon) connects outer edges of displacement-time observations. (b) Sampling distributions for single event displacement (normal) and recurrence interval (exponential). (c) Five hypothetical earthquake paths that successfully fit input data from in panel (a) using sampling distributions in panel (b, d, and e) Illustration of discretized earthquake paths at coarser (d) and finer (e) temporal spacing. Short line segments connecting blue circles (discretization points) are colored red, yellow, and green for relatively slow, moderate, and fast slip rates, respectively. (f) Hypothetical slip rate distribution shown as a histogram in the background and as a Kernel Density Estimation (KDE) in foreground. Median, 68% confidence interval (CI), and 95% CI are labeled.  case, we explore the use of a lognormal distribution in RI with a CV ∼0.7 at the Wallace Creek site in the Discussion Section 4.3 in this manuscript.
Once the sampling distributions of SED and RI are defined, the STEPS approach samples randomly from these distributions to pair the step height (SED) and length (RI) of the earthquake paths. The earthquake paths begin at a location randomly selected within the most recent event displacement-time observation box, and sequentially adds earthquakes further back in time. Once an earthquake path has reached a threshold value in time, the STEPS approach performs a spatial query to test if the included displacement-time observation boxes have been intersected by a given earthquake path. If all the displacement-time observation boxes have been intersected by an earthquake path, then that path is considered successful (Figure 2c). If the earthquake path does not satisfy the displacement-time observations, that particular path is removed from further analysis. The STEPS approach continues iterating until a preset number of earthquake paths has been reached. In the case of results presented in this manuscript, the target number of successful earthquake paths is 600, explored in the subsequent section.
If the timing and displacement of the MRE is unknown, the origin may be used as a substitute in lieu of the MRE. We do not combine the open interval and MRE as constraints given that the open interval may be significantly shorter than the RI sampling distribution would imply-the RI sampling distribution represents closed seismic intervals, and the time interval between the origin and MRE represents an open seismic interval.
The earthquake paths generated by the STEPS approach stray outside the straight-line swaths connecting the displacement-time observations (Figure 1a). The deviation of the earthquake paths from the straight lines connecting the edges of the displacement-time observations is due to the random ordering of samples from the predefined distributions of SED and RI ( Figure 2b). As is further described in the discussion of Figure 3, the deviation from the swath is a natural outcome of exploring the full range of permissible time-displacement combinations and a direct result of random ordering of the picked SED and RI values (e.g., Cowie et al., 2012).
Once a suite of successful paths is defined, the paths are discretized and subsampled at intervals to calculate slip rates (Figures 2d and 2e). The paths are discretized at time intervals as a multiplier of the input mean recurrence interval, and the displacement coordinate along the path is interpolated at these equal-interval time locations. Discretizing an earthquake path at longer time intervals (Figure 2e) yields fewer straight-line segments of constant slope, and therefore a more uniform slip rate through time. Averaging over shorter time intervals ( Figure 2d) yields a larger number of straight-line segments approximating the earthquake path's shape and will more reliably match the path and perhaps yield a wider distribution of slip rates. The discretization values utilized in this study are 5*mean RI and 10*mean RI. We select these discretization values following the observations by Styron (2019), Nicol et al. (2016), and Biasi (2014), which suggest that slip rate variability and uncertainty tends to decrease after 5-10*mean RI. While we have selected 5* and 10*mean RI for this analysis, these intervals can be modulated based on the time horizon of interest (such as 2% in 50-year recurrence of a threshold value). Once successful paths are subsampled for slip rates, all modeled slip rates for a given time interval discretization are fit with a kernel density estimation (KDE) (Figure 2f). From this slip rate distribution, we report the median (50th quantile), 68% confidence interval (16-84th quantile), and 95% confidence interval (2.5-97.fifth quantile).
The STEPS approach provides a self-contained, internally consistent methodology to synthesize different sets field observations across different time scales. While the STEPS approach distills existing information into a package intended for understanding slip rate uncertainty through time, the STEPS approach does not provide a substitute for geologic field-based observations. Instead, the STEPS approach combines these critical observations into a distribution of slip rate uncertainty including earthquake cycle effects, represented in STEPS as random ordering of SED and RI picked from predefined input distributions. Such earthquake cycle effects are accounted for in various ways in applications such as the Uniform California Earthquake Rupture Forecast, version 3 (UCERF3)  and physics-based rupture simulations (e.g., Dieterich and Richards-Dinger, 2010;Milner et al., 2020;Shaw et al., 2018). The STEPS approach provides an independent and geologically based understanding of simulated earthquake histories, and therefore slip rate variability through time that can then be compared to other existing methods that utilize different approaches to understand the same phenomena.

Parameter Exploration Within STEPS
In this section, we present basic parameter testing of the critical inputs to STEPS for illustrative purposes. We present results of 12 arbitrary test cases for a fault with a nominal mean slip rate of 1 mm/yr but vary displacement-time observation box size and SED and RI distribution shape and ranges. The results of these test cases from a hypothetical fault are shown in Figures 3 and 4, with common SED and RI sampling distributions in panels (a-f) in both figures but with smaller displacement-time observation boxes in Figure 3 and larger boxes in Figure 4. In Figures 3 and 4a, we select narrow uniform distributions for SED and RI. In contrast, panel (b) shows an arbitrarily broad uniform distribution for SED and RI. Panels (c and d) explore the coupling of a narrow, uniform distribution in one quantity with a broad, uniform distribution in the other. Finally, panels (e and f) test the use of more nuanced distribution choices in RI (exponential) and SED (normal), but explore the influence of increasing CV from panel (e to f). The goal of this set of comparisons is to understand the relative influence of input data (sampling distributions of SED and RI) and model constraints (displacement-time box sizes) on output slip rate distributions ( Figure 5).
Two basic observations emerge across the 12 model runs (Figures 3 and 4). The first observation is that the modeled earthquake paths overlap to form a distribution that is centered on the shortest connection between displacement-time boxes. The heat maps illustrate this intuitive but basic point. The second basic observation is that the earthquake paths extend beyond the straight-line swath connecting the edges of the displacement-time boxes, regardless of the box size or the sampling distributions. The extension of earthquake paths outside the swath, regardless of parameter choice, is a central difference between the STEPS approach and traditional mean slip rate estimate methods.
The deviation of paths from the straight-line swaths appears, in general, to be insensitive to sampling distributions chosen (Figures 3 and 4). However, if narrow SED and RI sampling distributions are selected with large displacement-time observation boxes (e.g., Figure 4a), the earthquake paths do not deviate from the swath, simply because the (in this case, hypothetical) uncertainty in offset and age measurements is larger than the uncertainty of the input sampling distributions. Figure 4a presents perhaps a rare instance where the distribution of SED and RI is better determined than the dated offset measurement uncertainty. This relationship is best observed between Figures 3c and 3d and Figures 4c and 4d. In Figure 3, the heatmaps for C and D have similar width, despite different input values for SED and RI (sub-panels ii-iii). The impact of changing of input sampling distribution becomes apparent when the displacement-time boxes increase in size (Figures 4c and 4d). Here, we see the earthquake paths in Figure 4c preferentially occupy the 'top' of the swath, and Figure 4d preferentially occupies the 'bottom' of the swath, due to the differences in their input sampling distributions. Sampling larger SED in Figure 4c leads to a preferentially higher rate for the majority of the earthquake paths whereas sampling longer RI in Figure 4d leads to a preferentially lower rate for the earthquake paths.
While sampling distribution choice and displacement-time box size influence the output slip rate distribution shape and associated uncertainty bounds, the median values tend to be similar ( Figure 5). Figures 5a and 5b compares the kernel density estimates (KDEs) of slip rates calculated from the six test cases shown in Figures 3 and 4, respectively. In this case, slip rates are calculated every 5 kyr (e.g., 5*mean RI). The median slip rates values with small displacement-time boxes are tightly clustered about 1 mm/yr (range of median values ∼0.8-1.1 mm/yr), whereas there is a greater spread of median values for slip rates calculated with larger displacement-time boxes (range of median values ∼0.7-1.4 mm/yr). Increasing the displacement-time observation box size unsurprisingly increases the scatter of slip rate distributions when comparing across different SED and RI sampling distributions. However, within either case of small or large displacement-time box size, sampling distribution affects the range of 68% confidence interval of the STEPS calculated slip rates (see shaded regions in Figures 5a and 5b). Distributions with the broadest ranges of SED and RI (e.g., B and F) yield broader 68% confidence intervals, and distributions with narrow ranges of SED and RI (e.g., A) yield narrow distributions of slip rate. Increasing displacement-time box size also produces a skew observed in slip rates for test cases C and D, increasing the slip rate distribution of C far to the right (84th percentile ∼1.7 mm/yr), and decreasing the slip rate distribution D to the left (sixteenth percentile ∼0.4 mm/yr). Further exploration of the parameter space within the STEPS method is the subject of future work.
Slip rate distributions across the 12 hypothetical test cases remain similar when sample size increases. Figures 5a and 5b show the probability density functions of slip rates calculated from 600 successful earthquake path simulations, and Figures 5c and 5d shows the same result, but from 1500 successful earthquake path simulations. Based on this sensitivity analysis, we target 600 successful simulations in our subsequent analysis to optimize computing time while maintaining model fidelity.

Data Selection
We present results from two example faults, the high strain rate San Andreas fault (California, USA) and the low strain rate Toe Jam Hill fault (Washington, USA). Both faults have prior studies in which at least two displacement-time observations have been made at the same site, in addition to observations of the most recent event.

San Andreas Fault Data
The San Andreas fault test case relies on field measurements made by Sieh and Jahns (1984) at the Wallace Creek site (Figure 6a(i)). At this locality, two displacement-time (DT) observations were made, each yielding a mean slip rate averaged from that observation to the origin: a channel was displaced 128 ± 1 m since 3.68 ± 0.155 kyr B.P. (DT1 on Figure 6a(i)), yielding a slip rate of 33.9 ± 2.9 mm/yr, and an older channel was displaced 475 ± 3 m (DT2 on Figure 6a(i)) since 13.25 ± 1.65 kyr B.P., yielding a slip rate of 35.8 + 5.4/−4.1 mm/yr (Figure 6a; values as reported by Sieh & Jahns, 1984). Given these results, the mean slip rate at Wallace Creek is reported as consistent over millennial time scales. The timing of the most recent event on this section of the San Andreas, the 1857 Fort Tejon event, is constrained by historical records and interpretation of high-resolution digital topography  (MRE on Figure 6a(i)). Single event displacement and recurrence interval distribution constraints for Wallace Creek are derived from field studies at the Bidart Fan, <10 km to the southeast. Ludwig et al. (2010) showed 15.9 m of offset accrued in the last five earthquakes at the nearby Bidart Fan site, producing a mean displacement of 3.2 m. Using this mean value, we construct a Gaussian distribution with a mean of 3.2 m and standard deviation of 1.9 m (CV of SED distribution ∼0.6; Hecker et al., 2013;Lin et al., 2020;Nicol et al., 2016) evaluated over 0.1 and 10 m. Truncating the distribution to prevent negative SED values induces a right-hand skew to the distribution while honoring the field observations incorporated in the pre-truncated parameterization of the SED distribution. Therefore, the SED sampling distribution mean and standard deviation are adjusted prior to truncation to maintain a mean of 3.2 m in the SED sampling distribution after truncation (CV of final distribution remains ∼0.6). The mean recurrence interval used to generate the exponential sampling distribution (CV of RI = 1) is 88 years , also measured at the Bidart Fan site, and we institute a minimum recurrence interval of 15 years. We assume that the Bidart Fan and Wallace Creek sites experience the same number of coseismic ruptures and very similar displacements in those ruptures, given that there is no obvious structural complexity between the two sites along the mature San Andreas fault. A lower CV of RI (CV = 0.7) is explored in the Discussion section of this manuscript.

Toe Jam Hill Fault Data
The Toe Jam Hill fault test case relies on field measurements presented by Nelson et al. (2003) along the length of this fault, a north-dipping backthrust of the Seattle fault (Figure 6b(i)). Specific event ages and per-event displacements utilized in this study were measured by Nelson et al. (2003) at the Crane Lake site. At this site, a paleoseismic trench exposure provides evidence for three, and possibly four, earthquakes between 1 and 2.5 ka (MRE and DT1 on Figure 6b(i)). These events are noted at additional trench sites along the fault; however, per event displacement is best preserved at the Crane Lake site of ∼1.5 m per event. In addition to this late Holocene cluster of earthquakes observed along the Toe Jam Hill fault, Nelson et al. (2003) also documented an older event at their Saddle trench with 0.3 ± 0.1 m of fault slip dated at >16 ka, coupled with previously observed 8-14 m offset of glacial stratigraphy across the fault zone at Alki Point, ∼6 km east of the Crane Lake trench site (DT2 on Figure 6b). Though the paleoearthquake record may or may not be complete, very little displacement has accumulated between DT1 and DT2, implying either numerous small displacement events or a small amount of moderate displacement events. The mean RI used to generate the exponential distribution is 0.8 kyr (Nelson et al., 2003), and we aim to match field-measured mean and standard deviation of SED of 1.5 nd 1 m, respectively (Nelson et al., 2003) (Figure 6b) but limit the SED range to 10 cm and 5 m. Just as in the San Andreas example, we fit these constraints by creating a truncated normal distribution with an adjusted mean and CV to avoid biasing the sampling distribution after truncation. A minimum recurrence interval was implemented to truncate the RI sampling distribution at a minimum value of 50 years. The Toe Jam Hill fault is a part of the complex Seattle fault system; thus, earthquakes recorded on the Toe Jam Hill fault may not be representative of the entire fault system. However, the Toe Jam Hill fault presents a well-dated record of earthquake occurrence along a slow fault with documented changes in slip rate behavior, which makes for a useful test case to test the STEPS method.

Results
In the following section, we describe results of the STEPS methodology using earthquake paths plotted through time, density plots of earthquakes through time, and slip rate probability plots.

San Andreas Fault Test Case
The earthquake paths that successfully fit the field-based DT observations from the Wallace Creek site are plotted in Figure 7a. The DT observations, with the swath connecting the edges of the boxes, is shown for reference. We highlight two randomly selected earthquake paths from the set of paths to easily view their shapes. Both highlighted paths are representative of the total population of earthquake paths in that these paths contain both earthquake lulls (long interevent times) and earthquake clusters (periods with short interevent times).
The mean and CV calculated for each path independently (Figure 7a) are similar because they both result from independent samples of the same parent distributions (i.e., Figure 7a inset histograms; 600 successful paths from 3,033 total attempted paths). Despite the similarity of the CV of two highlighted paths of Figure 7a, in detail, each path is different along their length. For example, at c. 11-12 ka, the highlighted paths diverge with ∼50 m between the paths. The reason for this divergence is these two highlighted paths host earthquake clusters and lulls that are offset in time resulting from the random ordering of SED and RI picks (e.g., Cowie et al., 2012). This divergence of paths with the same CV highlights the importance of random ordering on SED and RI picks from reasonable SED and RI distributions. Assuming independence of displacement and time yields this random ordering, which could be refined in future iterations of STEPS.
We utilize a density plot to spatially analyze all earthquake paths simultaneously (Figure 7b). The density plot shows the counts of the earthquakes within bins specified here at 100 years and 1 m intervals. Between the most recent event (MRE) and DT observation box at c. 4 ka, the densest overlap of the earthquake paths  (Sieh & Jahns, 1984) and STEPS-derived Kernel Density Estimations (KDEs). Inset table reports values at specific quantiles for 5*mean RI (0.44 kyr) and 10*mean RI (0.88 kyr). occurs at the observation box locations, as expected from the model requirements. As such, the densest parts of the earthquake path distribution extend beyond the swath between these younger two boxes. Following this box, until c. 5 ka, the earthquake paths remain densely clustered, due to the presence of the DT observation box. The density plots highlight that where DT boxes are tightly constrained, the earthquake paths around the displacement-time box are also limited. Compare, for example, the tightly constrained density plot near the 4 ka DT box, which has an uncertainty of ±1 m, with the wide density plot prior to the older (13 ka) DT box. As shown in Figures 3 and 4, however, the slip rate variability defined in the model parameters is recovered in the time periods between the two displacement-time boxes.
To estimate slip rates, all paths that obey displacement-time constraints were discretized at two different time intervals representing 5x and 10x the mean RI (e.g., 0.44 and 0.88 kyr) (Figure 6c). The Kernel Density Estimations (KDEs) of these slip rates from the Wallace Creek example are shown in Figure 7c, plotted with the uniform distributions of rates reported in the literature by Sieh and Jahns (1984).
The mean geologic slip rate and STEPS-derived estimates of slip rate variability have similar median values, but the STEPS method results in broader uncertainty estimates reflecting the earthquake-to-earthquake variability (Figure 6c). Assuming the Sieh and Jahns (1984) reported rates represent the 68% confidence interval (CI), the literature reported rate ranges of 31-36.8 mm/yr (short-term record) and 31.7-41.2 mm/yr (long-term record) are similar to, but are smaller than, the STEPS-derived 68% CI ranges of 17.9-55.2 mm/ yr over 5*mean RI.
The KDE of STEPS slip rates tightens slightly when averaging slip rates over 10*mean RI (0.88 kyr) instead of to 5*mean RI (0.44 kyr). The peak of the 10*RI distribution is taller, the median value is shifted slightly higher, and the 68% and 95% confidence interval bounds are narrower. These observations indicate that slip rates are more tightly clustered about the median value when averaged over more earthquake cycles, as expected. Despite the tightening of the distributions when averaging over 10*mean RI vs 5*mean RI, the distributions remain broad and encompass a wide range of possible, though relatively unlikely, slip rates observed only when interpolating between the displacement-time observation boxes with earthquake paths as opposed to straight-line approximations. Both the 5 and 10*mean RI STEPS rate distributions have a slight right-hand skew, perhaps due to the random ordering of SED and RI picks from the sampling distribution. This is further explored in the Discussion section.
To compare the STEPS over longer term, we average over 50 and 100 mean RI (4.4 and 8.8 kyr) intervals. The output slip rate distributions become narrower at increasing discretization intervals as the 50 and 100*mean RI intervals of the earthquake paths become more similar, and therefore approximate the means of the input sampling distributions. This decreasing in width of STEPS-derived slip rate distributions with increasing discretization interval highlights the effect of either including (5-10*mean RI) or excluding (50-100*mean RI) the effects of variable RI and SED on slip rate (i.e., the possibility of sampling only an earthquake cluster or lull is increasing when calculating slip rate over smaller time intervals vs longer time intervals).

Toe Jam Hill Fault Test Case
Results for the Toe Jam Hill study are presented in Figure 8. In addition to the post-glacial record constrained by Nelson et al. (2003), we plot simulation results of 30 additional older events, extending the record to c. 40-50 ka (Figure 8b). These events are randomly sampled from the predefined sampling distributions of SED and RI (Figures 6b, 8a, and 8b inset histograms). The Toe Jam Hill record is extrapolated for this additional 30 events, after fitting the field observations, to average slip rates over a longer record of slip (Figures 1c, 8b, 8d, and 8f). As with the Wallace Creek example, we plot the density maps (Figures 8c  and 8d) and slip rate distributions (Figures 8e and 8f) for the STEPS results at Toe Jam Hill for the post-glacial record and extrapolated record, respectively. The goal in extending the model record to c. 40-50 ka is to explore the behavior of earthquake paths in the absence of long-term geologic constraints on a low strain rate fault.
Focusing first on the post-glacial record (post c. 16 ka) in Figure 8a, we observe a broadly distributed population of earthquake paths that satisfy the field-observed data (600 successful paths from 323,952 total attempted paths). The Toe Jam Hill fault requires many more attempts to fit the displacement-time observations, implying that there may be a change in slip per event and/or recurrence interval between the younger, faster portion of the record compared to the older, slower portion of the record, as was reported by Nelson et al. (2003). For the purposes of this study to contain all available DT markers in a single slip rate uncertainty distribution, we continue with this analysis while noting this marked change in fault behavior between the MRE and DT1 vs. DT1 and DT2. Two selected paths are highlighted in this post-glacial portion to illustrate key features (Figure 8a). Both highlighted paths have similar mean recurrence intervals (RI) (c. 1.2-1.4 kyr) and similar single event displacements (SED) (∼1.2 m) along their respective lengths. Despite these similarities, these two earthquake paths do not match in number of events nor absolute location in the areas between the displacement-time observation boxes. These differences in earthquake path shape highlight that parameters typically used to describe earthquake behavior through time may satisfy a basic understanding of earthquake frequency and size, but do not provide a unique solution for event-by-event behavior. We do not impose event-by-event behavior from paleoseismology chronologies on the earthquake paths construction here to allow for events where the paleoseismic record may be incomplete or imprecise. The successful picks of SED and RI of the Toe Jam Hill test case (Figure 8a inset) show that the short-term mean RI and SED do not match the earthquake paths required by the DT boxes. The SED is preferentially selected from the smaller displacement portion of the of the sampling distribution data (mean value of output = 1.1 m vs. mean value of input distribution = 1.5 m). Likewise, the long right-hand tail of the exponential RI distribution is sampled more frequently than what was imposed with the input data (mean value of output = 1.63 kyr vs. median value of input = 0.8 kyr). These observations suggest that the STEPS method may have a role to play in determining better refined distributions of SED and RI, especially with the addition of paleoearthquake constraints and nuanced distributions applied to the displacement-time boxes. This parameter exploration remains a topic of future investigations.
The density map of the post-glacial record (Figure 8c) highlights that the majority of earthquake paths travel through and above the straight lines (swath) connecting the displacement-time observation boxes. The skew of earthquake paths to favor cumulative displacements greater than that predicted because few aspect ratios of paths can exist given the requisite low slip rate between DT boxes 1 and 2. As was observed in the Wallace Creek example, the density of earthquake paths is greatest, with the distribution of earthquake paths at its narrowest, near the observation boxes due to the model requirements.
Like the Wallace Creek example, slip rates are discretized by 5 and 10 times the mean recurrence interval; in the Toe Jam Hill example, this yields intervals of 4 and 8 kyr, respectively (Figure 8e). These slip rate distributions are plotted as KDEs, with the literature-reported rates of Nelson et al. (2003) slip rates plotted as uniform distributions. Here, we observe that the STEPS-derived slip rate distributions are relatively broad, and overlap both the slower, long-term and the faster, short-term slip rates of Nelson et al. (2003) at the 68% confidence interval (CI). For example, the 68% CI is 0.23-1.33 mm/yr for the 4 kyr STEPS rate. The 68% CI is slightly narrower when slip rates are averaged over a longer period of time (8 kyr), but generally represent similar ranges of possible slip rates. Both cases of STEPS rates have a right-hand skew in their distributions.
In Figure 8b, the same post-glacial record is shown, with the addition of the earthquake paths traversing the additional 30 events. Here, we observe that, without a displacement-time observation box constraining the older end of these extrapolated earthquake paths, the dispersion of the earthquake paths is much greater compared to the well-constrained post-glacial record. Additionally, the successful SED and RI picks mirror the input sampling distributions in absence of another displacement-time observation box to further filter the successful boxes (Figure 8b inset histograms). The density map (Figure 8d) confirms this broad scatter of earthquake paths. Despite the broad area of earthquake paths, there is a distribution of earthquake paths, with the densest section of paths occurring near the middle of the path distribution, again highlighting that the output values reflect the input sampling distributions of SED and RI. The slip rate distributions (Figure 7f) are much flatter and have more of a right-hand skew to the distribution in comparison to the post-glacial record (Figure 7e) as expected given the absence of constraints for the oldest 30 events. The extrapolation of earthquake paths highlights the potential prevalence of earthquake clustering when sampling over a longer time interval than what is available in the post-glacial record alone along a low strain rate fault such as the Toe Jam Hill fault. The wide range of the slip rate paths is a result of the combination of SED and RI when unfiltered by a DT box. However, the median values of slip rate compared between the post-glacial and extrapolated examples are largely similar and reflect the influence of the input SED and RI sampling distributions on the output slip rate distributions. In the former case, the mean (and CV) of RI and SED are 1.68 kyr (1.1) and 1.21 m (0.7), compared to the latter case of 1.02 kyr (1.1) and 1.53 m (0.6). In both cases, the mean RI and SED is greater than the mean of the RI and SED sampling distributions (0.8 kyr and 1.5 m), suggesting that, given the model setup, a longer RI and smaller SED is required to satisfy the earthquake cluster-lull pattern observed by Nelson et al. (2003).
In both the post-glacial and extrapolated earthquake path sets, we utilize the STEPS approach to calculate the slip rate from the initial to final point along each earthquake path to simulate a long-term geologic slip rate. Unlike the Wallace Creek example, the Toe Jam Hill fault record cannot be reliably averaged over 50 or 100*mean RI, given that the Toe Jam Hill fault is a much slower fault and records fewer earthquakes over the same time scale as the San Andreas fault. The purpose of this analysis is to show the relationship between slip rate distribution size and the number of mean recurrence intervals used in the discretization. By averaging over more events, STEPS slip rate distribution becomes tighter. This is particularly well expressed in the post-glacial STEPS-derived slip rate distributions (Figure 8e). Here, the end-to-end (longest possible term, c. 14 kyr) slip rate distribution is very tight and narrow, representing a very likely slip rate of ∼0.7 mm/yr. In the extrapolated record, the longest-term STEPS-derived slip rate has a slightly broader distribution, reflecting the scatter of earthquake paths in the extrapolated record without the additional displacement-time box to further constrain the record (Figure 8f). In either case, these end-to-end slip rate distributions highlight the long-term slip rates, while the shorter-term STEPS-derived slip rates do not. The STEPS-derived slip rates averaged over 5 and 10 mean recurrence intervals represent a short-term distribution of slip rate that subsamples this longer-term history, which is recoverable by averaging over the entire earthquake path.

Discussion
In this discussion section, we compare our results to previous studies, explore statistical tests to quantify the significance of the synthetic slip rate distributions, provide best practices for field geologists, discuss some limitations of the STEPS approach, and highlight potential applications of this methodology in seismic hazard analysis.

Statistical Comparisons of Slip Rate Estimates
Following Styron (2019), we employ the two sample Kolmogorov-Smirnov (KS) test to identify potential differences between the STEPS-derived slip rate distributions. The KS test is an appropriate tool to assess the goodness of fit of different distributions by quantifying the statistical difference or similarity across two distributions (e.g., Massey, 1951). This non-parametric test provides a quantification of the absolute distance between two cumulative distribution functions (CDF). By definition, the two sample KS test compares the aspect ratio of the cumulative distribution functions. The null hypothesis of this test is that the two distributions tested could be drawn from the same distribution at some threshold of significance. The null hypothesis is rejected if the maximum absolute distance between the CDFs is greater than a critical value, Dcrit. Dcrit is dependent on the significance level testing and is inversely proportional to number of positions over which the CDF-to-CDF distance is calculated. These critical values are based on tabulated values at specific significance levels and sample sizes. If the null hypothesis is rejected, the two distributions are considered statistically different from one another.
We first compare the two STEPS-derived slip rate CDFs from the San Andreas fault example (Figure 9a). We plot the CDFs of the 0.44 kyr (5*mean RI) and 0.88 kyr (10*meanRI), discretized at 200 increments between 0 and 100 mm/yr. At each of these increments, we take the absolute difference between the CDFs (Figure 9b). In the San Andreas example, the difference between the CDFs never exceeds the critical value at the 5% significance level (Dcrit = 0.136) meaning the null hypothesis cannot be rejected and the two CDFs could be sampled from the same parent distribution with the highest confidence possible with the two sample KS test. This indicates that, although changes between the two STEPS-derived slip rate distributions are slight for the San Andreas example, there appears to be no statistical difference between results using the two discretization values.
We next compare the STEPS-derived slip rate CDFs for the Toe Jam Hill fault example (Figures 10a-10c(i)). In Figure 10a, we compare the slip rate distributions for the post-glacial portion of the record, constrained on both ends by displacement-time observation boxes. Here, we see again that the null hypothesis cannot be rejected, indicating that the 4 kyr (5*mean RI) and 8 kyr (10*mean RI) could be derived from statistically similar distributions at the 5% significance level (Dcrit = 0.136) (Figure 10a(ii)). Similarly, when comparing the slip rate distributions including the extrapolated record beyond the post-glacial record, we cannot reject the null hypothesis and must accept that the two distributions are statistically similar at the 5% significance test results. (a) Cumulative density functions (CDFs) for 5*mean RI and 10*mean RI slip rates. Distance between two CDFS shown with vertical gray lines; (b) Absolute distances from panel i plotted against slip rate, with color of marker indicating which CDF is greater/has a higher probability at a given slip rate. Lack of markers above the Dcrit marker indicates the null hypothesis of the KS test cannot be rejected. level (Figure 10b(ii)). However, in comparing the CDFs between Figures 10a and 10b, we observe that the CDFs including the extrapolated record on Figure 8b are far broader with respect to the post-glacial record in Figure 10a. To formalize this comparison, we perform the KS test between the 5*mean RI slip rate distributions for the post-glacial and for the extrapolated records in Figure 7c. Here, we observe that the differences between the CDFs exceeds Dcrit at the plotted 5% significance level, as well as at the 1% significance level (Dcrit = 0.163), leading to rejection of the null hypothesis (Figure 9c). These distributions of STEPS-derived slip rates are therefore statistically different.
The difference between the CDFs of the post-glacial record and of the extrapolated record is due to the lack of displacement-time observation constraining the 'wander' of the earthquake paths further back in time. This may manifest as an overrepresentation of slow slip rates and fewer earthquakes in the post-glacial period compared to the possibility of faster slip rates and potentially more earthquake clusters previous to the 16 ka-present record. Without an older displacement-time observation box, there is no ability to filter out the earthquake paths that may not be representative of reality in this extrapolated record. However, this extrapolation process is not solved satisfactorily in this current iteration of STEPS and will be further explored in future work, perhaps in a full Bayesian framework.

Implications for Geologists
The model-derived slip rates from STEPS offer a more complete picture of shorter-term geologic slip rate variability than captured by mean geologic slip rate, but the numerical model can only exist with field-based slip rate studies to constrain the simulations. Field-based data are the basis of this model; compiling as many measurements of dated offsets at a site or reach of fault as possible will help calibrate and constrain numerical simulations such as these.
Determining the paleoearthquake chronology and single-event displacements along a fault is also critical to the STEPS approach, because these observations can inform the distribution of interevent times and coseismic displacements. These values are of clear importance for determining the sampling distributions used to set up the STEPS problem (Figures 3, 4, 8a). However, even with these sampling distributions informed by field measurements, the ordering of events is not necessarily predicted by CV of sampling distributions, nor of individual paths, alone (Figures 7a and 8a). Previous compilations of single event displacement suggest that single-event displacement is less variable than recurrence interval with a global CV of ∼0.5-0.6 (e.g.,  Figure 10. (a) 5*mean RI vs 10*mean RI for the post-glacial record only. (b) 5*mean RI vs 10*mean RI for the post-glacial record and extrapolated record combined. (c) 5*mean RI post-glacial record vs. 5*mean RI post-glacial and extrapolated record. Lack of markers above the Dcrit marker indicates the null hypothesis of the KS test cannot be rejected (panels a(ii) and b(ii)); markers exceeding the Dcrit marker reject the null hypothesis (panel c(ii)). Hecker et al., 2013;Lin et al., 2020;Nicol et al., 2016), but understanding these relationships through time can prove challenging with poorer and less certain recovery of per-event displacement for older and older events. With more detailed and longer records of both recurrence and single-event displacement collected along faults and within fault systems, a potential relationship may emerge to aid in understanding the true independence of recurrence and single-event displacement through time (e.g., Shimazaki & Nakata, 1980;Weldon et al., 2004).
Utilizing parameters from high-resolution, long-term paleoseismic chronologies that may yield a non-Poissonian recurrence (CV ≠ 1) embeds the assumption that the recent behavior is representative of the longterm record. This might be true on mature or fast-moving faults. For example, data along the San Andreas fault provide good bounds on recurrence interval mean, standard deviation, and coefficient of variation, suggesting a CV near 0.7 over 30 events (e.g., Scharer et al., 2010) and similar CV values have been found worldwide for most fast-moving faults (Yuan et al., 2018). However, as Salditch et al. (2020) found utilizing numerical simulations, the CV of recurrence interval may not be stable throughout a multi-millennial record of slip. Understanding CV of recurrence beyond the past few millennia is challenging with field methodologies alone. Although these records are typically hard to come by in field studies, particularly for low strain rate faults, more records of both short-and long-term earthquake behavior, especially those spanning clusters and lulls, may improve our understanding of the temporal consistency of CV for recurrence interval. Adding displacement-time markers along faults throughout the landscape evolution along a given fault, where possible, especially within fault systems, may lead to a deeper understanding of slip rate variability through time within an interrelated zone or region of faulting (e.g., Wallace, 1987). More complex modifications to STEPS could implement likelihood gains when similar displacement-time ratios have been documented over multiple time periods, as this implies slip rates are fairly stable.
We present an example of employing a recurrence distribution rooted in paleoseismic (short-term) observations in the STEPS (long-term) modeling approach utilizing a lognormal RI sampling distribution for the Wallace Creek example (Figure 11). Here, we observe that, in comparison to Figures 6 and 7, which utilized an exponential RI distribution, applying the lognormal sampling RI distribution (Figure 11a) limits the width of successful earthquake paths (Figure 11b). Notably, the earthquake paths in the lognormal framework favor the older end of the oldest displacement-time observation due to the change in shape of the RI sampling distribution. The lognormal RI distribution has less density at the shortest RI in comparison to the exponential distribution, and therefore should skew the selected RI to be slightly higher and less random. Despite this tighter distribution of earthquake paths and apparent preference to utilize the older portion of box DT2, the resultant slip rates are similar to selecting the exponential sampling distribution when sampled at 5*mean RI (0.44 kyr) (Figure 11d). This preference to sample the older portion of box DT2 reduces the upper limit of the slip rate distribution, reducing the 97.fifth quantile slip by about 10% compared to the exponential case (Figure 7c vs. 11d). Additionally, the 10*mean RI (0.88 kyr) STEPS rates produces a reduced right-hand skew, with the median value more centrally located within the distribution compared to the exponential case (Figure 7c vs. 11d). The median STEPS rates are also slightly slower in the lognormal case compared to the exponential case. Overall, the successful earthquake paths have a mean RI of 0.11 kyr (CV = 0.8) and mean SED of 3.75 m (CV = 0.6). The output values here favor a longer recurrence interval, consistent with the observation of sampling the older end of box DT2, and a larger and more periodic single event displacement, suggesting that a longer recurrence interval needs to be balanced by a greater single event displacement in order to satisfy the DT box constraints. This example shows that, if the assumption is made that the recent record is emblematic of the entire Holocene-late Pleistocene history, further refinements to the earthquake path distributions and subsequent slip rate estimates can be made.
An additional refinement can be made to this example by adding a paleoearthquake chronology to the suite of displacement-time constraints. In Figure 11c, we add the paleoearthquake ages of events that have paired displacements at the Bidart Fan site Ludwig et al., 2010). The paleoearthquake chronology at Bidart Fan is a well-resolved record of the penultimate five events following the MRE, but is this recent record a unique indicator of the longer earthquake record? Can the paleoearthquake record help limit the number of earthquake paths by providing an independent constraint on the slip history at this site? Akçiz et al. (2010) established evidence of five ground rupturing earthquakes since ca. 1400 C.E. and Ludwig et al. (2010) showed that 10 m of slip accrued in the last two earthquakes, whereas the previous three earthquakes accrued a total of 5.9 m of slip across a nearby channel. For the purposes of this analysis, we combine displacement measurements reported by Ludwig et al. (2010) for the penultimate event (1631-1823 C.E.) and fifth event back (1450-1475 C.E.) reported by Akçiz et al. (2010).
The successful paths described in the previous section were then re-interrogated with these additional displacement-time observations. Using these additional constraints, an earthquake path is considered successful if it satisfies: (a) all Wallace Creek constraints (Figure 7a), and (b) both additional paleoearthquake  Figure 11c also satisfies the observations of the Bidart Fan paleoearthquake chronology (∼180 earthquake paths). The earthquake paths that fit the Bidart Fan paleoearthquake record limit the area over which earthquake paths travel in the most recent 600 years (Figure 10c inset). After this point in the slip history, though, the Bidart Fan successful earthquake paths appear to be randomly distributed within the general population of earthquake paths. The resultant slip rate distributions here very closely resemble the slip rate distributions utilizing the 600 earthquake paths with a lognormal RI sampling distribution (Figure 11e vs 11f).
These additional analyses show that refined sampling distribution assumptions (e.g., going beyond the Poissonian assumption of earthquake recurrence) may narrow the resultant STEPS slip rate distributions. Additional closely spaced displacement-time observations, such as the Bidart Fan paleoearthquakes, though, only provide reduction of earthquake paths in close proximity to the boxes. This relationship was also observed in the post-glacial record of Toe Jam Hill fault with the closely spaced displacement time observations ( Figure 8).
An especially useful observation is the timing and displacement of the most recent earthquake (MRE), as this helps mitigate the origin/open interval problem (Figure 1b). In general, estimating the distribution of permissible slip rates with the STEPS approach benefits from collecting as much geologic data (geomorphic offsets and ages, paleoearthquake ages, and displacements) as possible.
Observations important for simulations such as ours, but that are typically not reported in geologic slip rate studies, are constraints on the "null space." We define the "null space" as the region within a displacement-time plot that earthquake paths are unlikely to occupy, based on any reasonable observations (geomorphologic, geochronologic, paleoseismic, hydrologic, climatic, or similar) (e.g., a 200-ka geomorphic surface with a single event fault scarp). Identification of the null space may not be encompassed in displacement-time measurements alone, but in understanding the landscape and fault zone as a whole. This effort may allow for constraining allowable displacement-time space within which earthquake paths can reasonably travel, and perhaps trimming of the extent to which the earthquake paths deviate from the straight-line connections between boxes.
The STEPS method may provide a way to narrow and refine understandings of SED and RI distributions, particularly through time. Salditch et al. (2020) suggested that mean and CV of RI may not be constant through time. This is perhaps true in the Toe Jam Hill fault example explored in Section 3.2 of this study. While condensing disparate observations to a single estimate of slip rate and uncertainty is useful for seismic hazard modeling, STEPS may provide an opportunity to untangle nuances in slip rate variability through time. The STEPS method can help test hypotheses about changing parameters (e.g., RI and SED) through time, coupled with implementation of earthquake cycle assumptions and physics, especially within fault systems.

Comparisons to Prior Work
By focusing on a suite of paths derived from individual earthquakes occurring between geological constraints, the STEPS approach is fundamentally different than methodologies that quantify slip rates by focusing primarily on offset features. This work builds from Gold and Cowgill (2011) and subsequent applications of their method (e.g., Burgette et al., 2020;Zinke et al., 2017Zinke et al., , 2019, which utilized straight line constant slope approximations between displacement-time observations. The Cumulative Offset-Based Bayesian Recurrence Analysis (COBBRA) method by Fitzenz et al. (2010) uses a Bayesian approach to refine RI probability distributions based on displacement-time estimates. They apply this approach to the Dead Sea fault to test if time-dependent recurrence is preferred over a random distribution. In this sense, the COBBRA method is focused on recurrence parameters, whereas STEPS was formulated to understand slip rate uncertainty. Similar to STEPS and COBBRA, Styron (2019) also utilized summed step functions to represent earthquake behavior through time. The STEPS model currently utilizes a graphical, forward-modeling approach compared to the Bayesian approach of Fitzenz et al. (2010) and leverages the Styron (2019) approach with the addition of input data based on field-based displacement-time observations. The STEPS results support the finding of Styron (2019) that averaging over a minimum of five earthquake cycles may be sufficient to initially minimize variability in slip rate, with a 10-20 earthquake cycle recommended to reach a stable minimum of slip rate uncertainty (Figure 6c). In addition to the minimum of five mean recurrence intervals to obtain a stable slip rate, the STEPS approach suggests that increasing the discretization of slip rate calculation from five to ten mean recurrence intervals continues to refine the probability distribution. Despite this incremental, visual change in the distributions observed in Figures 7c,  8e, and 11e, we do not observe a statistical difference between these distributions in either the San Andreas or Toe Jam Hill fault examples in the KS test outcome (Figures 9 and 10). This relative consistency between results given variable slip rate sampling interval may be dependent on recurrence interval sampling distribution CV; a more clustered CV (CV >> 1) may yield significantly different results across the 5 and 10*mean RI slip rate intervals.  (2019) also find that the distribution shapes of SED and RI have a modest effect on the simulated slip rate. Our analysis, for example, shows modest changes to slip rate uncertainty when a boxcar versus random distributions for recurrence interval is used with small displacement-time observation boxes (e.g., Figure 5a). Additionally, when exploring the influence of CV and recurrence interval sampling distribution shape within the San Andreas fault example, we observe a small reduction in the width of the 68 and 95% confidence intervals of slip rates measured at both the 5 and 10*mean RI intervals in the lognormal RI case compared to the exponential RI case (Figure 11d). Regardless of sampling distribution choices explored here, the calculated CV of recurrence interval or single event displacement for each path is slightly different than the CV of the sampling distribution, because values for each path are sampled independently and without replacement. Fitzenz et al. (2010) suggested that, simply because no distribution of recurrence performed significantly better in COBBRA, this does not mean that they all yield the same behavior, and hazard models should instead utilize a combination of recurrence assumptions. The influence of sampling distribution choices, as well as a possible combination/logic tree approach of recurrence models, on resultant slip rate distributions could provide additional questions for applications of the STEPS approach.
Like Styron (2019), we also observe a right skew in the slip rate distributions of the San Andreas and Toe Jam Hill fault examples. In the STEPS approach, this right-hand skew may be due to the overrepresentation of events in an earthquake cluster given the discretization method applied. Averaging over increasingly large time intervals (e.g., >20*mean RI) would perhaps reduce this right-hand skew in the slip rate distributions. Further reducing the CV of the recurrence interval sampling distribution, beyond the slight reduction explored herein between CV = 1 ( Figure 7) and CV = 0.7 (Figure 11) in the San Andreas fault example, may continue to minimize this right skew. By the same token, increasing the CV of the recurrence interval sampling distribution reflecting clustered earthquake behavior may extend the right-hand skew and/or result in a bimodal slip rate distribution. Additionally, introduction of earthquake cycle assumptions (i.e., slip, time, or strain predictability) may reduce (or exacerbate) this right-hand skew in output slip rate distributions.

Limitations of the Methodology
As with all modeling approaches, the STEPS method has limitations. One set of limitations is that the utility of the model is dependent on the quality of the field-based observations. The displacement-time observations (dated offsets) used to constrain earthquake paths must be accurate and precise. For the construction of earthquake paths, the estimates of mean recurrence and single event displacements, as well as their distributions, must be reasonable. This is perhaps the hardest data to collect because very few faults have paleoearthquake records sufficiently detailed to confidently characterize the distributions of offset and interevent times. Finally, a displacement-time constraint on the most recent event is required by the model as currently structured, to avoid averaging the earthquake paths through the open interval (Figure 1b).
Poorly studied fault zones that are thought the host earthquake clusters and likely have earthquake recurrence intervals with CV >>1, such as the New Madrid fault zone, may make a poor candidate for the STEPS methodology, given the lack of understanding of mean recurrence interval and the distribution of recurrence interval.
The current formulation of the STEPS approach does not impose an earthquake renewal model, magnitude-frequency distribution, rheology, or any parameterization of the earthquake cycle beyond displacement and recurrence. Therefore, the earthquake paths are created only by random selection and ordering of SED and RI picks, and not governed by any evolution of the earthquake paths themselves. By incorporating variability in RI and SED, though, we implicitly incorporate reasonable surficial observations of these higher-level components of earthquake recurrence. Once these factors are incorporated into STEPS, we can begin to explore the influence of these underlying behaviors on slip rate uncertainty in terms of the sequence of these SED and RI picks. We would expect that implementation of other constraints would likely impact the ordering of RI and SED picks. Implementation of additional physics and renewal assumptions will be the topic of future exploration.
Finally, the STEPS approach only provides slip rate uncertainty at a given site and not along the entire length of a given fault. Geologic slip rates are typically extrapolated along sections of fault; importantly, slip rate uncertainty is therefore only valid at a specific location given that the displacement-time observation constraints come from a specific location.

Implications for Hazard Analysis
The STEPS approach presented here may prove useful for integration into seismic hazard analysis. First, geologic slip rates are critical inputs in determining the activity and recurrence of a given fault in seismic hazard modeling. Geologic slip rates are often used as either benchmarks or direct inputs into crustal deformation models. A geologic deformation model may be produced in conjunction with geodetic deformation models, which are based primarily on GPS data. In the Uniform California Earthquake Rupture Forecast, version 3 (UCERF3), four deformation models were implemented with varying usage of geologic input data (Parsons et al., 2013), which incorporated per-earthquake uncertainty of long-term geologic slip rates. Additionally, two geodetic-based deformation models were implemented across the western U.S. in the U.S. National Seismic Hazard Model (NSHM) of 2014 (Petersen et al., 2015). Yet, geologic slip rate uncertainties were not explicitly incorporated outside of California into the 2014 or 2018 versions of the NSHM. The incorporation of multiple branches of the slip rate logic tree provides uncertainty in the determination of the slip rates, but does not provide a quantification of, for example, slip rate variability through time. Work by Zeng (2018) provides a path forward for understanding different distribution shapes using existing minima and maxima of geologic slip rate uncertainty, but again does not capture the uncertainty associated with slip rates that may be expressed beyond reported measurement uncertainty.
The STEPS method provides a way to combine multiple and disparate observations of geologic slip rate into a single estimate of slip rate uncertainty. For example, the U.S. NSHM is currently a time-independent calculation and utilizes a single value of slip rate with associated uncertainty. Therefore, this does not allow for incorporation of multiple, temporally variable slip rates. Incorporation of multiple, disparate geologic slip rates in a single estimate of slip rate uncertainty was explored in Section 3.2 of this manuscript with the Toe Jam Hill example. In this example, while there could easily be two distinct distributions of SED and RI between the younger, faster slip rate time period compared to the older, slower slip rate time period, the STEPS method provides an estimate of uncertainty in geologic slip rate that encompasses both of these histories. We intend this initial iteration of STEPS to begin a process of sensitivity testing within seismic hazard analyses to understand how, when, and over what time intervals to potentially implement these expanded slip rate uncertainties, The STEPS approach can result in slip rate uncertainties that are significantly broader than those derived from the ages and offsets of geologic features alone, depending on the length of time intervals over which STEPS-derived slip rates are calculated. How should the broader STEPS-derived uncertainties be interpreted? Geologically derived slip rates and STEPS-derived slip rates can be nearly identical, both in terms of median and range, if the STEPS-derived slip rates are calculated over sufficiently wide time intervals, regardless of fault strain rate and model sampling distributions (Figures 7c, 8e, and 11e). In the San Andreas fault test case, for example, the STEPS-derived slip rate distributions match the median values of Sieh and Jahns (1984), but the width of the 68 and 95% confidence intervals (CI) depend on how the path is sampled (e.g., 50 and 100*mean RI slip rate distributions on Figure 7c more closely approximate the Sieh and Jahns (1984) values).
The expanded STEPS-derived slip rate distribution calculated over geologically shorter time intervals (e.g., 5*mean RI and 10*mean RI) allows for the possibility that the slip rate appropriate at probabilities of interest for seismic hazard (e.g., 2% in 50 years) could be significantly faster or slower than the median, long-term purely geologically derived value. These model-derived slip rate distributions account for the possibility of earthquake bursts and lulls, which would skew a slip rate faster or slower than the narrowly defined uncertainty bounds inferred from long-term geologic offsets. That is, the STEPS-derived slip rates capture the possibility of a given fault being 'in cluster' or 'out of cluster' (e.g., Salditch et al., 2020). Such a possibility is reasonable, given our limited knowledge of earthquake behavior through time. By constructing earthquake paths from reasonable distributions of recurrence intervals and slip and focusing on the possible paths between offset geologic features, the STEPS method captures possible slip rate variation that can be lost in simple straight-line reconstructions of offset geologic features. A utility of the STEPS method is that slip rate uncertainty can be estimated over any time interval of interest and is not reliant on the time intervals explored in this manuscript (e.g., 5 and 10*mean RI). If a seismic hazard analysis period of interest is 2% in 50 years vs 10% in 50 years recurrence of a shaking threshold, STEPS can provide the appropriate estimate of slip rate uncertainty for the given hazard metric of choice.
Currently, the possibility of being in or out of cluster is explicitly accounted for in some logic trees for specific fault sources across the central and eastern United States in the 2014 and 2018 release of the NSHM (e.g., New Madrid fault zone), and indirectly accounted for in UCERF3 through Poissonian distributions of the recurrence interval. Care should be taken, however, to avoid any double counting with aleatory variability applied at the forecasting stage (which will assume a Poisson or other distribution of recurrence intervals). Integration of the STEPS method into hazard forecasting will require a careful accounting of how potentially variable faulting behavior is modeled.
One application of the STEPS approach within seismic hazard modeling is to define a geologically based estimate of the timing of the next future event. We highlight such a calculation in Figure 12. Here, we allow successful earthquake paths to continue from the most recent event (MRE) box into the future. The earthquake path must obey the open interval constraint, meaning a successful future event will only occur in the future time, future displacement quadrant (Figure 12a). The next future event is calculated by randomly selecting an additional single event displacement and recurrence interval from the values that make up the successful 600 stairs that have already satisfied the given displacement-time observation constraints and continues to assume independence between SED and RI values as a combination and compared to their previous values. We plot the results of these calculations for the San Andreas fault in Figure 12b and for the Toe Jam Hill fault in Figure 12c. In both cases, the majority of the successful future events occur close to the origin (histograms inset within Figures 12b and 12c). In particular, the San Andreas fault example mimics the exponential distribution of recurrence interval used to generate the earthquake paths Figure 12b inset). This distribution shape of next future event timing may be due to the current length of the open interval being well in excess of the applied mean recurrence interval (88 vs 164 years; Biasi & Scharer, 2019). However, in addition to the tight clustering of potential events near the origin of Figures 12b and 12c, some future events plot as far into the future as 300-500 years in the case of the San Andreas and 4,000 years in the case of Toe Jam Hill (see inset histograms). Related analysis may be helpful for deriving the probability of displacement in the next following event in probabilistic fault displacement hazard analysis.

Conclusions
We present a numerical approach called STEPS (Slip Time Earthquake Path Simulations), which uses fieldbased observations of data offsets and returns distributions of fault slip rates that encompass variability in recurrence interval and single event displacement. We discuss challenges facing geologic slip rate determination, and address potential paths forward using the STEPS approach, using a hypothetical set of test cases and two previously published slip rate sites along the Carrizo section of the San Andreas fault and the Toe Jam Hill fault within the Seattle fault zone. We find that using the STEPS approach enables a fuller understanding of permitted slip rate uncertainty compared to literature values. The STEPS approach accounts for earthquake cycle variability through time, providing, for the first time, a methodology to define slip rate uncertainty beyond dated offset measurement error.