The Shear Deformation Zone and the Smoothing of Faults With Displacement

We use high‐resolution earthquake locations to characterize the three‐dimensional structure of active faults in California and how it evolves with fault structural maturity. We investigate the distribution of aftershocks of several recent large earthquakes that occurred on continental strike slip faults of various structural maturity (i.e. various cumulative fault displacement, length, initiation age and slip rate). Aftershocks define a tabular zone of shear deformation surrounding the mainshock rupture plane. Comparing this to geological observations, we conclude that this results from the re‐activation of secondary faults. We observe a rapid fall off of the number of aftershocks at a distance range of 0.06‐0.22 km from the main fault surface of mature faults, and 0.6‐1.0 km from the fault surface of immature faults. The total width of the active shear deformation zone surrounding the main fault plane reaches 1.0‐2.5 km and 6‐9 km for mature and immature faults, respectively. We find that the width of the shear deformation zone decreases as a power law with cumulative fault displacement. Comparing with a dynamic rough fault model, we infer that the narrowing of the shear deformation zone agrees quantitatively with earlier estimates of the smoothing of faults with displacement, both of which are aspects of fault wear. We find that earthquake stress drop decreases with fault displacement and hence with increased smoothness and/or slip rate. This may result from fault healing or the effect of roughness on friction.

The above definitions allow us to distinguish two types of damage zones: the "dilatant damage zone" dominated by volumetric strains, and the "shear deformation zone" dominated by shear strains. The tensile (Mode I) cracks in the dilatant damage zone are dilatant cracks that align parallel to the maximum compression direction and perpendicular to the minimum principal stress. Hence the orientation of cracks provides evidence for the several different mechanisms responsible for their formation (Wilson et al., 2003). The shear deformation zone is characterized by secondary faults (Modes II and III cracks) and hence are oriented parallel to the maximum Coulomb stress. For example, in the case of a strike-slip fault, this zone is defined by a conjugate set of secondary faults (Little, 1995), in which one set is parallel to the primary fault.
The evolution of these three zones in Figure 1 defines what is called fault maturity. The three zones can be viewed as regions controlled by wear processes, and the fault structural maturity can hence be measured by its degree of wear, which depends primarily on the net fault displacement, assuming a constant slip-vector throughout the fault's history. However, previous studies have shown that, in the absence of data on net fault displacement, several other fault parameters such as the fault initiation age and the geological slip rate can be also used as a proxy of net displacement in evaluating the overall maturity of the fault (e.g., Choy et al., 2006;Choy & Kirby, 2004;Dolan & Haravitch, 2014;Hecker et al., 2010;Ikari et al., 2011;Manighetti et al., 2007;Niemeijer et al., 2010;Perrin, Manighetti, Ampuero, et al., 2016;Stirling et al., 1996;Wesnousky, 1988). As these parameters increase, the fault grows and becomes more "mature." Prior studies have suggested that the structural maturity may have a strong impact on earthquake behavior, such as magnitude, stress drop, distribution of slip, rupture velocity, ground motion amplitude, and number of ruptured segments (e.g., Cao & Aki, 1986;Dolan & Haravitch, 2014;Hecker et al., 2010;Malagnini et al., 2010;Manighetti et al., 2007;Perrin, Manighetti, Ampuero, et al., 2016;Radiguet et al., 2009;Stirling et al., 1996;Wesnousky, 1988).
The widths of the fault core and dilatant damage zones saturate at fault lengths comparable to the seismogenic thickness, due to the switch from crack like to pulse like earthquake propagation (Ampuero & Mao, 2017). For larger faults, the evolution of fault maturity involves only changes in the shear deformation zone. In this paper, we are concerned with the scaling of large faults (i.e., which breach the brittle seismogenic width) and their associated large earthquakes, so we are only concerned with the shear deformation zone. There is evidence that indicates that large faults become smoother with net displacement (Stirling et al., 1996;Wesnousky, 1988). This smoothing is probably the prime attribute of fault maturity. Here we show that the width of the shear deformation zone of large faults decreases with fault displacement, as a consequence of this smoothing.
Precise earthquake locations can be used to image the internal structure of faults and the zone of brittle deformation, often at a resolution similar to field observations (e.g., Hauksson, 2010;Powers & Jordan, 2010;Valoroso et al., 2014). Powers and Jordan (2010) studied the association of small earthquakes with large faults in California. They found that the frequency of small earthquakes is highest in a narrow region surrounding faults and then falls off as a power law at greater distances. They modeled this behavior with a fault model with rough (fractal) topography (Dieterich & Smith, 2009), showing that such a rough fault model would produce high stresses near the fault that could account for the seismicity. The earthquakes they used were from the interseismic period of the faults. Although stacked profiles from many fault sections show the association of seismicity with the faults, individual sections and map views, both in Powers and Jordan (2010) and Hauksson (2010) show that such a tight correlation with the fault is not typical for individual fault segments, which often display a wide variability in distribution. As Hauksson (2010) observed, this variability is likely due to many factors, such as heterogeneity in lithology, the effect of nearby faults both mapped and unmapped, and fault offsets and bends.
As Hauksson (2010) observed further, aftershocks of large earthquakes, in contrast, mostly show a tight clustering around the main fault. At depth, they are mainly distributed at the edges of regions experiencing high coseismic slip and in surrounding areas of low coseismic slip (e.g., Mendoza & Hartzell, 1988). Also, earthquakes that occur during interseismic periods often delineate the edges of stuck patches sometimes surrounded by creeping areas, and potential rupture surfaces at depth (e.g., Perrin et al., 2019;Schaff et al., 2002;Waldhauser et al., 2004). In both cases, the majority of aftershocks of large earthquakes sample a portion of the main fault zone structure at depth at the edge of the rupture surface, and are often used to delineate it. In an earthquake model with a smooth rupture surface, stress drop increases to a maximum near the fault surface, thus precluding aftershocks in the near-fault region (Kostrov & Das, 1984). Near-fault aftershocks therefore are an indication of a rough fault, as near-fault stresses are generated by dynamic slip on rough topography. They are greatest right after the mainshock, after which they are relaxed by aftershocks and other relaxation mechanisms. In this study, in order to estimate the roughness of active faults, and how they evolve with displacement, we use high-precision aftershock locations of large earthquakes on faults with different net displacements.

