Major California faults are smooth across multiple scales at seismogenic depth

Surface traces of earthquake faults are complex and segmented on multiple scales. At seis-mogenic depth the detailed geometry of faults and earthquake rupture is mainly constrained by earth-quake locations. Standard earthquake locations are usually too diffuse to constrain multi-scale fault geometry, while differential-timing relocation mainly improves finest scale precision. NLL-SSST-coherence, an enhanced, absolute-timing earthquake location procedure, iteratively generates traveltime corrections to improve multi-scale precision and uses waveform similarity to improve fine-scale precision. Here we apply NLL-SSST-coherence to large-earthquake sequences and background seismicity along strike-slip faults in California. Our relocated seismicity at seismogenic depth along major fault segments and around large-earthquake rupturesoftendefinessmooth,planarorarcuate,near-verticalsurfacesacrossthesub-kmto10’sofkmscales. Theseresultsshowthatmulti-scalesmoothfaultsegmentsarecharacteristicofmajor,strike-slipfaultzones andmaybeessentialtolargeearthquakerupture.Ourresultssuggestthatsmoothnessandcurvatureoffaults influencesearthquakeinitiation,rupture,rupturedirectionandarrest,andcandefineearthquakehazard.The resultscorroboratethatsurfacetracesofstrike-slipfaultzonesreflectcomplex,shallowdeformationandnot directlysimpler,mainslipsurfacesatdepth,andsupportuseofplanarorsmoothlycurvedfaultsformodeling primaryearthquakerupture. Non-technical summary Surfacetracesofearthquakefaults,likemanyfeaturesofthenaturalland-scape, are irregular and complex across many scales, from microscopic to many kilometers. However, there are conflicting views on the geometrical complexity and smoothness of earthquake faults deeper in the Earth atseismogenicdepth,wheremostearthquakesoccur-about4to15kmformanylargefaultsinCalifornia. But knowing fault geometry at seismogenic depth is key to understanding earthquakes, including the initiation, growth, and hazard of large earthquakes. The detailed geometry of earthquake faulting at seismogenic depth is mainly constrained by seismicity—earthquake locations derived from earthquake waves on seismograms. But the precision of earthquake locations from routine seismic network monitoring is insufficient for imaging the detailed geometry of faults, while modern, high-precision location procedures mainly improve locations on only the finest, sub-kilometer scales. Here we apply a new, multi-scale, high-precision earthquake location procedure to earthquake activity along major strike-slip faults in California. Our relocations reveal main fault zones at seismogenic depth as smooth, planar or arcuate, near-vertical surfaces across the sub-km to 10’s of km scales. These results suggest that multi-scale smooth faulting is characteristic of major strike-slip faultzones,likelyinfluencesearthquakeinitiation,rupture,rupturedirectionandarrest,andmayevenbenec-essaryforlargeearthquaketooccur.Theresultscanaidinmappingearthquakehazard,andunderlinethat surface fault traces mainly reflect complex, secondary, shallow deformation and not directly simpler, large earthquake slip surfaces at depth.

