Spatial Relationships Between Coseismic Slip, Aseismic Afterslip, and On-Fault Aftershock Density in Continental Earthquakes

Damaging aftershock sequences often exhibit considerable spatio-temporal complexity. The stress changes associated with coseismic slip and aseismic afterslip are commonly proposed to drive aftershock sequences, but few systematic studies exist and do not always support strong, universal driving relationships. To investigate the roles that these two sources of stress changes may play in driving aftershocks

Coseismic elasto-static Coulomb stress change is the most established candidate mechanism for triggering aftershocks (King et al., 1994;Stein, 1999).However, aftershock models based on coseismic Coulomb stress change do not always perform well (e.g., Cattania et al., 2018;Hardebeck, 2021;Woessner et al., 2011), possibly due to some methodological shortcomings (Cattania et al., 2014;Hainzl et al., 2010;Hardebeck, 2021;Mancini et al., 2019Mancini et al., , 2020) and because of the potential that mechanisms such as afterslip may also influence aftershock triggering.On fault, it has been proposed that aftershocks commonly occur where no to little coseismic slip has occurred (Dreger, 1997;Mendoza & Hartzell, 1988;Woessner et al., 2006), and that aftershocks may effectively fill in some regions that have not already slipped.Such claims therefore imply that coseismic slip is an important spatial control on where on-fault aftershocks occur.
Predominantly case-study-based evidence has led to inferences that afterslip may generally drive aftershocks (e.g., Avouac, 2015), yet the few extant efforts to systematically test the relationship do not support this.For example, Lange et al. (2014) compared afterslip and aftershock distributions for three large subduction earthquakes and concluded that afterslip was not the sole driver of the aftershocks, indicating that the stress change due to the mainshock and loading from the interseismic period are (as, or possibly more) important factors.Cattania et al. (2015) incorporated afterslip into Coulomb Rate-and-State models following the 2004 M w 6.0 Parkfield and the 2011 M w 9.1 Tohoku earthquakes, finding that it was also not the principal driver of the aftershocks and that secondary triggering (i.e., subsequent generations of aftershocks being triggered by the evolving stress state caused by earlier aftershocks) may play a more important role.In a systematic analysis of 41 mainshocks, Churchill et al. (2022b) showed that relative afterslip moment  (  aslip rel ) does not correlate discernibly with key characteristics of aftershock sequences, including the productivity of M w ≥ 4.5 aftershocks.To reconcile case study evidence that supports strong links between afterslip and aftershocks with systematic testing that appears to contradict this, Churchill et al. (2022b) asked: is the role of afterslip in driving aftershock sequences simply spatio-temporal in nature and/or specific to certain settings only?
In this study, we systematically explore the spatial relationships between distributions of coseismic slip, afterslip, and on-fault aftershock density through time, for seven M w 6.0-7.6 continental-setting mainshocks (in California & Northern Mexico, Turkey, and Italy, outlined in Section 3).Whilst well-studied individually, these continental-setting earthquakes have not been subject to systematic afterslip-aftershock analysis, such as in the subduction-setting study by Lange et al. (2014).For these earthquakes, high-quality coseismic slip models, afterslip models, and regional aftershock data are available, and the earthquakes constitute several key case studies analyzed by Churchill et al. (2022aChurchill et al. ( , 2022b) (e.g., in terms of outlying relative afterslip moment or relative aftershock productivity).For simplicity, we test relationships on a single best-fitting plane to the fault geometry and assume that aftershocks close to the fault plane (within a cut-off distance) occur on it rather than in the surrounding volume.Analyzing on-fault distributions allows us to test in a simple manner the spatial relationships between distributions of slip and aftershock density without considering stress in 3D, thus avoiding the complexity and 10.1029/2023JB027168 3 of 26 uncertainty inherent to 3D (tensorial) analyses such as Coulomb stress change modeling, which is uncertain in the near-field (specifically where we wish to test), and would require numerous additional assumptions (e.g., receiver fault orientations) to be made (e.g., Cattania et al., 2014;Hainzl et al., 2010;Mancini et al., 2019).We also test whether slip gradients (which might be viewed as a proxy for new stress concentrations) correlate spatially with aftershock distributions.This study aims to provide pragmatic, empirical constraints on where on-fault aftershocks could be expected given distributions of coseismic slip and afterslip (considering the uncertainties and resolution of data and models available).We also discuss the relative importance of coseismic slip and afterslip in controlling the spatio-temporal distribution of on-fault aftershocks and what our findings may imply about the fault zone.

Hypotheses
In this section, we describe our expectations for each of the tests undertaken in this study, based on theoretical considerations and previous observations.We test the spatial relationships between: (a) coseismic slip and afterslip, (b) coseismic slip and on-fault aftershock density, (c) afterslip and on-fault aftershock density, and (d) total slip (the sum of coseismic slip and afterslip) and on-fault aftershock density.In addition to this, we substitute the first slip variable in each of these tests with its respective slip gradient, to explore whether slip gradient provides more discernible relationships with the second variable.We assume that slip gradient can be used as a proxy for new (shear or total Coulomb) stress concentrations given the link between slip distribution and stress drop on the fault (e.g., Kaneko & Shearer, 2014, 2015;Noda et al., 2013).
First, we hypothesize, as per a very common assumption, that close to the mainshock rupture, the spatial distributions of coseismic slip and afterslip should be anti-correlated.Mechanically, one could argue that areas which have undergone significant coseismic slip have undergone a significant stress drop (Jaeger et al., 2009) and do not require significant afterslip.It is also generally understood that afterslip is confined to velocity-strengthening regions within the fault zone, whereas coseismic slip occurs in velocity-weakening regions (Bürgmann, 2018;Dieterich, 1979;Marone et al., 1991;Perfettini & Avouac, 2007;Perfettini et al., 2018), and that afterslip is triggered by the stress concentrations adjacent to coseismic slip patches (Bürgmann, 2018;Perfettini et al., 2018).Therefore, coseismic slip and afterslip should not overlap significantly, particularly when considering simple layered models of the crust (e.g., Hillers et al., 2006;Marone et al., 1991).This expected relationship has been supported by some previous observations (e.g., following the 2009 M w 6.3 L'Aquila earthquake: D 'Agostino et al., 2012;Gualandi et al., 2014).However, overlap of afterslip and coseismic slip distributions may be possible, through phenomena such as dynamic slip weakening (Noda et al., 2013) and slip reversals (e.g., Jiang et al., 2021a), and apparent overlap in inverted slip models may be due to unresolved frictional heterogeneity and/or conditional slip stability within the fault zone (Boatwright & Cocco, 1996;Bürgmann, 2018;Bürgmann et al., 2002;Helmstetter & Shaw, 2009;Scholz, 1998).
Second, we hypothesize that the spatial distributions of coseismic slip and aftershocks should also be anti-correlated.Again, as stress drops are large where changes in displacement (coseismic slip) are large (Jaeger et al., 2009), one could assume that high slip regions do not require further slip to be accommodated by aftershocks.However, unlike with coseismic slip and afterslip, coseismic slip and aftershocks are not confined to separate parts of the fault according to the arguments above, with both expected to occur within velocity-weakening material.Several studies testing the spatial relationships between coseismic slip (and associated factors such as moment release and stress drop) and aftershock density have indicated that aftershock density is either entirely anti-correlated with coseismic slip, and/or correlates with the edges of coseismic slip (Dreger, 1997;Mendoza & Hartzell, 1988;Wetzler et al., 2018;Woessner et al., 2006).Despite these points, however, analyses by Das and Henry (2003) suggested that there is no universal relationship between spatial distributions of coseismic slip and aftershocks across different case studies.Interestingly, Neo et al. (2021) proposed that mainshock rupture area may serve as a rough proxy for the aftershock zone (due their to spatial overlap at low resolutions), but also indicated that the ratio of aftershock zone area to mainshock rupture area is typically greater than one (and up to ∼5, indicating that a significant proportion of aftershocks may occur outside the mainshock rupture zone).In Section 5.5, we discuss how different choices of domain over which to test relationships may result in contradictory findings in some cases (essentially due to collider bias).
Third, we hypothesize that the spatial distributions of afterslip and aftershocks should be correlated.Whilst the two phenomena are expected to occur in different frictional stability regimes, existing observations indicate that afterslip and aftershocks migrate together with time, particularly shortly after the mainshock (e.g., Feng et al., 2020;H. Huang et al., 2017;Jiang et al., 2021a;Peng & Zhao, 2009).This is also supported by some models (e.g., Kato, 2007;Perfettini et al., 2018), which suggest that afterslip expands as a halo around the rupture with time, triggering aftershocks on isolated velocity-weakening fault asperities.Additionally, some studies claim that afterslip and aftershocks cluster together at the edges of coseismic slip (e.g., D'Agostino et al., 2012;Ozawa et al., 2004;Perfettini et al., 2010).This expected relationship is theoretically and observationally consistent with the two previous expected relationships (i.e., afterslip and aftershocks are spatially correlated, and both are anti-correlated with coseismic slip).
Given our contrasting expectations regarding the expected spatial relationships between distributions of coseismic slip and aftershocks, and afterslip and aftershocks, it is more difficult to form a single expectation for the spatial relationship between total cumulative slip and aftershocks.As coseismic moment is generally, but not always, greater than afterslip moment (Churchill et al., 2022a), one could posit that coseismic slip is the more important factor, thus the expected spatial relationship between total slip and aftershock density should be anti-correlated, albeit perhaps more weakly.We test this relationship specifically to explore whether it may be a useful tool for understanding the distribution of on-fault aftershocks, rather than to gain any specific mechanical insight.

Case Studies and Data
We provide background to the seven continental earthquake mainshocks analyzed in this study and outline the corresponding sources of data (summarized in Table 1).Estimates of relative afterslip moment (   aslip rel , afterslip moment divided by coseismic moment) correspond to those calculated by Churchill et al. (2022a) or approximated subsequently, in which values of ∼10%-30% are considered typical.These values are subject to some data and methodological uncertainties, but real differences (from earthquake-to-earthquake) are likely driven by the available fault area with frictional characteristics favoring afterslip (Rolandone et al., 2018).Estimates of relative aftershock productivity correspond to those calculated by Churchill et al. (2022b) or approximated subsequently using similar methods.Churchill et al. (2022b) defined relative aftershock productivity (n rel ) by the number of aftershocks observed (above a given magnitude) divided by an expectation based on the Utsu-Seki productivity scaling law (Utsu, 1970).They found values in the range ∼0.1-10 (where values ∼1 are considered typical) and that n rel does not correlate with   aslip rel .Similar variability in aftershock productivity has been observed by other studies, which are also unable to identify a dominant driving factor (e.g., Dascher-Cousineau et al., 2020;Marsan & Helmstetter, 2017).
We compiled most afterslip and coseismic slip models by contacting the original authors, except where resources could be found in original Supporting Information and through the Earthquake Source Model Database (Mai & Thingbaijam, 2014).Although slip models are increasingly available open access (e.g., Mai & Thingbaijam, 2014), data on afterslip distributions were often only available by contacting authors directly.We chose the afterslip models analyzed here based on factors such as the depth of the analysis and exploration of modeling space in ) , relative aftershock productivity (n rel ), a conservative (common) estimate of completeness (M c ), and the approximate number of (≥M c ) on-fault aftershocks in the 18 months following the mainshock (N on-fault , discussed in Section 4.1) are given for the earthquakes studied.For cases such as the Ridgecrest and El Major Cucapah earthquakes, in-faultzone may be a more appropriate term than on-fault.

Table 1
A Summary of Coseismic Slip, Afterslip, and Aftershock Data Analyzed, Ordered by Region the original study, quality of data, number of time-steps, and crucially, availability of a digital model.Afterslip distribution models are geodetically constrained, given the aseismic nature of afterslip, but often make use of a seismologically constrained fault geometry.Conversely, coseismic slip distribution may be resolved using direct seismological and geodetic constraints.Where possible, we analyze the coseismic slip model produced by the afterslip model study, but where this is not possible, we use a suitable coseismic slip model from the SRCMOD database (Mai & Thingbaijam, 2014).
There are substantial uncertainties in slip distribution models, including uncertainty and non-uniqueness that is inherent to numerical inversions (Menke, 2018;Scales & Tenorio, 2001).Limitations to the scales of resolvable slip (Lohman & Simons, 2005) is mostly a function of geodetic data quality and density (in slip models, the slip-resolution does not equal the model cell-size, which is generally set to a much finer scale and then smoothed to some degree).Additional assumptions, such as constraints on topography and assumed rheological structure, have significant effects on the final slip model (e.g., Hearn & Bürgmann, 2005;Marchandon et al., 2021;Sun & Wang, 2015).It is also possible that signal from early afterslip may be inadvertently and erroneously included within the coseismic slip model and/or some coseismic slip may be included in the very first afterslip time-step.Some studies specifically attempt to account for this (Churchill et al., 2022a), but it exists as a minor source of uncertainty in models, particularly those with infrequent time-steps of their afterslip models.Churchill et al. (2022a) discuss afterslip model uncertainties in more detail.We present a summary of models used in Table 1.We use regional seismic catalogs to investigate the spatial distributions of aftershock sequences.Regional catalogs generally have better location accuracy than global catalogs, as well as a greater data density due to lower completeness magnitudes, M c .The specific catalogs used are described in the following sections.We consider the following: (a) as M c can vary between catalogs and through space and time within the same catalog, we determine a conservative value of M c for each aftershock sequence, and establish a common M c when comparing different sequences, (b) whilst regional catalogs typically contain a mixture of magnitude scales including local magnitudes (M L ), we assume that these are sufficiently calibrated to M w between ∼ M2 and M6 (the approximate magnitude range they were designed for, Richter, 1935;Deichmann, 2006), (c) we note that location uncertainties are inherent to seismic catalogs and are typically worse in the vertical direction and with increasing depth (Kagan, 2003;Waldhauser & Schaff, 2008), however, the limits of acceptable location uncertainties are effectively constrained by the length-scales of features (i.e., individual slip patches) in the coseismic slip and afterslip models which we compare to, and (d) we analyze aftershocks up to 18 months (maximum) after the mainshock to ensure that other potential sources of seismicity do not become dominant in our study regions (given their temporal decay).Additional case-specific points are discussed in the following sections and data limitations are discussed at greater length within the Methods and Discussion sections (Sections 4.1 and 6).For context, Figure 1 shows a map for each earthquake, highlighting fault model geometry and the seismicity in the 18 months that follow the mainshock.
The 2014 M w 6.0 South Napa earthquake was followed by a moderate-to-high relative afterslip moment (   aslip rel ∼ 30%-40%: Floyd et al., 2016;Premus et al., 2022).These models suggest that afterslip occurred as several patches, most of which were shallow (0-5 km), although limited afterslip also occurred at greater depths (5-15 km).Afterslip was rapid (Wei et al., 2015), and likely >95% complete within a year compared to only ∼70% of afterslip within a year following Parkfield (Lienkaemper et al., 2016).The (hosting) West Napa fault  Van, and (g) 2009 M w 6.3 L'Aquila earthquakes.Panels show the (above M c ) seismicity in the 18 months following each mainshock (catalogs are given in Table 1), with color indicating time since the mainshock and size indicating magnitude.The afterslip model fault geometry (see Table 1) is given in black in each case, except for in panel (c), where we show the foreshock fault plane as well, and panel (d), where we include both the (spatially-separated) coseimsic slip model too, for context.exists within the broader Contra-Costa shear zone, and is more structurally complex than the Parkfield section of the San Andreas fault (Brocher et al., 2015).We analyze coseismic and afterslip models (12 time-steps between days 1 and 67) by Floyd et al. (2016).
Both earthquakes occurred within the authoritative region of the Northern California Earthquake Data Center (NCEDC, 2014), hence we analyze aftershock data from this catalog.The catalog has suitably low hypocentral location uncertainty, ranging from hundreds of meters to a few kilometers at worst (Waldhauser & Schaff, 2008).Preliminary analysis of seismicity in the approximate space-time windows of the respective aftershock sequences (within 18 months and within 50 km of the hypocenter) indicates low completeness thresholds.We analyze M c using three methods (throughout this study): maximum curvature, a goodness-of-fit test (Wiemer & Wyss, 2000), and a b-value stability test (Cao & Gao, 2002), which indicate that following the Parkfield earthquake: M c = 1.1, 1.0, and 1.2, respectively, and following the South Napa earthquake: M c = 1.2, 1.5, and 1.5, respectively.We use the most-conservative common estimate of M c = 1.5.The Parkfield aftershock sequence had an average-to-high productivity (Peng & Zhao, 2009;Shcherbakov et al., 2006), with n rel ≈ 1.8 (Churchill et al., 2022b), and links with afterslip have been proposed (e.g., co-migration in the first few days: Peng & Zhao, 2009;Jiang et al., 2021a).Conversely, the South Napa aftershock sequence had a low productivity, particularly on-fault (Hardebeck & Shelly, 2016;Llenos & Michael, 2017;Toda & Stein, 2015), which has been attributed to the spatially-limited distribution of afterslip and a lack of stick-slip patches to facilitate aftershocks (Hardebeck & Shelly, 2016;Premus et al., 2022).

Southern California and Northern Mexico
The 2019 M w 7.1 Ridgecrest earthquake and its M w 6.4 foreshock 2 days prior were followed by a low-to-moderate relative afterslip moment (   aslip rel ∼ 10% within 6 months: Yue et al., 2021).The M w 7.1 right-lateral strike-slip earthquake occurred in a heavily fractured zone (Barnhart et al., 2019), on a plane near-orthogonal to that of the (left-lateral) foreshock, resulting in some complexity in coseismic slip models (e.g., Feng et al., 2020;Jin & Fialko, 2020;Pollitz et al., 2020;Qiu et al., 2020;Ross et al., 2019;Yue et al., 2021).Afterslip occurred on both planes, between the surface and depths of approximately 30 km (Pollitz et al., 2022;Qiu et al., 2020;Yue et al., 2021).As multiple afterslip time-steps are not available from a single afterslip model, we analyze the coseismic and afterslip model (at 5 months) by Yue et al. (2021), and a currently-unpublished afterslip model at 2 years (a follow-up to the study by K. Wang and Bürgmann (2020)), which share a similar fault geometry.We remove the foreshock fault segment (shown for context in Figure 1c) from these models for projection and analysis.
The 2010 M w 7.2 El Mayor Cucapah was followed by a moderate-to-high relative afterslip moment (   aslip rel ∼ 20%-70%: Rollins et al., 2015;Tang et al., 2020).Lower estimates of   aslip rel also exist but are based on time-scales of weeks not several years (Gonzalez-Ortega et al., 2014;Ross et al., 2017).The fault geometry in the region is quite complex: the main rupture had a primarily strike-slip focal mechanism (M.-H.Huang et al., 2017;USGS, 2017;Wei et al., 2011), but involved two major northwest-southeast striking fault segments that dip in opposite directions at dip angles of approximately 60 and 75° (Hauksson et al., 2011;Rollins et al., 2015;Wei et al., 2011), highlighted in Figure 1d.Afterslip mostly occurred beneath the coseismic rupture on both of these major fault segments (Hines & Hetland, 2016;Rollins et al., 2015;Tang et al., 2020).A smaller, shallower fault cutting across these two segments is believed to have hosted an initial, normal mechanism sub-event (M.-H.Huang et al., 2017;Wei et al., 2011).In our study, we analyze the coseismic slip model by Wei et al. (2011) and the afterslip model (at 3 years) by Rollins et al. (2015), who uses down-dip and along-strike extensions of the geometry by Wei et al. (2011).We remove two smaller sub faults which cross-cut the main two segments from these models for projection and analysis.
For these earthquakes, we analyze data from the California Integrated (CI) catalog by the Southern California Seismic Network and Southern California Earthquake Data Center (SCEDC, 2013).Whilst the Ridgecrest earthquake is well within this catalog's authoritative region, the El Mayor Cucapah earthquake is close to the edge, where coverage begins to drop off.This catalog's location uncertainties are thought to be reasonably low: consistently less than 3 km post-1970s and typically much better than this (Hutton et al., 2010).Preliminary analysis of seismicity in the approximate space-time windows of the respective aftershock sequences (within 18 months and within 100 km of the hypocenters) indicates that following the Ridgecrest earthquake: M c = 1.1, 1.2, and 1.9 (by different methods, respectively), and following the El Mayor Cucapah earthquake: M c = 1.5, 2.2, 2.8 (by different methods, respectively).By visual inspection, the M c by b-value stability estimate for the El Mayor Cucapah aftershock sequence appears overly conservative, and in the interest of maximizing available data to analyze, we use M c = 2.2 as the common M c for analyzing these two sequences.Both aftershock sequences were moderately productive (Feng et al., 2020;Jin & Fialko, 2020;Kroll et al., 2013;Pollitz et al., 2022;Shelly, 2020), with n rel ≈ 1.5 following the El Mayor Cucapah earthquake (Churchill et al., 2022b).In both cases, there have been suggestions that early afterslip may have influenced aftershock distributions (e.g., Feng et al., 2020;Ross et al., 2017).
The 2019 M w 7.1 Ridgecrest and 2010 M w 7.2 El Mayor Cucapah earthquakes present some challenges to our analysis.Given the relative geometric complexity of these earthquakes (see Figure 1), the term in-faultzone may be more appropriate than on-fault when we consider their aftershocks (in Section 4), however, we use the term on-fault for consistency.In Section 4.1, we outline how we only select aftershocks close to the main fault plane for analysis.Given the geometric complexity of these two events, it is possible that some aftershocks associated with other close-by, minor fault strands may be included in our analysis.This limitation is difficult to avoid when analyzing data at such resolutions with such spatial uncertainties, but we argue that its effect will be relatively small given that we only consider up to a short cut-off (inclusion) distance from the main fault (detailed in Section 4.1).We also assume that every earthquake that follows the M w 7.1 Ridgecrest mainshock (near its fault plane) is an aftershock of it, and not the M w 6.4 foreshock.We argue that this is a reasonable assumption, given the much greater size of the mainshock than the foreshock, the time already elapsed between the two events, and the generally much closer proximity of the mainshock fault plane to the aftershocks we are analyzing.
For these earthquakes, we analyze data sourced from the International Seismological Center (ISC, 2022).For the two Turkish events, we use the Kandilli Observatory And Earthquake Research Institute (KOERI) catalog (KOERI, 1971), and for the L'Aquila earthquake, we use the Istituto Nazionale di Geofisica e Vulcanologia (INGV) catalog (INGV, 2005).The location uncertainties for the Izmit aftershocks are generally greater than for the Van aftershocks, given that the KOERI network had fewer stations in 1999 (Çıvgın & Scordilis, 2019).Additionally, aftershocks of the Izmit earthquake appear to be at specific depths (default/likely limited by resolution) and aftershocks are missing in the first 90 min.However, analysis of the KOERI catalog is still meaningful following the Izmit earthquake, particularly because the Izmit rupture was so elongate that analysis in the horizontal/along-strike direction will be less affected.Location uncertainties following the L'Aquila earthquake (INGV catalog) are relatively low (up to a few kilometers), with few default depths assumed (Scudero et al., 2021).Preliminary analysis of seismicity in the approximate space-time windows of the respective aftershock sequences (within 18 months and within 150 km of the hypocenter following the Izmit earthquake, and within 100 km following the Van and L'Aquila earthquakes) indicates that following the Izmit earthquake M c = 2.7 (by all methods), following the Van earthquake M c = 2.6, 2.5, 2.6 (by different methods, respectively), and following the L'Aquila earthquake: M c = 1.8, 1.9, 2.0 (by different methods, respectively).We assume the conservative common M c estimate of 2.7.
The Izmit aftershock sequence had a low productivity (n rel ≈ 0.5) (Churchill et al., 2022b), particularly on-fault (Bouchon & Karabulut, 2008;Ozalaybey et al., 2002).This has been linked to the super shear nature of the rupture -possibly because more of the fault was activated in the first place, or because dynamic (triggering) stresses were higher during the rupture (Bouchon & Karabulut, 2008).Conversely, the Van and L'Aquila aftershock sequences were highly productive (Altiner et al., 2013;Chiarabba et al., 2009;Chiaraluce et al., 2011), with n rel ≈ 3.4 and 5.6, respectively (Churchill et al., 2022b).We have not found discussions of strong links between afterslip and aftershocks following the Van earthquake in the literature, but such links are well-established following the L'Aquila earthquake, with strong spatial and temporal correlations reported by several studies (D'Agostino et al., 2012;Gualandi et al., 2014;Yano et al., 2014).

Data Projection and Regridding
We project slip models and aftershock locations orthogonally onto planes in three-dimensional space.We first convert all data to local Cartesian co-ordinate systems: we use a geodesic measure (WGS-84 ellipsoid) which is accurate to round-off (Karney, 2013) and use a local origin-point in each case, so that all conversions are over short distances.We then find the best-fitting plane for each case study: the best-fitting plane is a single plane that best fits the fault geometry of the given afterslip model (determined by least squares regression).We use the afterslip model geometries because these are the same as (or a spatially extended versions of) the coseismic slip model geometries for the Parkfield, South Napa, Ridgecrest and El Mayor Cucapah earthquakes.For the Izmit, Van, and L'Aquila earthquakes, the afterslip and coseismic slip model geometries are sufficiently similar that this is a reasonable approximation.In the case of the Ridgecrest earthquake, the orthogonal foreshock plane is removed from coseismic and afterslip models, and in the case of the El Mayor Cucapah earthquake, the small cross-cutting fault segments are also removed from models.Residual distances between the best-fitting plane and afterslip and coseismic slip models are generally small (<2 km in most cases).However, due to geometric complexity in the El Mayor Cucapah earthquake case, residual distances between the parts of the original fault models and the best-fitting plane may be up to several km, although these points tend to be far along-strike and down-dip parts of the fault geometry.
The use of a best-fitting plane allows us to conduct simple spatial relationship analyses between slip distributions and aftershock density without considering three dimensional stress changes.This approach also allows us to regrid projected slip data and estimate aftershock density over regularly spaced grids, allowing for simple statistical comparisons of variables.Using homogeneous cell-size and a regular grid also allows us to account for overlapping fault segments when projecting onto the best-fitting plane, and effectively helps to bring slip models closer to a similar level of smoothing (and means that calculated gradients are more comparable).This does not homogenize the level of smoothing in slip models, which is decided within the original modeling.
We consider the distribution of distances between aftershock locations and the best-fitting plane to establish a maximum cut-off distance for aftershocks to be included as on-fault.Figure 2 shows the distance distributions between aftershocks and the best-fitting plane in each case: a narrow distribution (e.g., for the Parkfield earthquake) suggests that a significant proportion of aftershocks occurred on-fault, whilst a broad distribution (e.g., for the Izmit earthquake) could indicate a significant proportion of off-fault aftershocks, a more complex fault geometry (meaning a less representative best-fitting plane), and/or uncertainties in the seismic locations themselves.Given the approximate location uncertainties in both the seismic catalogs and in the fault model geometries (including further spatial uncertainties that arise due to projection), we assume that aftershocks that occur close to the best-fitting plane, actually occur on it.We use a maximum cut-off distance of 2.5 km from the plane for the (M w 6.0) Parkfield and South Napa earthquakes, and 5 km for the remaining earthquakes, which are bigger in size (and in the Ridgecrest and El Mayor Cucapah cases, involve complex geometries).Specifically regarding the 2019 M w 7.1 Ridgecrest and 2010 M w 7.2 El Mayor Cucapah earthquakes: whilst Figure 2 indicates that many aftershocks are more than 5 km from the best-fitting plane in either case, we argue that it is reasonable to exclude these aftershocks from our analysis, as many of these aftershocks are not associated with the mainshock plane (as Figure 1 shows, many aftershocks are associated with the foreshock plane in the Ridgecrest case, and many occur along strike in the El Mayor Cucapah case).This therefore suggests that our methodology is appropriate.Maximum aftershock projection distances of up to 5 km are common in the literature (e.g., D' Agostino et al., 2012;Jiang & Lapusta, 2016;Mendoza & Hartzell, 1988;Yano et al., 2014;Yue et al., 2021).Subsequent references to aftershocks in this study specifically refer to assumed on-fault aftershocks (i.e., those which are projected onto the best-fitting plane from within the maximum cut-off distance).
Data are projected orthogonally onto the best-fitting plane and regridded.Figure 3 show examples for coseismic slip, afterslip, and aftershocks in the Parkfield case study (the equivalent Figure is given for the other six earthquakes studied in Figures S1-S6 of Supporting Information S1).In most instances, regridding is achieved by the linear interpolation of data on the plane.However, for the South Napa and Ridgecrest model geometries (which include overlapping segments), we interpolate the projected slip contributions from each overlapping segment onto the grid and then sum these contributions to conserve the slip and moment of the original model.We also assume that no coseismic slip or afterslip occurs outside the given coseismic and afterslip models, even if slip models have non-zero slip near their edges.We do this because additional smoothing may not conserve moment or necessarily be physically meaningful, as slip models can contain coseismic slip or afterslip which terminates relatively sharply at certain boundaries or asperities.We therefore choose to analyze the compiled slip models as presented by the original publications.For all cases except the Izmit earthquake, we use a 2 by 2 km gridding, which retains the major features of each slip model and is a reasonable spatial binning size for aftershocks given their typical location uncertainty.However, in the Izmit case, we use a 5 by 5 km binning due to depth uncertainty The distance-distribution between aftershocks and the best-fitting plane for each mainshock.These histograms include all events up to 18 months after the mainshock, above the common M c used for analysis: 1.5 for the Parkfield and South Napa earthquakes (b), 2.2 for the Ridgecrest and El Mayor Cucapah earthquakes (b), and 2.7 for the Izmit, Van, and L'Aquila earthquakes (c).As the best-fitting plane is extended through the entire initial selection area, these selections may include some along-strike events away from coseismic rupture and afterslip.In each case, the projection distance cut-off is given by a vertical line. in the aftershock data (as the rupture is very elongate, this still allows us to test relationships reasonably well, albeit at a lower resolution).In short, our methods result in projected and regridded slip models which retain the major features of the original slip models and conserve their moment (e.g., Figures 3a and 3c).

Slip Gradients and Aftershock Density
We calculate the gradient magnitude field across our regridded slip models.We use a scalar measure of gradient, so that it can be easily compared with other scalar fields (e.g., afterslip, coseismic slip, aftershock density).At each point in a slip model, we calculate the along-strike and downdip components of the gradient vector (∇D, where D is slip) using a finite difference approach.The gradient magnitude (‖∇D‖) is then simply given by: Figures 3b and 3d show illustrative results for the Parkfield earthquake.
We use a Gaussian kernel approach to estimate aftershock densities on the best-fitting plane.Following the projection (and binning) of aftershock locations onto the best-fitting plane, we make a 2D density estimate using a Gaussian kernel approach (Scott, 2015;Silverman, 2018), which is a relatively well-established means of estimating aftershock densities (e.g., Bai et al., 2022Helmstetter & Werner, 2014;Helmstetter et al., 2006;).This method smooths the raw aftershock density values, which may mitigate the impact of location errors on our analyses to some extent, but is sensitive to the chosen bandwidth (Scott, 2015).The units of aftershock density is normalized count per unit area (i.e., this method does not consider moment, only aftershock number), for a given time-step.An example of a Gaussian kernel aftershock density estimate is given in Figure 3f.

Statistical Analysis
We consider two different methods for designating the spatial domain over which relationships are tested in this study.We explore spatial relationships between (a) coseismic slip and afterslip, (b) coseismic slip and aftershock density, (c) afterslip and aftershock density, and (d) total cumulative slip and aftershock density (and each of these with the respective slip gradient substituted for the first variable) over the spatial grids described in the previous section, through time.In our main tests, we consider all grid cells which contain non-zero values (i.e., if either of the tested variables is non-zero in a given cell, that cell is included in the analysis).Therefore, these main tests effectively analyze spatial relationships across the whole plane: a positive correlation implies that the two tested variables occur in similar regions with correlated magnitudes, whereas an anti-correlation implies that the two tested variables either occur in similar regions with anti-correlated magnitudes, and/or occur in entirely different regions on the fault.Results of our main tests are presented in Sections 5. 1-5.4.
In an additional set of tests, we explore spatial relationships over a reduced spatial domain, only considering cells that contain at least 10% of the maximum value of the first test variable (at each given time step).For example, when exploring the relationship between coseismic slip and afterslip, we would analyze the relationship where there has been significant coseismic slip.Because these tests only examine where the first variable is significant: a positive correlation implies the two variables have spatially correlated magnitudes in this region and an anti-correlation implies spatially anti-correlated magnitudes in this region.We explore whether these tests produce any results which are substantially different from those produced by the main tests and comment on this.
Results of additional tests are presented in Section 5.5.
Throughout all analyses, we use Spearman's rank correlation coefficient as a simple and intuitive test for monotonic relationships across cell values (Dodge, 2008).A rank-based approach helps mitigate the impact of outlying data-points on inferred relationships and is less sensitive to the exact method we use to estimate aftershock densities.We treat relationships as near-zero if the absolute magnitude of correlation is 0-0.2, weak if the magnitude of correlation is 0.2-0.4,moderate if the magnitude of correlation is 0.4-0.8, and strong if greater than 0.8.
Analyses involving aftershocks presented and discussed in the main text of this study are conducted above the conservative, common M c values (estimated in Section 3).However, in Figures S7-S9 of Supporting Information S1, we provide results based on using all available aftershocks in the respective catalogs (i.e., including aftershocks below the M c ).The results of these are similar to those presented in the main text, but have less confidence, given that aftershocks at incomplete magnitudes are included, and the argument could be made that distributions are systematically missing certain aftershocks by magnitude.

Coseismic Slip and Afterslip
Despite expectations that coseismic slip and cumulative afterslip should be spatially anti-correlated, relationships vary considerably between case studies (Figure 4a).Correlations are moderately positive for the 2004 M w 6.0 Parkfield earthquake, weakly positive for the 2019 M w 7.1 Ridgecrest earthquake, near zero for the 2014 M w 6.0 South Napa, M w 7.6 Izmit earthquake, and M w 6.3 L'Aquila earthquakes, and moderately negative for the 2010 M w 7.2 El Mayor Cucapah and 2011 M w 7.1 Van earthquakes.Most relationships are temporally quasi-stable, but in the Izmit case the relationship becomes increasingly negative, reflecting the strong migration of afterslip along-strike and down-dip away from coseismic slip (Bürgmann et al., 2002).This null result implies no universal relationship between distributions of coseismic slip and afterslip.
Figure 4b shows the relationships between distributions of coseismic slip gradient and cumulative afterslip.These relationships appear broadly similar to those in Figure 4a, but are weaker in magnitude.This implies that coseismic slip gradient is not a suitable indicator of afterslip distribution at these coarse resolutions (e.g., using 2 km cells).Assuming that afterslip is triggered by the stress concentrations at the edges of coseismic slip (Bürgmann, 2018;Perfettini et al., 2018) and that coseismic slip gradient is a suitable proxy for these stress concentrations, it is interesting that coseismic slip gradient does not correlate with afterslip distribution.We discuss the implications of these findings in Section 6.

Coseismic Slip and Aftershock Density
Despite expectations that coseismic slip and cumulative aftershock density should also be spatially anti-correlated, relationships are weak-to-moderately positive in five out of seven case studies (Figure 5a).These tests therefore suggest that coseismic slip may be a potentially useful indicator of where aftershock distributions will occur at these coarse (∼kilometer-scale) resolutions.Figure 5a also indicates that relationships between coseismic slip and cumulative aftershock density are typically more stable through time than relationships between coseismic slip and afterslip (including for the near-zero Izmit relationship, and the weakly negative South Napa relationship).This suggests that in these case studies, aftershocks may have been less migratory than afterslip.Similarly to in the previous section, the relationship between coseismic slip gradient and cumulative aftershock density (Figure 5b) does not appear to imply that coseismic slip gradient is a useful indicator of aftershock distribution.These findings are discussed Section 6.

Afterslip and Aftershock Density
Despite expectations that afterslip and cumulative aftershock density should be spatially correlated, relationships vary considerably between case studies (Figure 6a).Whilst relationships are positive for the Parkfield, El Mayor Cucapah, and L'Aquila earthquakes, they are near zero for the Ridgecrest, Izmit, and Van earthquakes, and moderately negative for the South Napa earthquake.Again, the Izmit case is the least temporally stable relationship becoming increasingly negative through time.This behavior reflects how afterslip and aftershock distributions became increasingly separated with time: afterslip migrated to depth as well as along-strike (Bürgmann et al., 2002) whilst aftershocks appear to migrate along-strike broadly constrained to the upper 15 km of crust.This null result implies no universal relationship between distributions of afterslip and aftershocks, despite the expectations outlined in Section 2. The relationships shown in Figure 6b also suggests that afterslip gradient is also not a useful tool for understanding the spatial distribution if aftershocks at these kilometer-scale resolutions of study.We discuss the implications of these findings in Section 6.

Total Slip and Aftershock Density
The spatial relationship between total cumulative slip (the sum of coseismic slip and afterslip) and cumulative aftershock density is weak-to-moderately positive for six out of seven of our case studies.Figure 7a indicates that only in the case of the South Napa earthquake is aftershock density anti-correlated with the distribution of total modeled slip.We discuss why this may be the case for the South Napa case study in Section 6. Results from this test suggest that total modeled slip is a potentially useful indicator of aftershock density distribution, at these kilometer-scale resolutions of study.Figure 7b indicates that spatial relationships between total slip gradient and aftershock density are broadly similar to those between total slip and aftershock density.However, in the case of the Parkfield and South Napa earthquakes, correlation coefficients are near zero.

The Effect of Test Domain on Results
In Figure 8, we present results from our additional set of tests.As outlined in Section 4.3, these additional tests explore relationships in cells which contain at least 10% of the maximum value of the first test variable at each given time step, which is shown schematically in Figure 8a.Most of the results presented in Figures 8b-8e are broadly similar to those found in the main tests (Sections 5.1-5.4),however, there are some noteworthy differences.For example, following the L'Aquila earthquake, the relationship between coseismic slip and aftershock density is weak-to-moderately positive when considering where both phenomena occurred (Figure 5b, the main tests), which implies that aftershocks broadly occurred where there was coseismic slip.However, this relationship is negative when only considering regions in which significant coseismic slip occurred (Figure 8c, the additional tests), implying that more aftershocks occur where coseismic slip is low.Additionally, the relationship between afterslip and aftershock density appears generally positive across case studies, when only considering regions which underwent significant afterslip (Figure 8d).This includes the South Napa earthquake, for which the relationship was weakly negative in the main tests (Figure 7a).However, whilst interesting, these results offer no constraints on the spatial distribution of the second variable outside the spatial distribution of the first variable.These findings also highlight how seemingly-contradicting claims can be made regarding the spatial relationship between two distributions based on differences in the chosen tested domain.All results are summarized in Table 2.

Discussion
Our tests imply that at the spatio-temporal resolutions considered here, on-fault aftershock density is generally spatially correlated with total modeled slip (in six out of seven cases), which may be pragmatically useful for estimating where aftershocks may broadly occur.However, our overall findings are not consistent with our initial hypotheses, most notably: there is no universal spatial relationship between distributions of coseismic slip and afterslip despite an expectation of universal anti-correlation, spatial relationships between coseismic slip and aftershock density are generally (in five out of seven cases) positive despite expectation of universal anti-correlation, and there is no universal relationship between distributions of afterslip and aftershock density despite an expectation of universal positive correlation.In this section, we therefore discuss why in some case studies: (a) coseismic slip and afterslip can appear to be spatially correlated, although the former implies a velocity-weakening host patch and the latter implies a velocity-strengthening host patch (Bürgmann, 2018;Dieterich, 1979;Perfettini & Avouac, 2007), (b) coseismic slip and aftershocks can appear spatially correlated, although the former implies that the fault patch has already slipped, and (c) afterslip and aftershocks can appear anti-correlated, despite models proposing that they co-migrate (e.g., Kato, 2007;Perfettini et al., 2018).
Our results could be seen as evidence that additional mechanisms such as pore pressure changes (e.g., Jónsson et al., 2003;Piombo et al., 2005), dynamic triggering (e.g., Brodsky & van der Elst, 2014;Pollitz & Sacks, 2002) or secondary triggering (e.g., Cattania et al., 2015;Meier et al., 2014) play an important role in (at least on-fault) aftershock triggering.Whilst the theoretical considerations outlined in Section 2 may still be sound, additional processes may cause the spatio-temporal distributions of aftershock sequences to deviate from expectations that are based solely on coseismic slip and afterslip as the assumed drivers.If this is the case, then current aftershock models based on coseismic slip and/or afterslip, such as Coulomb rate and state (Dieterich, 1994) and afterslip expansion models (e.g., Helmstetter & Shaw, 2009;Kato, 2007;Perfettini & Avouac, 2007;Perfettini et al., 2018) may not be suitable for all case studies and settings.
Our results may also, to some extent, be explained by temporally evolving frictional properties and regions of conditional slip stability in the fault zone.Regions of conditional frictional stability may be able to accommodate either seismic and aseismic slip depending on applied slip velocity (Boatwright & Cocco, 1996;Bürgmann, 2018;Bürgmann et al., 2002;Helmstetter & Shaw, 2009;Scholz, 1998).This could theoretically explain the positive spatial correlation between coseismic slip and afterslip observed following some earthquakes.Frictional properties may also change with time; for example, Veedu and Barbot (2016) suggested that the frictional slip stability of specific asperities in the Parkfield fault zone may be altered by temporal changes in pore-fluid pressure.However, it is still unclear why afterslip would occur on fault patches which have already slipped recently (i.e., coseismically).
Another possible explanation is that the (kilometers-scale) resolution of the slip models used here is too coarse to image fine scale complexity in slip distributions, which are most likely a function of fine scale geometric fault heterogeneity (Milliner et al., 2016) and/or material, rheological, and frictional heterogeneity (Fagereng & Beall, 2021) on the fault.This could theoretically explain why afterslip sometimes appears co-located with coseismic slip and why aftershocks can appear co-located with both coseismic slip and afterslip distributions.Mai and Beroza (2002) proposed that (coseismic) slip distributions are fractal in nature, meaning that asperities likely also exist at all scales.In exhumed faults and shear zones, fractal distributions of compositional asperities are also reported (Fagereng, 2011;Grigull et al., 2012;Kirkpatrick et al., 2021).There is evidence that earthquakes can leave behind unruptured velocity-weakening patches that can rupture later as aftershocks (Das & Aki, 1977), and are models that predict that aftershocks occur on isolated, velocity-weakening asperities surrounded by velocity-strengthening material undergoing afterslip (Kato, 2007;Perfettini et al., 2018), perhaps even as repeating events if loaded to failure more than once (H.Huang et al., 2017;Perfettini et al., 2018;Uchida & Bürgmann, 2019).Considering therefore that coseismic slip, afterslip, and aftershock distributions may occur close to one another on fine spatial scales, it is highly likely that these phenomena may appear co-located in relatively coarse analyses such as this one.This is even more convincing when considering the spatial uncertainties, which are likely on the scale of hundreds of meters to a few kilometers.These reasons also likely explain why our analysis of relationships involving slip gradients do not produce strong, universal correlations, because our calculation of gradients does not reflect fine-scale features, only large features such as the edges of major slip patches.This phenomenon can be highlighted by artificially reducing the resolution of models (by increasing the grid size) in our analyses, which typically results in more positive relationships (results are shown in Figure S11 of Supporting Information S1).
Uncertainties in afterslip and coseismic slip models may have affected our results.As indicated in Section 3, slip modeling has inherent difficulties concerning uncertainty, non-uniqueness, and resolution (Lohman & Simons, 2005;Menke, 2018;Scales & Tenorio, 2001).In this study, efforts are made to analyze high-quality slip models, but although we chose well-studied earthquakes, limitations in the published models (particularly for afterslip) limit the scope of our analyses.Model resolution is a key issue in slip modeling, and whilst the distribution of slip at low (or even moderate resolutions) may be relatively well constrained, fine-scale detail is limited.Abercrombie et al. (2020), for example, highlights difficulties ascertaining whether the coseismic slip of the 2004 M w Parkfield earthquake ruptured through small regions of repeating earthquakes on the fault.To highlight the effect that the choice of model may have on our results (given their variability and uncertainty), we re-conduct analyses involving coseismic slip models for the 2004 M w 6.0 Parkfield earthquake.We reanalyze Figure 8. Spatial relationship results based on using a different test domain to the main tests.Panel (a) shows a schematic example distribution of two variables on the fault plane, and highlights the domains over which we would test our main analysis (green) and our second analysis (purple), which includes only regions with a significant quantity of variable 1. Panels then show the results of spatial relationship analysis using the second test domain: (b) the spatial relationships between coseismic slip and cumulative afterslip in regions where significant coseismic slip occurred, (c) the spatial relationships between coseismic slip and cumulative aftershock density in regions where significant coseismic slip occurred, (d) the spatial relationships between cumulative afterslip and cumulative aftershock density in regions where significant afterslip slip occurred, and (e) the spatial relationships between cumulative total slip and cumulative aftershock density in regions where significant total slip occurred.As there is only one afterslip time-step for the El Mayor Cucapah earthquake, these relationships are shown as a circle.2004), all hosted on the SRCMOD database (Mai & Thingbaijam, 2014).Analyses incorporating these alternative coseismic slip models produce generally similar results to those produced in the main analysis (results are presented in Figure S10 of Supporting Information S1), with some variation reflecting uncertainties in the models themselves.Overall, this gives confidence in our analyses, at the relatively coarse resolutions we are investigating.
More broadly, regarding the confidence in and value of our analyses, we wish to highlight that: (a) whilst aftershock catalogs are associated with location uncertainty (Section 3), our relatively coarse resolutions of analysis will mitigate the impact of this uncertainty on analyses, (b) our methodology makes a number of simplifying assumptions that may conceal nuance in relationships between coseismic slip, afterslip, and on-fault aftershocks (e.g., use of a single best-fitting plane, selecting aftershocks within some distance of this, and projecting aftershocks and slip models orthogonally onto this plane), (c) our study is limited to analyzing on-fault (or in-fault-zone) aftershocks and our results cannot be extrapolated to reflect aftershock behaviors off-fault (outside the fault-zone), (d) it is beyond the scope of this study to analyze additional postseismic mechanisms, such as pore-fluid processes, which may play a role in influencing aftershock distributions (e.g., Ross et al., 2017).This analysis specifically seeks to understand relationships between distributions of coseismic slip, afterslip, and on-fault aftershocks, given hypotheses that both coseismic slip and afterslip are key drivers of aftershock distributions.In the absence of identifying strong, universal relationships, we propose that unresolved fine-scale detail in slip distributions is likely, and that mechanisms such as pore pressure changes (e.g., Jónsson et al., 2003;Piombo et al., 2005) or dynamic triggering (e.g., Brodsky & van der Elst, 2014;Pollitz & Sacks, 2002) likely also play an important role in (at least on-fault) aftershock triggering.Technological and methodological innovation and refinement may provide the means for future studies to analyze slip models and aftershock catalogs of greater resolution and confidence.Potentially significant advances in this space include newer satellites for geodesy, such as Sentinel-1 (Geudtner et al., 2014), and new high resolution seismic catalogs produced through new template-matching or machine learning techniques (e.g., Ross et al., 2018).

Notes on Individual Case Studies
Regarding individual case studies, there are several points worth noting.The 2011 M w 7.1 Van earthquake appears to be the only case-study in which there is an obvious, apparent separation of velocity-strengthening (an assumption due to the spatial distribution of afterslip) and velocity-weakening (an assumption due to the spatial distribution of coseismic slip and aftershocks) parts of the fault (this is clearly highlighted in Figure S5 of Supporting Information S1).In this case, afterslip occurred mostly updip of the coseismic slip and aftershocks (Trasatti et al., 2016), reducing the shallow slip deficit (Dogan et al., 2014).This case study is perhaps most in keeping with simple frictionally-stratified models of the crust, with an upper layer of velocity-strengthening material overlying a velocity-weakening seismogenic zone (e.g., Hillers et al., 2006;Marone et al., 1991;Scholz, 1998), although significant afterslip is not modeled beneath the seismogenic zone.
Observations following the 2011 M w 7.1 Van earthquake directly contrast the 2004 M w 6.0 Parkfield earthquake.In the Parkfield case, coseismic slip, afterslip, and aftershocks occur seemingly in the same place (highlighted in Figure 3).This earthquake occurred adjacent to the creeping section of the San Andreas fault, which appears to contain a high proportion of weak but velocity-strengthening phyllosilicates (Ikari et al., 2011;Moore & Rymer, 2007;Niemeijer, 2018).Here, even small amounts of these materials and/or temporal evolution of pore pressure may significantly alter the frictional slip stability of fault patches (French & Zhu, 2017;Segall & Rice, 1995;Veedu & Barbot, 2016;C. Wang et al., 2017), potentially facilitating afterslip closely intertwined with coseismic slip.There is also evidence of small and isolated velocity-weakening patches in regions understood to otherwise be velocity-strengthening on this section of the fault (Nadeau & Johnson, 1998), implying that afterslip may have surrounded and loaded some patches to rupture (i.e., possible "repeaters" or repeating aftershocks).It is therefore likely that following the Parkfield earthquake, coseismic slip and afterslip distributions were heterogeneous and potentially anti-correlated at finer scales than could be resolved.As an interesting point of note, some previous studies have used aftershock locations to better constrain the coseismic slip distribution itself (e.g., Ziv, 2012).Whilst almost all aftershocks of the Parkfield earthquake appear very close to the main fault plane (as highlighted in Figures 1 and 2a), there is evidence of multiple fault strands separated by up to several hundred meters (Thurber et al., 2004;Zoback et al., 2011), which also could explain why some aftershocks appear 10.1029/2023JB027168 20 of 26 co-located with coseismic slip and afterslip.However, identifying this phenomenon in real-time, at typical study resolutions such as those presented here, would be very challenging.
The 2014 M w 6.0 South Napa earthquake is an outlier throughout our analysis, as aftershocks appear to occur spatially separated from coseismic slip and afterslip (highlighted in Figure S1 of Supporting Information S1).
In addition to this, relative aftershock productivity is low, which some studies have attributed to the spatially confined afterslip distribution and the lack of triggered velocity-weakening patches (Hardebeck & Shelly, 2016;Premus et al., 2022).Considering also that the M w 6.0 South Napa earthquake was relatively high in slip over a much smaller area than the Parkfield earthquake (Floyd et al., 2016), it seems likely that confined and relatively high slip (both coseismic or afterslip) resulted in few aftershocks being triggered close by.Wetzler et al. (2016) also noted a link between high stress drop (megathrust) earthquakes and lower aftershock productivities, proposing that a relatively small rupture area (for a given magnitude) results in a reduced aftershock triggering by coseismic stress changes.

Conclusions
Coseismic slip and afterslip are commonly proposed as drivers of aftershock sequences, but synthesis studies have yet to show systematic evidence of universal driving relationships.In this study, we investigate how the density of on-fault aftershocks correlates with distributions of coseismic slip, afterslip, their sum, and their respective gradients, following seven M w 6.0-7.6 well-studied continental earthquakes.We compile high quality coseismic slip models, afterslip models, and regional seismic data, and project these at short distances onto a best-fitting plane, so that we can undertake spatial correlation analysis.We find that: 1.The spatial relationships between the distributions of afterslip and coseismic slip, and between the distributions of afterslip and aftershock density differ between earthquakes, and between different studies for individual earthquakes.2. Aftershock density is weak-to-moderately correlated with coseismic slip in five out of seven cases, and with total cumulative slip (the sum of coseismic slip and afterslip) in six out of seven cases.Therefore, we propose that the broad distribution of on-fault aftershock density can be approximated by the distribution of total cumulative slip at coarse (kilometers-scale) spatial resolutions.3. Our results imply that the consideration of additional triggering mechanisms such as pore pressure changes and dynamic triggering may be vital to understand aftershock distributions in some settings.4. Additionally, conditional slip stability and/or fine-scale geometric/material/frictional heterogeneity (resulting in unresolved coseismic slip and afterslip heterogeneity at the scale of slip models) within the fault zone likely explains why distributions of coseismic slip, afterslip, and aftershocks may appear co-located.As a likely result, we also find that the roles that slip gradients (treated as proxies for new stress concentrations) play in driving aftershock distributions are not clear.5. Different choices of the spatial domain over which relationships are tested can result in seemingly-contradictory results (e.g., aftershocks may broadly occur located with coseismic slip, but aftershock density and coseismic slip may be anti-correlated when only looking at regions that undergo significant coseismic slip).
Future works may attempt to support the hypotheses outlined in Section 2, which we argue are still theoretically sound.However, higher slip model resolutions and more accurate aftershock locations may be required, which are significant methodological challenges.

Data Availability Statement
The codes and data used in this study are provided through the repository: Churchill (2023).Several of the slip models used are also available online through: Jiang et al. (2021b), Floyd et al. (2016), Yue et al. (2021), http:// seismo.berkeley.edu/∼burgmann/RESEARCH/TURKEY/turkey.htmlTrasatti et al. (2023), and the SRCMOD database (Mai & Thingbaijam, 2014).Additionally, Chris Rollins (El Mayor Cucapah) and Daniele Cheloni and Nicola D'Agostino (L'Aquila) allowed us to put their afterslip models (which were not stored online) in our repository.The Ridgecrest afterslip model by Kang Wang (kangwang@berkeley.edu) is not publicly available (as of October 2023) but will be published in the near future.The seismic data utilized are available through the Northern California Earthquake Catalog (NCEC) (Waldhauser, 2009;Waldhauser & Schaff, 2008), the California We thank Kelian Dascher-Cousineau, an anonymous second reviewer, Rachel Abercrombie, and the Associate Editor for reviews which also improved this manuscript.

Figure 1 .
Figure 1.Maps showing seismicity in the 18 months following each mainshock and a given fault geometry model for each of the seven earthquakes studied: (a) the 2004 M w 6.0 Parkfield, (b) 2014 M w 6.0 South Napa, (c) 2019 M w 7.1 Ridgecrest, (d) 2010 M w 7.2 El Mayor Cucapah, (e) 1999 M w 7.6 Izmit, (f) 2011 M w 7.1 Van, and (g) 2009 M w 6.3 L'Aquila earthquakes.Panels show the (above M c ) seismicity in the 18 months following each mainshock (catalogs are given inTable 1), with color indicating time since the mainshock and size indicating magnitude.The afterslip model fault geometry (see Table1) is given in black in each case, except for in panel (c), where we show the foreshock fault plane as well, and panel (d), where we include both the (spatially-separated) coseimsic slip model too, for context.

Figure 2 .
Figure2.The distance-distribution between aftershocks and the best-fitting plane for each mainshock.These histograms include all events up to 18 months after the mainshock, above the common M c used for analysis: 1.5 for the Parkfield and South Napa earthquakes (b), 2.2 for the Ridgecrest and El Mayor Cucapah earthquakes (b), and 2.7 for the Izmit, Van, and L'Aquila earthquakes (c).As the best-fitting plane is extended through the entire initial selection area, these selections may include some along-strike events away from coseismic rupture and afterslip.In each case, the projection distance cut-off is given by a vertical line.

Figure 3 .
Figure 3.An illustrative example of the regridding of coseismic slip and afterslip models, the calculation of slip gradients, and a Gaussian kernel density estimate of aftershock locations following the 2004 M w 6.0 Parkfield earthquake: (a) shows the coseismic slip model by Jiang et al. (2021a) projected and regridded (using a 2-by-2 kilometer grid) onto the best-fitting plane, and (b) shows the normalized coseismic slip gradient magnitude field calculated from panel (a).Panels (c and d) show the same as a and b, but for afterslip at approximately 1 month (for illustrative purposes) (Jiang et al., 2021a).Panel (e) shows the projected location of aftershocks within 2.5 km of the best-fitting plane at 1 month (for illustrative purposes) and (f) shows a Gaussian kernel density estimate of e, on the same common grid as the coseismic slip and afterslip models.Distances across the x-axes are relative to an arbitrary along-strike reference point.This figure is provided for the other six earthquakes studied in Figures S1-S6 of Supporting Information S1.

Figure 4 .
Figure 4. Spatial relationships between the distributions of coseismic slip and cumulative afterslip through time: (a) shows the Spearman's rank correlation coefficient between the distributions of coseismic slip and cumulative afterslip and (b) shows the Spearman's rank correlation coefficient between the distributions of coseismic slip gradient and cumulative afterslip.As there is only one afterslip time-step for the El Mayor Cucapah earthquake, this relationship is shown as a circle.The expected relationship is shown as a gray box in panel (a).

Figure 5 .
Figure 5. Spatial relationships between the distributions of coseismic slip and cumulative aftershock density through time: (a) shows the Spearman's rank correlation coefficient between the distributions of coseismic slip and cumulative aftershock density and (b) shows the Spearman's rank correlation coefficient between the distributions of coseismic slip gradient and cumulative aftershock density.The expected relationship is shown as a gray box in panel (a).

Figure 6 .
Figure 6.Spatial relationships between the distributions of cumulative afterslip and cumulative aftershock density through time: (a) shows the Spearman's rank correlation coefficient between the distributions of cumulative afterslip and cumulative aftershock density and (b) shows the Spearman's rank correlation coefficient between the distributions of cumulative afterslip gradient and cumulative aftershock density.As there is only one afterslip time-step for the El Mayor Cucapah earthquake, this relationship is shown as a circle.The expected relationship is shown as a gray box in panel (a).

Figure 7 .
Figure7.Spatial relationships between the distributions of cumulative slip (the sum of coseismic slip and afterslip) and cumulative aftershock density through time: (a) shows the Spearman's rank correlation coefficient between the distributions of cumulative slip and cumulative aftershock density and (b) shows the Spearman's rank correlation coefficient between the distributions of cumulative slip gradient and cumulative aftershock density.As there is only one afterslip time-step for the El Mayor Cucapah earthquake, this relationship is shown as a circle.
Integrated catalog by the Southern California Seismic Network (SCSN) (SCEDC, 2013), and the International Seismological Centre (ISC) Bulletin (ISC, 2022).R. M. Churchill and M. J. Werner were supported by the European Union H2020 program (Grant 821115, Realtime earthquake rIsk reduction for a reSilient Europe [RISE]).R. M. Churchill acknowledges the support of a NERC GW4+ Doctoral Training Partnership studentship [NE/L002434/1].J. Biggs and M. J. Werner were also funded by the NERC Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET, http://comet.nerc.ac.uk), a partnership between UK Universities and the British Geological Survey.J. Biggs is also funded by a Leverhulme Prize (PLP-2018-362).A. Fagereng received funding from the European Union H2020 programme, ERC Starting Grant MICA [No 715836].We thank the following individuals who provided, or helped find, coseismic and afterslip models: Roland Bürgmann, Daniele Cheloni, Nicola D'Agostino, Semih Ergintav, Wangpeng Feng, Elizabeth Hearn, Junle Jiang, Fred Pollitz, Chris Rollins, Elisa Trasatti, Kang Wang, Sam Wimpenny, and Han Yue.We thank Robert Myhill and Jessica Hawthorne for discussion that improved this work.
Note.Summary of results across the main and additional sets of tests.Columns designate the relationship being tested, the type of test, and the approximate results for the 2004 M w 6.0 Parkfield, 2014 m w 6.0 South Napa, 2019 M w 7.1 Ridgecrest, 2010 M w 7.2 El Mayor Cucapah, 1999 M w 7.6 Izmit, 2011 M w 7.1 Van earthquake, and 2009 M w 6.3 L'Aquila earthquakes, respectively.P, N, and 0 correspond to broadly positive, negative, and (near) zero relationships, respectively, and w and m stand for weak and moderate, respectively.As results are time dependent, these are approximations.

Table 2 A
Dreger et al. (2005)5)oss the Two Sets of Testscorrelations between distributions of coseismic slip and afterslip, coseismic slip and aftershock density, and total cumulative slip and aftershock density, using alternative coseismic slip models byCustódio et al. (2005),Dreger et al. (2005), and Ji (