Width of the Shear Deformation Zone From Aftershock Locations
We use high-resolution earthquake catalogs available in the literature to analyze the aftershock distribution of seven large (M w ≥ 6) continental strike-slip earthquakes that occurred on faults with a wide range of net displacement: the 1984 Morgan Hill, 2004 Parkfield, and 2014 South Napa earthquakes in northern and central California (Waldhauser, 2009;Waldhauser & Schaff, 2008Superstition Hills, 1992Landers, 1999Hector Mine, and 2010 El Mayor Cucapah earthquakes in Southern California/Mexico (Hauksson et al., 2012). These earthquakes were selected because high-precision aftershock catalogs, relocated by the double difference method (Waldhauser & Ellsworth, 2000) and with relative lateral location errors derived from bootstrap analysis that are typically less than ≤0.03 km on average ( Figure S1). In order to perform a homogeneous analysis of all sequences, we selected all aftershocks with M L ≥ 1 that occurred within 2 months of each mainshock.
For each aftershock sequence, we determine the three-dimensional fault geometry by applying a principal component analysis (PCA) to all events within boxes that are between 3 and 10 km long along strike and between 3 and 20 km wide across strike (i.e., centered on the surface fault trace), stepping at 1 km intervals along the fault trace (Perrin et al., 2019). For each box, we obtain a plane that best fits the locations of aftershocks. For simplicity, we assume a constant dip of the calculated planes in each box (see also Perrin et al., 2019). In general, the calculated fault plane parameters are in good agreement with existing fault parameters derived from source inversion models (Mai & Thingbaijam, 2014 and references therein; see details in Table S1 and Figure S2).
We then compute the orthogonal distance between each hypocenter and the calculated fault plane segment, count the number of aftershocks as a function of distance from the fault plane within bins of 0.03 km, and normalize the results by the total number of aftershocks in each box (gray lines in Figures 2b, 2d and 2e, 2f, 2g and Figure S3). Finally, we compute the mean of these normalized counts for each distance bin (black line in Figure 2 and Figure S3) to describe the smoothed distribution of aftershocks as a function of distance from the fault plane for each of the seven earthquake sequences.
We use a power law function to fit the mean normalized aftershock distribution (red curve in Figure 3 and Figure S4) to find the two parameters that describe the width of the shear deformation zone. W S1 is the distance from the fault plane where the fall-off in number of events exhibits the maximum curvature (i.e., maximum in second derivative). W S2 is the distance from the fault plane where the number of events departs from the power law fit and either flatten out (see Figure 3a) or drop toward zero (see Figure 3e). Both parameters represent significant changes in event density, with W S2 indicating the transition from the zone of active faulting to a largely undeformed crustal host rock (see Section 4.3 for a comparison with Powers & Jordan, 2010). This transition generally occurs when less than 20% of all boxes used to count aftershocks include events (gray circles in Figures 3 and S4). We use this criterion to define W S2 for all aftershock distributions considered here.   Dokka (1983), Dokka and Travis (1990), Garfunkel (1974), Jachens et al. (2002), Rockwell et al. (2000), and Rubin and Sieh (1997) +0 Atwater and Stock (1998), Critelli and Nilsen (2000), Crowell (1979), Graham et al. (1989), Liu et al. (2010), Matthews (1976), Revenaugh and Reasoner (1997) Gurrola and Rockwell (1996), Hudnut and Sieh (1989), Kirby et al. (2007), and Lutz et al. (2006) New Zealand -0.12-1.83 (field measurements) ∼2.80 (field measurements) Little (1995) and references therein Uncertainties in W S1 and W S2 correspond to size of bin (0.03 km) used in this analysis. For multiple broken fault sections (Hector Mine, Landers, El Mayor Cucapah), W S1 and W S2 are mean values, when possible, and the uncertainties represent the minimum and maximum range of values (see detailed measurements in Figure S3). Field measurements of the Awatere fault from Little et al. (1995) are also indicated (see text for details). a cumulative slip deduced from scaling relations in Scholz and Lawler (2004) (West Napa fault). Initiation age of the West Napa fault is estimated fromthe age of surrounding faults in California. See text for details.  that a single value could be assigned to each aftershock sequence. The variance in the parameters obtained for individual segments is typically small. Table 1 lists W S1 and W S2 we obtained for the seven earthquakes.
Also listed in Table 1 are data from the immature Awatere fault. This right-lateral strike-slip fault is a first order splay of the Alpine fault of New Zealand. It crosses the coast at a sea-cliff that offers an almost complete exposure of the entire fault zone from the fault core through the shear deformation zone. This was mapped by Little (1995), who showed that the shear deformation zone consists of a set of right-lateral and conjugate left-lateral secondary faults that decreased in frequency with distance from the primary fault in the same manner as the aftershocks do in our study. The right-lateral set of the secondary faults is nearly subparallel to the main fault and the left-lateral set about 60° from that. Values of W S1 and W S2 obtained PERRIN ET AL.
10.1029/2020JB020447 6 of 18  Figure S4 for all cases). Blue circles represent the mean hypocenter distances from each fault section within 0.03 km wide bins (see black curve in Figure 2 and Figure S3). Gray circles represent the mean distances of events that are only included in less than 20% of all boxes. The red curve is the best (power law) fit to the blue circles, taking into account the power law singularity on the fault (see Figure S4 for details). The vertical gray dashed lines labeled W S1 and W S2 point out the locations where the power law fit shows the maximum curvature and where the data (blue circles) deviate from the power law fit by dropping or flattening (i.e., the transition between seismically active fault damage zone and the undeformed crustal medium), respectively.
from that data are indicated in Table 1. An exposure at mid-crustal depths of a strike-slip fault in Austria shows just the same orientation of secondary faults (Frost et al., 2009). These observations provide the "ground truth" for the structures upon which the near-fault aftershocks occur.