The roughness of natural faults has been idealized as fractal or self-affine (Aviles et al., 1987;Power et al., 1988;Scholz, 2019). Fault surfaces on smaller scales (e.g., up to~100 m) are often close to statistically selfsimilar but acquire scale dependence as a result of wear, at least along the direction of slip (Power et al., 1988;Renard et al., 2006;Sagy et al., 2007). However, the magnitude and sense of any roughness scaling depends on the definition of roughness, how it is measured, and if re-ferring to single or multiple fault strands (Beeler, 2023). Surface mapping and exhumed faults suggest faults are rough or corrugated at all scales, from sub-millimeter to hundreds of km, but with a scale dependence (Candela et al., 2012;Renard and Candela, 2017;Beeler, 2023). On scales of up to~100 m the roughness of exposed fault surfaces is found to decrease with total slip (Sagy et al., 2007). Larger scale surface mapping implies a reduction in fault system complexity with increasing geologic offset (Wesnousky, 1988;Stirling et al., 1996;Perrin et al., 2016;Manighetti et al., 2021). And fractal analysis of the main San Andreas, California fault trace indicates it is simple or planar at scales larger than about 1-2 km (Aviles et al., 1987).
The seismogenic depth, a brittle zone where most co-seismic slip and energy release occurs, is from 4 km or less to 10-15 km for many large faults in California (Sibson, 1982;Marone and Scholz, 1988). At these depths, standard, arrival-time relocations of background seismicity and aftershock sequences are usually too diffuse to constrain multi-scale and detailed fault geometry. For the~60 km Parkfield segment of the strike-slip, San Andreas fault in California, various high-precision, differential-timing relocations image a twisting (Thurber et al., 2006;Kim et al., 2016;Perrin et al., 2019), straight (Simpson et al., 2006) or predominantly planar (Thurber et al., 2006) surface, and, on the km scale, as multiple active fault patches offset by tens to hundreds of meters perpendicular to the overall fault surface (Nadeau et al., 1994;Waldhauser et al., 2004). However, the geometry of faults segments and the main rupture zones of larger earthquakes are usually modeled as (Savran and Olsen, 2020;Ramos et al., 2022), and often imaged as (Cockerham and Eaton, 1984;Schaff et al., 2002;Graymer et al., 2007;Waldhauser and Schaff, 2008;Lomax, 2020a) smooth, near-planar surfaces.
Thus, current understanding of the multi-scale geometrical complexity and smoothness of faults and earthquake rupture zones at seismogenic depth is based on conflicting results (see also Goebel et al., 2014). This shortcoming hinders better understanding of earthquake rupture physics and of the relations between faulting at seismogenic depth and surface traces, with important ramifications for earthquake hazard assessment.
Seismicity, the distribution of earthquakes in space, time, and size, is fundamental for understanding earthquakes and for earthquake hazard assessment and forecasting (e.g. Scholz, 2019). Seismicity can show the geometry and activity of faults, the stages of earthquake initiation, and the extent of large earthquake rupture. In particular, seismicity can provide detailed information on fault geometry at seismogenic depth, including on and around surfaces of main, co-seismic slip and energy release. However, relative to the needs of modern seismological study, standard, arrival-time based earthquake locations often have low accuracy and precision, where accuracy is closeness to a usually unknown ground-truth, and precision is relative location accuracy-the correctness of the relative positions of nearby hypocenters. Useful and unbiased determination of the geometry, complexity and smoothness of faults from seismicity requires earthquake location with uniformly high precision over multiple scales.
NLL-SSST-coherence (Lomax and Savvaidis, 2022), a recently developed, arrival-time earthquake location procedure, iteratively generates traveltime corrections to improve precision over many scales (e.g., from the size of a study area to~1 km) and then uses waveform similarity to further improve precision on the finest scales (e.g. sub-km). NLL-SSST-coherence provides multi-scale, high-precision earthquake location.
Here, extending the work of (Lomax and Henry, 2022), we apply NLL-SSST-coherence to relocate recent large-earthquake sequences and background seismicity on major strike-slip faults in and around California. The relocated seismicity at depth surrounding high-slip patches of large earthquakes and on long stretches of major fault zones generally defines planar or arcuate, near-vertical surfaces that are multi-scale smooth (i.e., have fractal or Hausdorff dimension~2.0 over a specified range of length scales). These results have implications for understanding of earthquake rupture physics, fault zone maturity, hazard, and maximum size, for the relationship of surface traces and paleo-seismic results to faulting at seismogenic depth, and for earthquake rupture modeling.
We obtain multi-scale high-precision earthquake relocations through the combined use of source-specific, station traveltime corrections (SSST) and stacking of probabilistic locations for nearby events based on inter-event waveform coherence (Lomax and Savvaidis, 2022). We use the NonLinLoc location algorithm (Lomax et al., 2000, which performs efficient, global sampling to obtain an estimate of the posterior probability density function (PDF) in 3D space for hypocenter location. This PDF provides a comprehensive description of likely hypocentral locations and their uncertainty, and enables application of the waveform coherence relocation. Within NLL, we use the equal differential-timing (EDT) likelihood function (Zhou, 1994;Font et al., 2004;Lomax, 2005Lomax, , 2008Lomax et al., 2014), which is highly robust in the presence of outlier data caused by large error in phase identification, measured arrival-times or predicted traveltimes. We use a finite-differences, eikonal-equation algorithm (Podvin and Lecomte, 1991) to calculate gridded P and S traveltimes for initial NLL locations.

Source-specific station term corrections
In a first relocation stage, NLL-SSST-coherence iteratively develops SSST corrections on collapsing length scales (Richards-Dinger and Shearer, 2000;Lomax and Savvaidis, 2022), which can greatly improve, multiscale, relative location accuracy and clustering of events (Pavlis and Hokanson, 1985b;Richards-Dinger and Shearer, 2000;Nooshiri et al., 2017). In contrast to station static corrections (Tucker et al., 1968;Ellsworth, 1975;Frohlich, 1979;Lomax, 2005Lomax, , 2008, which give a unique time correction for each station and phase type, SSST corrections vary smoothly throughout a 3D volume to specify a sourceposition dependent correction for each station and phase type. These corrections account for 3D variations in velocity structure and corresponding distortion in source-receiver ray paths. NLL-SSST uses smooth, Gaussian distance kernels for accumulating SSST cor-rections, while Richards-Dinger and Shearer (2000) use a fixed number of neighboring events, and  use fixed distance and shrinking-box approaches.
Spatial-varying, SSST corrections are most effective for improving relative locations on all scales when the ray paths between stations and events differ greatly across the studied seismicity, including when stations are inside the seismicity distribution, the extent of seismicity is large relative to the distance to the stations, or the depth range of events is large. SSST corrections can improve multi-scale precision when epistemic error in the velocity model is large, such as when a 1D, laterally homogeneous model or a large-wavelength, smooth model is used in an area with sharp, lateral velocity contrasts or smaller scale, 3D heterogeneities.

Waveform coherency relocation method
In a second relocation stage, NLL-SSST-coherence reduces aleatoric location error by consolidating information across event locations based on waveform coherency between the events (Lomax and Savvaidis, 2022). This coherency relocation, NLL-coherence, is based on the concept that if the waveforms at a station for two events are very similar (e.g. have high coherency) up to a given dominant frequency, then the distance separating these events is small relative to the seismic wavelength at that frequency (e.g. Geller and Mueller, 1980;Poupinet et al., 1984), perhaps less than about a quarter of this wavelength (Geller and Mueller, 1980;Thorbjarnardottir and Pechmann, 1987). A pair of similar events is a doublet and a set of similar events may be called a cluster, multiplet or family; these events all likely occur on a small patch of a fault with similar magnitude and source mechanism (Gedney, 1967;Hamaguchi and Hasegawa, 1975;Ishida and Kanamori, 1978;Geller and Mueller, 1980;Poupinet et al., 1982;Nadeau et al., 1994;Cattaneo et al., 1997;Ferretti, 2005). In a high-precision, microseismic study Goertz-Allmann et al. (2017) show for waveform windows spanning both P and S waves that correlation coefficients greater than about 0.7 indicate event multiplets locate within about 0.1 km, which is about a quarter wavelength for the typical dominant waveform frequency~20 Hz and wave velocity of~2.5 km/s shown in their study. The results of Goertz-Allmann et al. (2017) (their figs. 4 and 6) also show lack of clustering and large separation of event pairs for correlation coefficients less than about 0.5.
For detailed seismicity analysis, the precise hypocenter locations of events in multiplets can be assigned to a unique centroid point or coalesced in space through some statistical combination of the initial hypocenter locations (Jones and Stewart, 1997;Kamer et al., 2015). Alternatively, precise, differential times between likephases (e.g., P and S) for doublet events can be measured using time-or frequency-domain, waveform correlation methods. Differential times from a sufficient number of stations for pairs of doublet events allows high-precision, relative location between the events, usually maintaining the initial centroid of the event po-sitions (Nakamura, 1978;Poupinet et al., 1982Poupinet et al., , 1984Ito, 1985;Got et al., 1994;Nadeau et al., 1994;Waldhauser and Ellsworth, 2000;Matoza et al., 2013;Trugman and Shearer, 2017).
NLL-coherence uses waveform similarity directly to improve relative location accuracy without the need for differential time measurements or many stations with waveform data. The method assumes that high coherency between waveforms for two events implies the events are nearly co-located, and that all of the information in the event locations, when corrected for true origin-time shifts, should be nearly identical in the absence of noise. Then, stacking over probabilistic locations for nearby events can be used to reduce the noise in this information and improve the location precision for individual, target events. We measure coherency as the maximum, normalized cross-correlation between waveforms from one or more stations for pairs of events within a specified distance after NLL-SSST relocation (5-10 km in this study). We take the maximum station coherence between the target event and each other event as a proxy for true inter-event distances and thus as stacking weights to combine NLL-SSST location probability density functions (PDF's) over the events. In effect, this stack directly improves the hypocenter location for each target event by combining and completing arrival-time data over nearby events and reducing aleatoric error in this data such as noise, outliers, and missing arrivals. For a ground-truth test of NLL-SSST-coherence using controlled-source, explosion data from Finland, Lomax and Savvaidis (2022) estimated a relative horizontal location error of about 75 m. See Lomax and Savvaidis (2022) for more discussion and details, and Supplementary File S1 for NLL-SSSTcoherence processing parameters used in this study.
A representation through analogy of the improvement of location precision given by NLL-SSSTcoherence and by cross-correlation based, differentialtiming methods is shown in Figure 1. Relative to a set of true locations, standard catalog locations using arrival-time based location methods contain multiscale distortion primarily due to epistemic error in the velocity model, and smaller scale blurring primarily due to aleatoric error in the arrival-time data. NLL-SSST corrections remove epistemic error to improve multi-scale precision, and NLL-coherence removes aleatoric error to improve smaller scale precision. Differential-timing methods remove mainly aleatoric error in the arrival-time data to improve smaller and finest-scale precision, while, in practice (Waldhauser and Schaff, 2008;Trugman and Shearer, 2017), their time-difference formulation explicitly ignores and does not correct for larger scale, epistemic, velocity model error. These methods thus produce high, finest-scale precision, but do not remove larger, multiscale distortion introduced in underlying, standard catalog locations. No methods can directly improve absolute epicenter and depth shifts and distortions on the largest scale (e.g., the full study area), for which accurate velocity models and calibration with ground truth information is needed.

NLL-SSST-coherence relocations along strike-slip faults in California
We present and discuss NLL-SSST-coherence relocations for well recorded, recent, moderate to large earthquake sequences and background seismicity on major strike-slip faults in and around California ( Figure 2). We do not analyze sequences for very large, California strike-slip earthquakes such as 1993 Mw 7.3 Landers (Hauksson et al., 1993), 1999 Mw 7.1 Hector Mine (Hauksson, 2002) or 2019 Mw 7.1 Ridgecrest (Ross et al., 2019) primarily because aftershock seismicity for such large earthquakes occurs mainly away from main rupture surfaces (e.g. Das and Henry, 2003;Liu et al., 2003), but also, for the earlier Landers and Hector Mine events, due to sparsity of station distribution and lack of available digital waveform recordings.
In most cases we compare the results with highprecision, cross-correlation based, differential-timing relocations for the same time period and magnitude ranges, as available in the Northern California Seismic System HypoDD catalog (NCSS-DD Waldhauser and Ellsworth, 2000;Waldhauser and Schaff, 2008;Waldhauser, 2009). The number of available NCSS-DD events is always less than the number of NLL-SSST-coherence events since HypoDD relocates only events having a minimum number of high cross-correlation connections to nearby events. We analyze the results with a focus on the geometry and smoothness of apparent faulting as imaged or inferred by the multi-scale, NLL-SSST-coherence seismicity. It is important to note that even small heterogeneities in fault geometry, including roughness, kinks, bends, or offsets, can have a large effect on rupture physics (King and Nábělek, 1985;King, 1986;Dieterich and Smith, 2010;Fang and Dunham, 2013).
In the following relocations, we apply NLL-SSST with a smallest, Gaussian kernel smoothing length of 4 km or 2 km and apply NLL-coherence using waveforms up to 10 Hz or 20 Hz frequency (Supplementary File S1). We obtain formal NLL-SSST-coherence epicenter (errH) and depth (errZ) uncertainties as low as 100-200 m (Supplementary Table S1 and Datasets S1-5). This uncertainty range represents the relative locations accuracy (precision) of the NLL-SSST-coherence relocations, but not the absolute location accuracy which may be much larger. As with other location procedures, NLL-SSSTcoherence does not directly improve absolute epicenter and depth accuracy on the largest scale (e.g., the full study area), for which accurate velocity models and calibration with ground truth information is needed. Thus, in the following, "multi-scale" precision ranges from approximately sub-km (as low as 100-200 m) to the extent of each study area.

Smooth, planar faulting: the Parkfield segment of the central San Andreas fault
We first examine the 2004 Mw 6.0 Parkfield sequence and background seismicity along the Parkfield segment of the central San Andreas fault ( Figure 2). The Parkfield segment is at the southeastern end of an over 100

Figure 1
A photo analogy representing the multi-scale improvement in location precision of NLL-SSST-coherence versus the fine-scale improvement of differential-timing methods. NLL-SSST-coherence may not achieve the same precision at finest scales as differential-timing methods such as HypoDD or GrowClust, but can give a better representation of the true geometry of locations across other scales.

Figure 2 Study areas. Study areas in California and
Nevada for NLL-SSST-coherence relocations presented (magenta) and discussed (white) in this work. Green lines show faults from the USGS Quaternary fault and fold database for the United States.
km long, near-straight stretch of the San Andreas fault which exhibits surface creep, and just northwest of a long, locked stretch of the fault that hosted the 1857 M 7.9 Fort Tejon earthquake (Bakun et al., 2005;Langbein et al., 2005). It is generally accepted (e.g. Simpson et al., 2006;Thurber et al., 2006) that there is a single main fault surface at seismogenic depth for the 2004 rupture, and that this rupture falls not along the curved, main San Andreas surface trace, but instead along and under the straighter, Southwest Fracture zone ( Figure 3).
For the Parkfield area, NCSS-DD HypoDD relocations (NCSS- DD Waldhauser and Schaff, 2008;Waldhauser, 2009) show a fault geometry that is kinked and smallscale segmented (Figure 3ac Waldhauser et al., 2004), a common result for high-precision, differential-timing relocation. In contrast, NLL-SSST-coherence relocations ( Figure 3bd) show a much smoother, near-planar fault surface across scales from sub-km to the~50 km study extent. The difference in geometry for the two sets of relocations is emphasized by a singular value decomposition (SVD) fit of a single plane to each of the respective hypocenter sets (Figure 3cd); the mean absolute deviation of NLL-SSST-coherence hypocenters from their SVD plane (140 m) is 54% of that for NCSS-DD (260 m). These characteristics and differences between the two sets of relocations are accentuated in stretched views of the seismicity (Supplementary Figure S1) and in animated, rotating views along the San Andreas fault zone (Supplementary Movies S1 and S2).
The near-planar surface defined by the NLL-SSSTcoherence relocation follows the overall trend of surface faults on the largest scales, but not the smaller scale segmentation and complexity of these faults. At seismogenic depth the NLL-SSST-coherence relocations show no offsets or bends below the epicenters of the 1966 M 6 or 2004 Mw 6.0 Parkfield earthquakes. These two events ruptured nearly the same fault area but initiated at opposite ends of this area and propagated in opposing directions (Bakun et al., 2005); such differences in initiation point and rupture direction may be possible due to the planarity and smoothness of the fault at depth, i.e., fault complexity was not a controlling factor for initiation and other rupture characteristics (Bakun et al., 2005).

Smooth, arcuate faulting: the southern Calaveras fault zone
We next examine the 1984 Mw 6.2 Morgan Hill, California sequence and background seismicity along a 90 km stretch of the southern Calaveras fault zone (Figure 4). The southern Calaveras fault zone, which exhibits shallow creep, branches towards the north from the north end of the creeping section of the San Andreas fault (Watt et al., 2014). Since the installation of dense seismometer networks in the 1970's along this stretch of the Calaveras fault zone there has been abundant micro-  (Reasenberg and Ellsworth, 1982) and 1984 Mw 6.2 (Cockerham and Eaton, 1984) (Oppenheimer et al., 1990), and the 2007 M 5.4 Alum Rock earthquake. For this area, high-precision, NCSS-DD differential-timing relocations (Figure 4a; see also Schaff et al., 2002) again show a kinked and segmented character for the main lineation of seismicity. In contrast, NLL-SSST-coherence relocations (Figure 4b) form a smoother, arcuate lineation on intermediate and larger scales, especially along and around the 1984 Mw 6.2 aftershock zone (yellow events in figure). This arcuate lineation follows closely the circumference of a circle of radius 428 km centered to the south-southwest (Supplementary Figure S2), the significance of which we discuss later.
Neither set of relocations shows a clear relation of seismicity to the complex multitude of surface mapped faults, beyond similar, largest scale trends. And neither set shows a bend in the fault at seismogenic depth near the 1979 M 5.8 or 1984 Mw 6.2 hypocenters; such bends on surface fault traces have been proposed as related to the rupture initiation point for these and other earthquakes (Bakun, 1980;King and Nábělek, 1985). Both events ruptured to the southeast (Figure 4), with the 1984 Mw 6.2 rupture terminating to the southeast where both sets of relocations show clustered, shallow, offfault aftershock seismicity and a possible small offset or kink at depth, while the 1979 M 5.8 main rupture likely terminated at a right step in fault segments, with later aftershocks along the segment further to the southeast (Reasenberg and Ellsworth, 1982;Oppenheimer et al., 1990). The 2007 M 5.4 event also ruptured to the south-east (Oppenheimer et al., 2010).

Smooth faulting in oceanic crust: Mendocino triple-junction, California
Next, we consider seismicity from 1995-2022 around the Mendocino triple-junction, Norther California (Figure 2). This area hosts a fault-fault-trench triple junction with complex, 3D plate interactions: dextral Pacific -Gorda plate motion across the Mendocino fault zone (MFZ), oblique subduction of the east-dipping Gorda plate under the North American plate along the Cascadia subduction zone, and dextral North American -Pacific plate motion across the San Andreas fault (Smith et al., 1993). For the Mendocino triple-junction area, highprecision, differential-timing NCSS-DD relocations (Figure 5a) show diffuse lineations of seismicity offshore along the MFZ and in the underlying, subducting Gorda plate. In contrast, NLL-SSST-coherence relocations ( Figure 5b) show a smooth, narrow, gently curved distribution of hypocenters along the MFZ around 20 km depth (yellow events in figure) and, within the Gorda plate, show several narrow,~20-30km deep, NW-SE lineations of events suggesting smooth or linear sets of parallel fractures (Gong and McGuire, 2021;Lomax and Henry, 2022). For these relocations, NLL-SSST-coherence precision may only be 500 m to 1 km in some areas due to poor station coverage (Supplementary Table S1).

Smaller scale, smooth, planar faulting: the 1986 M 5.7 Mount Lewis sequence
We next examine on a smaller scale the 1986 M 5.7 Mount Lewis sequence and surrounding, 1984-1999 background seismicity (Figure 2). This sequence occurred just to the north of our southern Calaveras fault zone study in an~25 km, north-south area with no mapped surface faults. The 1986 M 5.7 mainshock had a north-south oriented, right-lateral strike-slip mechanism and a highly productive aftershock sequence within in a distinctive, north-south oriented "hourglass" shaped volume (Zhou et al., 1993;Dodge et al., 1996;Kilb and Rubin, 2002). For the Mount Lewis sequence, high-precision, differential-timing NCSS-DD relocations (Figure 6a), more clustered and organized, NLL-SSST-coherence relocations ( Figure 6b) and the results of Kilb and Rubin (2002) all define well the extensive, volumetric, hourglass form of seismicity to the north and south of the mainshock hypocenter. All sets of relocations also show a narrow, central,~2 km long (NCSS-DD) to~3 km long (NLL-SSST-coherence; see Supplementary Movie S5) north-south, tabular trend of foreshocks (blue) and early aftershocks (yellow) around the mainshock hypocenter, extending from about 2 km above to 1 km below the hypocenter (large yellow dot). Kilb and Rubin (2002) interpret this trend as a kinked mainshock rupture surface. Here, SVD analyses of events within a 1.2 km wide rectangular prism centered on the tabular trends shows that near vertical planes fit well (mean absolute deviation 66 m for NCSS-DD, 35 m for NLL-SSST-coherence) all foreshocks, mainshock and early aftershocks while covering an area similar to that inferred through teleseismic waveform analysis for the mainshock rupture (Zhou et al., 1993). For the NLL-SSST-coherence relocations, the aftershocks just north and south of the SVD plane are mainly offset east and west, respectively, from the strike of the plane, as expected for aftershocks concentrating in the extensional quadrant of a right-lateral, strike-slip event (Kim et al., 2004). These results suggest a simple, smooth, planar surface for the main M 5.7 rupture, while most aftershocks occur outside this surface on extended, complex secondary structures, including fault sets perpendicular to main rupture, indicating an immature fault system (Kilb and Rubin, 2002).

Indirect indication of smooth, planar faulting: southwest of San Francisco
We next consider background seismicity along~50 km of the San Andreas fault zone (SAFZ) to the south and west of San Francisco (Figure 2). The SAFZ to the west of San Francisco likely hosts the hypocenter of the M 7.9 1906 California earthquake (Lomax, 2008). NLL-SSSTcoherence relocations from Feb 1981 to April 2021 for this area are shown in Figure 7. A major part of the seismicity forms an~35 km long zone at about 5-11 km depth which is rotated about 5 • clockwise to and crosses under the surface expression of the SAFZ. This seismicity has predominantly extensional focal mechanisms and is associated with an extensional, right stepover between the onshore San Andreas fault and the offshore Golden Gate fault (Zoback et al., 1999;Parsons, 2002;Lomax, 2008).
To the south of San Francisco, the onshore surface trace of the San Andreas fault (SAFZ and cyan line in Figure 2a) is nearly linear and exhibited up to 4.5 m of rupture during the M7.9 1906 earthquake (Reid and Lawson, 1908). These relations and evidence from seismicity (Zoback et al., 1999) and reflection seismics (Hole et al., 1996) suggest that the active fault surface that hosted 1906 rupture at seismogenic depth may be represented by a vertical plane under and along the surface trace. Such a plane delimits well the northeastern boundary of the extensional seismicity along this segment (Figure 7ab), as also found by Zoback et al. (1999) and Lomax (2008). The truncation of ongoing extensional seismicity along a vertical fault may indicate a strong contrast in geologic structure across the fault (e.g. Liu et al., 2003) such that present-day background stress leads to, in this case, distributed, normal faulting to the southwest while the northwest side remain mainly aseismic.
Along and below the northeast boundary of the extensional seismicity there is a 30 km long set of seismicity clusters at around 11-13 km depth (Figure 7ac) with mainly strike-slip focal mechanisms (Lomax, 2008). This set of deep clusters is well fit by a linear trend rotated about 5 • clockwise to the SAFZ (white rectangle in Figure 7c) and which appears to connect the San Andreas and Golden Gate faults below and across the ex- tensional stepover. This linear trend of deep seismicity suggests either a linear, strike-slip fault structure around 12 km depth, or localized, synthetic, or antithetic faulting at the base of a brittle crust in response to dextral shear across an underlying, linear ductile zone (Lomax, 2008). An upwards extension of the linear trend of deep seismicity appears to delimit the northeastern boundary of shallower extensional seismicity (Figure 7a), but examination of the geometrical relations in 3D and structural considerations favor that this boundary is right-stepping with segments parallel to the main SAFZ.

NLL-SSST-coherence methodology
NLL-SSST and NLL-coherence together greatly increase precision (relative location accuracy) on multiple scales within a standard, arrival-time location framework (Fig-Figure 6 1984-1999 Mount Lewis sequence relocations. M ≥ 1.0 hypocenters from 1984-01-01 to 1999-12-31 for (a) 3343 NCSS-DD and (b) 3503 NLL-SSST-coherence event relocations in map view (upper panels) and lateral view from east (lower panels). Hypocenter color shows origin time (cyan events show foreshocks before the M 5.7 mainshock (large yellow dot), yellow events the first 7 days of aftershocks), symbol size is proportional to magnitude. Inverted pyramids show nearby seismic stations used for relocation. White rectangles show plane of best SVD fit to hypocenters in a 0.6 km wide zone around the rectangle, heavy line indicates top of the plane. Background topography image from OpenTopography.org. See also Supplementary Movie S5. ure 1; Lomax and Savvaidis, 2022). Building on and making use of the thorough, probabilistic global sampling and robust, EDT likelihood function of NLL, NLL-SSST improves multi-scale precision by iteratively removing common-mode traveltime residuals at available stations as a function of hypocentral position. This procedure reduces epistemic model errors and location bias between nearby events located with differing sets of stations or phase types. NLL-coherence location improves smaller scale precision by stacking probabilis- tic, NLL-SSST location PDF's of nearly co-located, multiplet events, as measured by waveform similarity. This stacking of PDF's effectively reduces aleatoric data error and suppresses outliers in the underlying arrival times, while filling in missing arrival time data across multiplet events, resulting in a data-driven, spatial coalescence of location for events with similar waveforms. In the following, as explained earlier, "multi-scale" precision ranges from approximately sub-km (as low as 100-200m) to the extent of each study area.
In contrast to the coherence-weighted stacking of PDFs for nearby events in NLL-coherence, crosscorrelation based, differential-timing methods such as HypoDD or GrowClust achieve high to very high, finescale precision through explicit, inter-event, differential location. This location involves over nearby event pairs of differences in distance along event-station ray directions, as constrained by all available arrival-time differences and the used velocity model. For relocation studies with good station coverage, and thus good ray coverage around the events, these differential-timing methods should achieve higher, finest-scale precision than NLL-SSST-coherence. However, for cases of poor station and ray coverage, NLL-SSST-coherence may retain higher relative location accuracy and better depth control than do cross-correlation based, differentialtiming methods, as indicated by our results for Mendocino triple-junction and Mount Lewis seismicity and previous results for the 2020 Mw 5.8 Lone Pine, California sequence (Lomax and Savvaidis, 2022). Furthermore, differential-timing methods such as HypoDD or GrowClust usually preserve or allow only slight change in the centroid of the starting locations for individual clusters of high similarity events (Waldhauser and Ellsworth, 2000;Waldhauser and Schaff, 2008;Trugman and Shearer, 2017) and thus, in general, will not improve larger, multi-scale precision beyond that of the starting locations (Figure 1). It is possible that applying cross-correlation based, differential-timing methods after NLL-SSST relocation would produce optimal multiand finest-scale location precision. And applying these methods after NLL-SSST-coherence may also have advantages, as NLL-coherence can collect noisy, outlier locations back into their correct clusters. The XCORLOC method (Neves et al., 2022) precedes cross-correlation based, differential-timing relocation with L1 and L2 norm SSST and so can improve multi-scale precision in a manner analogous to NLL-SSST-coherence.
There is no evident reason why the NLL-SSSTcoherence methodology might smooth hypocenters locations over larger distances as an artifact. Indeed, the iterative, multi-scale NLL-SSST procedure generates smaller-scale traveltime corrections that are independent over larger distances; this independence should preserve true roughness or offsets in alignments of seismicity. For the case of large error in the arrivaltime data, the independence of NLL-SSST corrections over larger distances would more likely produce artifact error and offset between clusters of hypocenters than artifact smoothing. For Mount Lewis, the NLL-SSST-coherence relocations match well detailed features of the high-precision NCSS-DD ( Figure 6) and Kilb and Rubin (2002) relocations, including definition of the main rupture surface, complex secondary structures and fault sets perpendicular to the main rupture, without exhibiting additional smoothing or smearing of these features that might be artifacts of the NLL-SSSTcoherence methodology.
To further illustrate this point, we compare NLL-SSST-coherence with high-precision, differentialtiming relocations based on a precursor to GrowClust (HYS Hauksson et al., 2012) for the complex, 2021 Mw 5.3 Calipatria sequence (Supplementary Figure S3). This sequence, in the Brawley Seismic zone at the southern end of the San Andreas fault system, ruptured 3 larger, near-orthogonal segments overs scales of~2-10 km and numerous smaller scale features over about 7 days (Hauksson et al., 2022). NLL-SSST-coherence and HYS relocations of the Calipatria seismicity (Supplementary Figure S3) closely reproduce the same features over all scales, including the larger scale, near-orthogonal planes, and smaller scale splay and clustered seismicity, while both set of relocations show similar depth distribution of these features, including shallowing of a sharp base of seismicity towards the north. These results show that, besides larger scale distortion due to different velocity models and station correction procedures, NLL-SSST-coherence does not oversimplify or smooth distributed, multi-scale, multi-fault ruptures compared to high-precision, differential-timing relocations.
Similar to NLL-SSST, 3D, tomographic, velocity model inversions implicitly generate station and source dependent traveltime corrections, and may involve collapsing scale lengths. However, differential-timing relocation in 3D tomographic models for Parkfield relocations (Thurber et al., 2006) do not show as smooth, vertical, or planar a surface as does NLL-SSST-coherence (Figure 3b). This difference may be related to the mapping in tomographic inversion of all traveltime residuals to unique, 3D, Vp and Vs velocity grids, while SSST maps the residuals to a large set of 3D, station Vp and Vs traveltime grids; this latter procedure retains many more degrees of freedom and so can preserve a greater amount of information from the residuals that may be useful for precise location.

Relocation results
With relocation of earthquake sequences and background seismicity along the San Andreas (Parkfield) and southern Calaveras strike-slip faults in California, we have shown that NLL-SSST-coherence relocated seismicity at seismogenic depth along major faults and surrounding large-earthquake ruptures often defines narrow, multi-scale smooth, planar (e.g., Parkfield) or arcuate (e.g. southern Calaveras), near-vertical surfaces across multiple scales. For Parkfield, the high-precision relocations of Thurber et al. (2006) (see interpretation of Simpson et al., 2006), and the XCORLOC relocations of Neves et al. (2022) also suggests, on the intermediate and largest scales, smooth, near-vertical faulting.
NLL-SSST-coherence relocations for the Mendocino triple-junction area show that such multi-scale smooth faulting also occurs for strike-slip faulting in oceanic crust, as found in other areas (e.g. Schlaphorst et al., 2023). NLL-SSST-coherence relocations to the southwest of San Francisco suggests a deep, linear fault or response to a deeper, linear shear zone over~30 km. For the 2019 Mw 7.1 Ridgecrest, California, sequence highprecision (e.g. Ross et al., 2019;Shelly, 2020) and standard NLL locations (Lomax, 2020a) define two planar, orthogonal faulting surfaces for main rupture of the Mw 6.4 foreshock, while Lomax (2020b) additionally determines that these planes are at different depths and nonintersecting.
On a smaller scale, for NLL-SSST-coherence relocations of the 1986 M 5.7 Mount Lewis sequence (Figure 6b), a near vertical,~4x4 km plane fits well the foreshocks, mainshock and early aftershock hypocenters and covers an area similar to that of mainshock rupture as inferred through teleseismic waveform analysis. An hourglass form of aftershocks to the north and south of this simple, planar mainshock rupture surface, sug-gest complex splay and cross faulting which falls mainly withing the dilatational quadrants of right-lateral mainshock rupture as delimited by the~4x4 km plane.
The 2020 Mw 6.5 Monte Cristo Range, Nevada sequence occurred along an immature, strike-slip fault zone, with no clear surface traces related to primary rupture (Koehler et al., 2021). At seismogenic depth, NLL-SSST-coherence relocations for this sequence (Lomax, 2020a) define two, en-echelon, smooth, planar faulting segments corresponding in size and orientation to expected and modeled mainshock rupture, as well as lateral and shallow secondary and splay faulting forming an extensive damage zone, in similarity to the Mount Lewis results.
For NLL-SSST-coherence relocations to the south of San Francisco, the northeast limit of diffuse, extensional seismicity corresponds to a multi-scale smooth, vertical, planar fault along the SAFZ. In this case, the seismicity does not directly fall on and define an active surface of faulting, but instead delimits the edge and depth limits of the likely 1906 faulting surface that is currently mainly aseismic. Preliminary analysis of NLL-SSST-coherence relocations for the 2014 M 6.0 South of Napa, California sequence and background seismicity shows similar, though less clear, indirect indications of planar faulting for the M 6.0 mainshock rupture.
These NLL-SSST-coherence relocations define multiscale smooth faulting over at least 10's of km for segments of mature, strike-slip fault zones, and over the likely rupture zones of large and moderate earthquakes along these faults and on less mature faults. These results suggest that multi-scale smooth (down to sub-km lengths), planar and arcuate faulting is characteristic of strike-slip fault zones at seismogenic depth and perhaps necessary for larger earthquake rupture. We next consider some important implications of these results.

Rupture physics
The smoothness and curvature of fault segments at seismogenic depth likely influences earthquake rupture physics-initiation, rupture, and arrest (Okubo and Dieterich, 1984;Ben-Zion and Sammis, 2003;Dieterich and Smith, 2010;Fang and Dunham, 2013), and perhaps enables the occurrence of larger earthquakes (Goebel et al., 2017(Goebel et al., , 2023. Earthquake initiation may be possible anywhere within smooth fault segments, though most likely at non-geometrical asperities-areas of stress concentration either within the segments due to previous rupture history or material heterogeneities and perhaps indicated by concentrations of microseismicity, or at limits of the segments including at kinks or stepovers (Das and Henry, 2003;Aki, 1979;King, 1986;Scholz, 2019). The arrest of earthquake rupture is likely favored at barriers such as kinks, steps, or other fault complexities at the limits of smooth segments at seismogenic depth, as well as within smooth segments at areas of stress relaxation due to previous rupture history (Sibson, 1985;King, 1986;Scholz, 2019). Laboratory experiments indicate that smooth faults have larger coseismic slip, lower residual stress and fewer aftershocks compared to rough faults (Goebel et al., 2023).
Most importantly, for faults that are smooth and planar or horizontally arcuate (especially when following closely the circumference of a circle, as is the case for the southern Calaveras; Supplementary Figure S2), there may be negligible geometrical interactions and resulting backstresses (Dieterich and Smith, 2010) impeding strike-slip rupture displacement. In this case, minimal energy is absorbed by off-fault inelastic deformation as fracture energy (Cocco et al., 2023), and a maximum of strain energy released during rupture would be available to further drive the rupture. Thus, sustained earthquake rupture, and perhaps even the occurrence of larger earthquakes, would be more likely over a smooth fault surface than a rough or fine-scale segmented surface (e.g. Dieterich and Smith, 2010;Fang and Dunham, 2013;Perrin et al., 2016).
Multi-scale smooth faults are also considered most likely to support sustained supershear rupture propagation (Bouchon et al., 2010;Bruhat et al., 2016), and should radiate relatively less high-frequency energy than rough faults (Madariaga, 1977;Shi and Day, 2013;Trugman and Dunham, 2014). For planar faulting, the relative displacement between crustal blocks would be planar shear, while the displacement of blocks across arcuate faults would include a rotational component.
In addition, for a multi-scale smooth, curved fault, there may be a preferred direction for rupture (Rubin and Gillard, 2000), perhaps due to the position of the curved fault surface forward of rupture relative to the dilatational and extensional quadrants of the rupture (Fliss et al., 2005). This is suggested by our results for the southern Calaveras fault zone where the 1984 Mw 6. 2, 1979 M 5.8 and 2007 M 5.4 events all ruptured to the southeast (Figure 4b). Given the sense of curvature of NLL-SSST-coherence seismicity and right-lateral slip, for southeastward rupture the fault forward of the rupture front is in the dilatational quadrant of strain from the current and previous rupture. This dilatational strain would decrease normal stress across the fault, producing dynamic "unclamping" in front of the rupture and facilitating further slip, in a manner analogous to a continuum of infinitesimal, extensional bends or stepovers (Poliakov et al., 2002;Oglesby, 2005;Oglesby and Mai, 2012;Parsons and Minasian, 2015). Under this mechanism, rupture in the opposite direction, northwest in this case, would be impeded, as with infinitesimal, compressional bends or stepovers, which may explain why 1979 M 5.8 rupture did not propagate or trigger slip to the northwest into the presumably well loaded, future rupture zone of the 1984 Mw 6.2 event. A general rule would be that rupture is promoted in the direction along which the fault is concave to the right for right-lateral slip and concave to the left for left-lateral slip. This mechanism could explain sense of rupture for other large earthquakes on smoothly curving segments of strike-slip faults, including for southeastwards rupture for the 2002 Mw 7.9 Denali fault, Alaska-Canada earthquake (Eberhart-Phillips et al., 2003), for primarily eastwards rupture of the 1943 Ms 7.7 Tosya, Turkey earthquake along the North Anatolian fault (Dewey, 1976;Barka and Kadinsky-Cade, 1988;Stein et al., 1997) and perhaps for at least the first half of westward rupture for the 1939 Ms 7.9 Erzincan, Turkey earthquake (Emre et al., 2021) where the North Anatolian fault appears mildly concave to the north (Emre et al., 2018). A preferred, eastward rupture direction due to fault curvature for the 1943 Tosya earthquake could in part explain why it ruptured in the opposite direction to that of other major events along the North Anatolian fault in the past century and why its epicenter is, anomalously, not in an area of increased stress from previous events (Stein et al., 1997). For a planar fault, this mechanism is not active and gives no preferred rupture direction, which is consistent with the 1966 M 5.5 and 2004 Mw 6.0 Parkfield earthquakes rupturing in opposing directions from opposite ends of a segment of the San Andreas fault that our NLL-SSST-coherence results show is smooth and planar (Figure 3b).
Most of these relations between rupture physics and smooth faulting may apply to individual fault segments of any size in a self-similar manner, so that, for example, the smaller is a segment of smooth faulting, the smaller is the largest rupture that can occur on that segment. Thus, while seismicity on secondary, splay and damage zone faulting, such as around the apparently smooth main rupture segments for the Mount Lewis and Monte Cristo (Lomax, 2020b) sequences, and aftershock seismicity in general may not exhibit larger scale patterns suggesting fault smoothness, individual events and localized clusters could follow these relations on the (small) scale of their rupture segments.
On larger scales, many complex and very large earthquakes involve rupture on a number of separate fault segments, each of which may be multi-scale smooth, as is the case for Monte Cristo, and possibly the case for the 2016 Mw 7.8 Kaikōura, New Zealand earthquake (Hamling et al., 2017), and for the 1993 Mw 7.3 Landers (Hauksson et al., 1993) and 1999 Mw 7.1 Hector Mine (Hauksson et al., 2022) earthquakes in California. However, aftershock seismicity for large earthquakes occurs mainly or almost entirely off of main rupture surfaces (e.g. Das and Henry, 2003;Liu et al., 2003;Goebel et al., 2023) making it difficult to define precisely the geometry of main rupture surfaces for very large earthquakes.
It is possible that on the largest scales, e.g., for the 1857 M 7.9 and 1906 M 7.9 rupture zones along the San Andreas fault, main rupture may occur on one or few long, smooth segments. In this case features of rupture physics of smooth faults discussed above combined with a possible strong locking of smooth faults due, for example, to efficient healing by cementation or other processes on thin, smooth fault surface (Muhuri et al., 2003;Williams and Fagereng, 2022), may be explanations for potentially long recurrence intervals and resulting large size of these events. Additionally, the seismicity patterns southwest of San Francisco, where upper crustal seismicity occurs mostly off the SAF with dominantly extensional focal mechanisms, is consistent with near complete release of shear stress after the 1906 earthquake. Thus, for a given length scale, a smooth segment may be expected to take more time than a rough segment to reload back to a critical state after rupture. This self-similarity likely extends to the concept of fault maturity, characterized by increased simplification including smoothing and reduced extent of lateral damage zones (Naylor et al., 1986;Wesnousky, 1988;Scholz, 2019). We find that sequences on immature fault systems such as Mount Lewis and Monte Cristo, though terminated by extensive damage and splay faulting, contain core segments of main rupture which may be planar, multi-scale smooth surfaces with little, lateral damage zone seismicity. Such main rupture surfaces may be mature, on their smaller length and age scales, in the same sense as are much longer segments of major fault systems such as the San Andreas on much larger scales. Such a scale invariance is found by Evans et al. (2000) for faults exhumed from seismogenic depth, which show similarity of shear-induced microstructures and deformation mechanisms for faults 10 m to 10 km long starting from an early age as inferred from total slip.

Earthquake hazard and maximum size
Our overall results suggest identification of stretches of smooth faulting would help identify zones of earthquake hazard and possibly help quantify maximum earthquake size. This identification might be made directly from background seismicity or aftershocks falling on smooth surfaces, or indirectly by the distribution and geometry of clustered, diffuse, or other seismicity which may delimit stretches of aseismic, smooth faulting, as suggested in our analysis of seismicity south of San Francisco. Difficulties arise due to the possibility that future large earthquakes may occur in areas of current seismic quiescence, including gaps or locked patches, for example as indicated by the sparsity of recent seismicity along the 1857 M 7.9 and 1906 M 7.9 rupture zones on the San Andreas fault (Jiang and Lapusta, 2016;Scholz, 2019). In this case, if delimited by stretches of smooth faulting at seismogenic depth, nearsilent segments of known or inferred fault zones might be identified as having elevated hazard. Also, surface mapped fault traces that are rough, multi-stranded or offset, but smooth on average over large length scales may be generated by smooth faulting at seismogenic depth, as is found with analogue fault modeling and interpreted for shallow natural faulting (Naylor et al., 1986;Klinger, 2010;Dooley and Schreurs, 2012), such surface features can therefore indicate elevated hazard.

Faulting at shallow versus seismogenic depth
Our NLL-SSST-coherence relocations along major strike-slip faults mainly concentrate on a single, smooth surface at more than a few km depth below zones with a multitude of surface traces, sometimes offset from these traces. These relations provide further evidence that surface traces and offsets of strike-slip fault zones reflect complex, shallow deformation, perhaps involving braided and upwards diverging fault structures (e.g. Christie-Blick and Biddle, 1985;Richard et al., 1995;Graymer et al., 2007), and not directly simpler and hidden slip surfaces at seismogenic depth (e.g. Michael, 1988;Oppenheimer et al., 1990;Schaff et al., 2002;Ponce et al., 2004;Graymer et al., 2007;Watt et al., 2014;Chaussard et al., 2015), where most co-seismic slip and energy release occurs. Furthermore, NLL-SSST-coherence relocations for the Mount Lewis and Monte Cristo (Lomax, 2020b) sequences, along smaller and immature strike-faults, also define at seismogenic depth one or more smooth faulting surfaces which correspond to probable mainshock rupture surfaces. However, for these cases there are either no mapped surface faults (Mount Lewis) or complex, mapped surface fractures showing little relation to the deeper mainshock rupture (Monte Cristo), again emphasizing an indirect relation of surface features to main rupture surfaces at seismogenic depth.
In general, shallow deformation associated with earthquake ruptures involves significant diffuse anelastic deformation (e.g. Antoine et al., 2023). Several processes may contribute to explain a broadening of the zone of deformation around faults from the seismogenic zone to the surface, forming, for instance, flower structures above strike-slip faults (e.g. Harding, 1985). These include transition from unstable to stable sliding, which limits co-seismic slip on fault surfaces toward the surface (Scholz, 1998), and increasingly distributed damage due to relatively weaker shallow materials and the geometrical effect of the free surface, as shown in numerical and analog models (e.g. McClay and Bonora, 2001;Finzi et al., 2009;Wu et al., 2009;Ma and Andrews, 2010). For basic understanding of large earthquake faulting and hazard it is important to better define the geometrical transitions and physical connections between complex shallow faulting and potentially simpler and smoother fault segments at seismogenic depth.

Earthquake rupture modeling
The occurrence of earthquake rupture on multi-scale smooth faults justifies and would require the use of planar or smoothly curved fault segments for kinematic or dynamic numerical modeling (Ramos et al., 2022) of primary rupture and energy release of an earthquake at seismogenic depth. However, modeling of secondary, splay and damage zone faulting, such as is apparent around the main ruptures for the Mount Lewis and Monte Cristo sequences, likely requires use of complicated and rough model fault geometries (Ramos et al., 2022) or may be better represented by continuum mechanics-based numerical modeling (Preuss et al., 2020) across 3D volumes. Additionally, an indirect relation of surface fault traces to main rupture surfaces at seismogenic depth may preclude simple, downward projection of these shallow traces for rupture modeling; available information from aftershock or background seismicity, and geophysical and geologic studies should always be considered for constructing fault segments at seismogenic depth.

Conclusions
Our NLL-SSST-coherence relocations for California along major, strike-slip faults and surrounding largeearthquake ruptures show narrow, planar, or arcuate, near-vertical, multi-scale smooth faulting at seismogenic depth across the sub-km to 10's of km scales. These results suggest that multi-scale smooth faulting may be a characteristic of segments of major, strikeslip fault zones, of large earthquake rupture within individual fault segments, and, in a self-similar manner, of earthquake ruptures of smaller sizes.
The smoothness and curvature of faults likely influences large earthquake initiation (possible anywhere within or at the limits of smooth fault segments), rupture (multi-scale smooth faults facilitate, and may be required, for large earthquake rupture; if the fault is curved, there may be a preferred direction for rupture), and arrest (favored at kinks, steps, or other nonsmooth fault complexities). Consequently, zones of earthquake hazard can be identified directly from planar and smooth alignments of seismicity and indirectly from patterns in clustered or diffuse seismicity.
Our findings provide further evidence that surface traces and offsets of strike-slip fault zones reflect complex, shallow deformation, and not may not correspond directly to simpler, smoother slip surfaces at depth where most co-seismic slip and energy release occurs. This relation has important implications for earthquake hazard assessment, and supports use of planar or smoothly curved surfaces, but not necessarily the complexity of surface rupture traces, for earthquake rupture modeling. ment and discussion, and Egill Hauksson for assistance with the high-precision, HYS catalog for the 2021 Calipatria sequence. Special thanks to the analysts, technicians and scientists who formed the highquality, NCEDC and SCEDC arrival-time and event catalogs which made the precision of this work possible. We use LibreOffice (https://www.libreoffice.org, last accessed April 2023) for word processing, spreadsheet calculations and drawings, and Zotero (https:// www.zotero.org, last accessed April 2023) for citation management. This work was performed outside any specific funding or grant.

Data and code availability
The supplementary information for this article includes Table S1 and Figure S1-3, and, in a repository (Lomax and Henry, 2023), Movies S1-7 showing animated views of NLL-SSST-coherence relocations presented in the main paper, Datasets S1-6 containing CSV format catalogs of NLL-SSST-coherence relocation results, and File S1 containing run scripts and related set-up, configuration and other meta-data files for locations cases presented in this paper.