Fault Parameters
We collect key parameters describing the degree of evolution of the faults that broke during the selected earthquakes. Since faults propagate laterally through time, their structural maturity varies also along strike Perrin, Manighetti, Ampuero, et al., 2016). Thus, it is necessary to use fault parameters that describe the local fault maturity where the rupture occurred. This is particularly true for long faults, such as the San Andreas Fault, for which fault initiation age varies greatly along the fault, from ancient fault sections in Central California (24-29 Ma at Parkfield; e.g., Atwater & Stock, 1998;Liu et al., 2010) to younger fault sections in Southern California (<12 Ma, e.g., Powell & Weldon, 1992;Sims, 1993), and which therefore have different local net displacements. Table 1 presents the fault parameters used in this study. The seven fault sections span a wide range of structural maturity, with various initiation ages (1.1-29 Ma), cumulative displacements (1-315 km) and geological slip rates (0.2-26 mm/yr). Some fault parameters are not fully constrained, especially for the fault associated with the South Napa earthquake. The West Napa fault corresponds to the northern extension of the Calaveras fault, which transfers its slip northwards via the Contra Costa shear zone (Brossy et al., 2010;Catchings et al., 2016;Unruh et al., 2002). We estimated its cumulative slip using a linear fault tip model (Scholz & Lawler, 2004). The West Napa Fault terminates about 35 km north of the South Napa earthquake, and the tip taper for interacting faults is 0.1-0.8 (Scholz & Lawler, 2004). Combining these values give a total slip range of 3.5-28 km ( Table 1). As the West Napa Fault is genetically linked to the Calaveras fault and situated at its northern tip, we assume that it formed after the formation of the main Calaveras fault (about 12 Ma; Wakabayashi, 1999) and before the younger surrounded Greenville and Green Valley faults situated to the east (about 6 Ma ago; Wakabayashi, 1999). The slip rate on the West Napa fault is estimated to exceed 1 mm/yr, based on geomorphic evidence (Wesling & Hanson, 2008) and to be 4 ± 3 mm/yr as derived from geodetic data (d'Alessio, 2005;Field et al., 2015). For clarity, the estimated geological parameters of the West Napa fault are represented as open symbols in Figure 5.

Results
We investigate the relationship between the aftershock distributions, as characterized by W S1 and W S2 , and the various, independently derived fault parameters (Table 1). Our first observation is that W S1 (i.e., the half distance from the calculated fault plane where the aftershock rate saturates) and W S2 (i.e., the full halfwidth of the shear deformation zone) correlate with each other (Figure 4), indicating that these two parameters are not independent and that the shape of the shear deformation zone grows in a self-similar way as In Figure 5, we plot W S1 and W S2 as a function of various fault parameters. When plotted against cumulative fault displacement (Figures 5a and 5b), they both indicate that the width of the shear deformation zone (as characterized by W S1 and W S2 ) decreases with net displacement as a power law with an exponent of −0.54 and −0.39 for W S1 and W S2 , respectively. We note that the field measurements across the Awatere fault (red symbols in Figure 5) are in good agreement with the trends highlighted from seismological data. This supports the interpretation that the shear deformation zone defined by the width of the aftershock zone corresponds to the width of the zone of active secondary faulting. Relations between W S1 and W S2 for the seven earthquakes analyzed in this study. Gray lines indicate two affine functions between W S1 and W S2 that bound our data. Uncertainties are reported in Table 1.
Furthermore, we observe a decrease of W S1 and W S2 with increasing fault initiation age (Figures 5c and 5d) and geological slip rate (Figures 5e and 5f), showing a wider shear deformation zone for younger fault sections and slower slip rates. This is also in good agreement with geological observations across the Awatere fault. We note that the power law exponents in Figure 5 should be used with caution given the occasionally large uncertainties in the fault parameters.  Table 1 and text for details). Power laws are indicated by gray lines. For comparison, red symbols indicate geological surface measurements along the Awatere fault (from Little, 1995). D is cumulative fault displacement, A is initiation age of the fault, SR is slip rate, k is a different constant in each relation.
narrow deformation zone at the scale of hundreds of meters. In contrast, intermediate to immature faults (cold colors in Figure 6) exhibit a wider deformation zone at the kilometer scale where events are more widely distributed within the surrounding medium (Figures 6a and 6b). The distributions remain proportional to one another, indicating that the shape of the shear deformation zone remains constant. Field measurements of the cumulative number of faults across the Awatere fault (black squares, Figure 6a; Little, 1995) show a similar trend and corroborate our seismological observations.

The Nature of the Shear Deformation Zone and the Smoothing of Faults
In this study we show that, for fault displacements greater than 1 km, the width of the shear deformation zone (defined by W S1 and W S2 ) decreases with fault displacement. So, what determines the width of the zone of near-fault aftershocks that define the shear deformation zone? If the fault is smooth, the mainshock would produce a stress shadow centered at the fault and there would be no near-fault aftershocks (Kostrov & Das, 1984). Slip on a rough fault, on the other hand, would result in high stresses in a narrow region adjacent to the fault, within which earthquakes could be generated. This was the explanation used by Powers PERRIN ET AL.
10.1029/2020JB020447 9 of 18 Interpretative cross section describing the structural makeup of immature (blue) and mature (red) faults. As fault structural maturity increases, inner and outer bounds of the shear deformation zone (W S1 and W S2 , respectively) decrease, as expressed by a decreasing width of the fault-normal aftershock distribution. See text for more discussion.
and Jordan (2010) to explain near-fault seismicity during the interseismic period. There they used the static rough fault model of Dieterich and Smith (2009).
Here, we use instead the model of Aslam and Daub (2018) which calculates the stresses resulting from dynamic ruptures propagating on a rough surface. They characterized the fault as a self-affine fractal with Hurst exponent H and fault roughness measured by the RMS of the ratios between the amplitude and the wavelength of the fault surface topography (γ). This is a fairly realistic rendition of the observed topography of faults (e.g., Candela et al., 2012). Their model predicts large changes in the Coulomb Failure Function (CFF = τ + μσ; of more than 20 MPa) within a well-defined narrow region, W nf close to the fault. The modeled receiver faults they assumed are parallel to the primary fault, which is consistent with the orientation of most of the secondary faults observed by Little (1995) and Frost et al. (2009). They find that the width of this near-fault region of high stresses is insensitive to H but decreases with γ. Assuming that our observed shear deformation zone represents the nearfault region of high stresses, this suggests that the observed decrease of W S2 with fault displacement is the result of smoothing of the fault with slip, similarly to observations made during laboratory experiments (Goebel et al., 2017) Aslam and Daub (2018) found that the half-width of the near-fault zone W nf is ∼2.7 and ∼0.9 km for profiles with fault roughness RMS values γ of 0.01 and 0.001, respectively. Comparing these results with the range of W S2 in Figure 5 indicates that more than an order of magnitude of roughness change will be required to explain these observations. This indicates a wear rate far greater than that observed for the roughness of individual fault segments measured in the outcrop scale range for fault displacements in the range 10 −1 to 10 3 m . Similar measurements made on fault segments in the high porosity and friable Navajo and Entrada sandstones of Utah (Dascher-Cousineau et al., 2018) indicate greater smoothing rates but predict fault roughness at our scales orders of magnitude smoother than indicated by our results as interpreted with Aslam and Daub (2018). However, faults are composed of many segments or subfaults at many scales (Ferrill et al., 1999;Klinger, 2010;Manighetti et al., 2015;Scholz, 1998;Wallace, 1989) and of nested subfaults offset from one another (Ben-Zion & Sammis, 2003;de Joussineau & Aydin, 2009;Segall & Pollard, 1980). The length distribution of subfaults follows a power law (Scholz, 1998) and they are offset by jogs that are self-similar (de Joussineau & Aydin, 2009). The combination of the latter relations would produce a fractal topography at a hierarchy higher than that of the roughness of the individual subfault. A fault smoothing relationship based on roughness at this supra-segment level, measured as fault step spacing, shows an approximately inverse linear reduction of fault roughness with fault displacement (Stirling et al., 1996). This is a much higher rate of smoothing than that of individual fault segment surfaces, and is comparable to that we infer from the narrowing of the shear deformation zone. Converting the roughness parameter step frequency of Stirling et al. (1996) to γ using the typical step width 1-2 km (Wesnousky, 1988) and their smoothing curve allows us to recast the two fiducial points of Aslam and Daub (2018) from (W nf , γ) to (W nf , D) coordinates (i.e., γ values of 0.001 and 0.01 corresponding to fault displacements "D" of 100 and 10 km, respectively). These are plotted in Figure 7 (blue circles) where they are compared with our W S2 measurements from Figure 5b (squares). The good agreement supports the hypothesis that narrowing of the shear deformation zone with fault displacement is the consequence of the smoothing of the fault with wear (see also discussion in Section 4.3).

The Effect of Fault Maturity on Earthquake Properties
Once a fault reaches the seismogenic thickness, with a length of ∼10 km and slip ∼100 m, its core and dilatant damage zones are fully formed and further maturation occurs by the fault becoming smoother through the processes of frictional wear. It follows that the characteristics of earthquakes may also evolve with the PERRIN ET AL.   Aslam and Daub (2018) that were calibrated to fault displacement with roughness-displacement relation of Stirling et al. (1996) (blue circles, dashed lines and shaded area). See text for explanation of the latter. Blue circles are for step widths of 1 and 2 km, which are typical step width range (Wesnousky, 1988). smoothing of the fault. It has been suggested, for example, that earthquake stress drop, apparent stress, and radiation efficiency vary with fault maturity (Choy et al., 2006;Hecker et al., 2010;Ross et al., 2018). These are not entirely independent parameters. The radiation efficiency, η R , is given by: where E R is radiated energy, E G is energy dissipated in damage,  Δ is (static) stress drop, σ a is apparent stress, M o is seismic moment, and μ is shear modulus. Here we focus on stress drop, which is the most commonly reported of these properties.
Five of the earthquakes in our study have three or more estimates of stress drop (  Δ ). Following the procedure of Hardebeck (2020), we assume that log 10 (  Δ ) has a normal distribution (Baltay et al., 2011) and calculate its mean and the standard error of the mean, as reported in Table 2. The corresponding values of mean  Δ are plotted versus the cumulative fault slip in Figure 8. The correlation  Δ ∝ D −0.23 , compared to our earlier finding that W S2 ∝ D −0.39 implies that stress drop scales with W S2 and hence, from Figure 7, that stress drop scales with fault roughness. Hecker et al. (2010) measured the maximum slip to length ratio of prehistoric earthquake scarps on intraplate dip-slip faults in the western U.S. They found that this ratio, an indicator of stress drop, tends to decrease with net fault displacement, consistent with our finding.
However, net fault slip usually correlates with slip rate (see Figure 5), so we also find a similar correlation between stress drop and that parameter, as shown in Figure 9a (the relation between stress drop and fault initiation age is shown in Figure S5). This is shown also in log-linear coordinates in Figure 9b to compare it with the form expected from fault healing. The correlation of stress drop with the log of slip rate or its inverse, the recurrence time, is well known (Kanamori & Allen, 1986;Scholz et al., 1986) as is its correlation with strain rate (Hauksson, 2015). The size of the effect seen in Figure 9b is about −2 MPa/decade in slip rate, in rough agreement with the other studies. Scholz et al. (1986) pointed out that this behavior could either reflect a rate effect due to fault healing as predicted by the rate-dependent and PERRIN ET AL.    Table 2. state-dependent friction laws (Dieterich, 1972) or a roughness effect in which friction increases with surface roughness, as indicated by several experimental rock friction studies (Biegel et al., 1992;Byerlee, 1967), or some combination of the two.
Observations of repeating earthquakes in which seismic moment increased with recurrence time promised a separation of these effects (Marone et al., 1995;Vidale et al., 1994). Simulations with one-dimensional spring-slider models with rate-state friction indicated that these data could be explained with fault healing at a rate of 1-3 MPa/decade. However, subsequent studies using a 3D model of a repeating earthquake occurring on a velocity-weakening patch imbedded within a creeping velocity-strengthening fault showed that the recurrence time-moment scaling could occur without any increase in stress drop (Chen & Lapusta, 2009). Observations at Parkfield, California, showed both positive and negative moment-recurrence time scaling of repeating earthquakes (Chen et al., 2010). Both types of scaling can be predicted with the Chen and Lapusta (2009) model by varying the ratio of nucleation to patch size.
In the Chen and Lapusta (2009) model, aseismic slip propagates from the external creeping region into the velocity-weakening patch such that only a fraction of its slip is seismic. The ratio of seismic to aseismic slip in the velocity-weakening patch is sensitive to parameters such as the radius of the velocity-weakening patch and its ratio to the nucleation radius. As a result, it behaves very differently from a simple velocity-weakening block slider model in which the slip is entirely seismic. The recurrence time-moment scaling can be reproduced by varying the patch radius without varying stress drop. In our problem, however, we are considering the fault to be fully coupled seismically so that its slip is entirely seismic. In that case, the simple spring-slider model should be more applicable. If so, the analysis of Beeler et al. (2001) finds that stress drop should increase logarithmically with recurrence time at about 0.5-1.5 MPa per decade, assuming typical laboratory values for the friction parameters a and b. This is in reasonable agreement with the slope in Figure 9b of −2 MPa/decade in slip rate, although, as we said earlier, the effect of fault roughness on stress drop cannot be ruled out.

Strengths, Limitations, and Comparison With Other Studies
In our study, we use aftershock distributions to study how the early postseismic events are reactivating the long-term damage zone around active faults. We define W S2 as the distance from the fault plane where the mean aftershock distributions start to flatten out or drop abruptly (Figures 3 and S4). Either the flattening or the abrupt termination of the number of events with distance away from the fault would indicate that there are no faults or fractures big enough to accommodate a large number of events. Powers and Jordan (2010) observed a flattening of the interseismic event distribution away from the fault which they attributed to background seismicity. Since we are only using 2 months long aftershock sequences, with significantly larger number of events compared to the interseismic period, we do not systematically observe the curve flattening and cannot reliably estimate the off-fault background seismicity level as Powers and Jordan (2010) did for interseismic periods. If significant off-fault seismicity exists during our observational period, then it may have occurred shortly after the mainshock, and went undetected due to the low signal-to-noise ratio PERRIN ET AL.
in the coda of the mainshocks. However, a few of our aftershock distributions show curve flattening (e.g., Parkfield, South Napa, Morgan Hill). These sequences occurred on intermediate to mature faults where creep has been observed (e.g., Harris & Segall, 1987;Schaff et al., 2002;Xu et al., 2018). Our hypothesis is that these aftershock distributions have fewer absolute number of events because they are related to smaller mainshocks (M ∼ 6), and also lower stress drops (see Figure 8). Moreover, the presence of creeping sections might increase the regional background seismicity level.
In this study, we selected aftershocks of M L > 1, thus including events with magnitudes that are below magnitudes of completeness, M c , levels for most aftershock sequences, especially in southern California. Assuming that small events are randomly distributed in the stress field, recovering and including all small events in our analysis would increase the number of events in all boxes. This would result in shifts along the y axis in Figure 6a, but would not change the shape of the curves. It might lead to slight increases in the W S2 values, which therefore should be considered here as lower bounds of our shear deformation zone widths. This is in good agreement with Figure 7, where some measurements corresponding to southern California sequences seem to be under-estimated compared to the expected theoretical models.
Our estimates of the shear deformation zone width might also be affected by the aftershock productivity, which depends on the magnitude of the mainshock. The number of aftershocks does increase with magnitude, but also linearly with fault area (Wetzler et al., 2016;Yamanaka & Shimazaki, 1990), and therefore the density of aftershocks per unit area is constant. This is what we show in Figure S6, where the mean number of events in boxes per km 2 are quite similar and do not show any magnitude dependency (Parkfield and Morgan Hill being slightly off because of the creeping areas that bound the rupture zones, which produce abundant interseismic earthquakes on the fault, contaminating the data at small distances). The apparent correlation between mainshock magnitude and W S2 in our study results rather from the happenstance that the larger earthquakes happen to have occurred on low maturity intraplate faults whereas the smaller ones occurred on mature faults.
As mentioned earlier, Powers and Jordan (2010) estimated the normal distance of the earthquakes away from a single vertical fault plane inferred from the surface trace of the fault. They defined the zone of shear deformation from the near-fault seismicity during the interseismic period and averaged over the entire fault length, possibly neglecting local variations in fault plane orientation and dip, leading to wider zones compared to our results. In addition, geological heterogeneity over long distances along the fault and other effects (such as triggering or regional stress release caused by other earthquakes) may play a role, as pointed out by Hauksson (2010). Finally, our approach of fitting a single fault plane to the local aftershock distribution as we step along the fault assumes that the rupture occurred mainly on one fault strand and that aftershocks are homogeneously distributed in the surrounded medium. While individual secondary faults are distinguishable around mature fault sections (see, for instance, oblique linear faults around Morgan Hill in Figure S3), it is more difficult to identify them around immature fault sections, within the broader zone of deformation, especially near the rupture tips (see the northern part of the El Mayor Cucapah earthquake for instance; Figure S3).
Previous studies have estimated damage zone width or low velocity zones across several faults in California using various seismic and/or geodetic methods (see Yang, 2015, for a review). We compare some measurements with our results in Table S2. First, we note that, different methods result in different estimates of "damage zone" widths. Second, all studies using seismic fault zone trapped waves (FZTW) techniques find a similar low velocity zone that is 100-200 m wide, regardless the fault structural maturity. The FZTW derived width correspond to the constant width of dilatant damage zones around large faults, already described by geological studies (e.g., Faulkner et al., 2011;Savage & Brodsky, 2011; also called the inner damage zone in Perrin, Manighetti, Ampuero, et al., 2016). However, other seismic and geodetic studies (sometimes used together; see Cochran et al., 2009), find greater damage zone widths, that are in general agreement with our results. This discrepancy between results from different methods was also pointed out by Spudich and Olsen (2001). We tested the relations found in Figures 5a and 5b to deduce the size of the shear damage zone of the San Jacinto fault, for which both a narrow damage zone using FZTW methods (e.g., Lewis et al., 2005) and wider damage zone using tomography (e.g., Allam et al., 2014) have been described. Assuming a cumulative slip of 24 km (Sharp, 1967), we obtain a half-width of the shear deformation zone ranging from 0.2 (i.e., W S1 ) to 1.5 km (i.e., W S2 ), which is in better agreement with tomographic studies (half-width of 1.5-2.5 km; Allam et al., 2014). Thus, we believe that the studies using FZTW are mainly sensitive to the most damaged inner part of the fault zone, also called the dilatant damage zone, which we do not seem to be able to resolve with our aftershock data. The other seismic and geodetic studies found damage zone widths that are between 2(W S1 ) and 2(W S2 ), which we define as the inner and outer bounds of the shear deformation zone.
A recent geodetic study from Scott et al. (2018) has shown that the dilatant damage zone is equivalent to slightly wider than the shear deformation zone. These results come from surface measurements (i.e., LiDAR) along the active fault zone responsible for the 2016 Kumamoto earthquake. It is possible that extensional strains are greater at the surface than at depth. As cracks close down rapidly under pressure, dilatant strains in the dilatant damage zone will decrease more rapidly with depth (Yoshioka & Scholz, 1989) than the shear strains. This may explain why their results from surface measurements differ from ours measured at depth.
In our study, we average the near-fault aftershock distributions assuming that the degree of fault maturity is locally constant along each broken fault sections we investigate. However, it has been shown that the fault maturity can vary along strike and this can affect the heterogeneity of fracture density around the fault core (Ostermeijer et al., 2020), the distribution of secondary faults at greater distances away from the fault  and finally the behavior of earthquakes (e.g., Cappa et al., 2014;Huang, 2018;Perrin, Manighetti, Ampuero, et al., 2016). Consequently, for long or multiple-fault earthquake ruptures, the fault-normal distribution of aftershocks could change along strike, and present a shear deformation zone that widens toward the most immature parts of the rupture. Such a behavior would be in good agreement with mapped faults at the surface showing that the off-fault damage zone widens in the direction of long-term fault propagation (e.g., Manighetti et al., 2001;. Future work is needed to relate such observations to the occurrence of seismicity.

Conclusion
Our study presents strong correlations between independent data sets that describe the near-fault distribution of aftershocks following large earthquakes and geological fault parameters such as cumulative displacement, initiation age and slip rate. We find that for large faults, defined as those that have ruptured the entire brittle crust, the zone of active shear deformation narrows as a power law with cumulative fault displacement, hence with fault maturity. This result is predicted by a dynamic rough fault model combined with an empirical roughness/displacement relation, indicating that the narrowing of the shear deformation zone is the result of fault smoothing with cumulative displacement. Earthquake stress drop also decreases with fault displacement and hence fault smoothness or slip rate.