A New Focal Mechanism Calculation Algorithm (REFOC) Using Inter‐Event Relative Radiation Patterns: Application to the Earthquakes in the Parkfield Area

Accurate earthquake focal mechanisms are essential for solving fault zone structure, estimating stress variations, and assessing seismic hazards. Small earthquakes' focal mechanisms are usually solved using P‐wave first‐motion polarities and/or S‐/P‐wave amplitude ratios, which are limited due to the low signal‐to‐noise ratio of small‐earthquake waveforms and a limited number of three‐component seismograms. To increase the number of high‐quality focal mechanisms, we develop a method that utilizes the inter‐event relative radiation patterns to perform a joint focal mechanism inversion of numerous clustered events (called REFOC). The method first uses P‐wave polarities and S‐/P‐wave amplitude ratios to constrain the initial solutions and then combines these solutions and the inter‐event P‐/P‐wave and S‐/S‐wave amplitude ratios to refine solutions. For example, we apply the method to 38,413 earthquakes in the Parkfield region, California. The REFOC outperforms traditional methods with 57% more solutions with <55° uncertainties (>70% of the catalog events) and 126% more solutions with <25° uncertainties (>40% of the catalog events), illuminating unprecedented fine‐scale rupture processes. Instead of rupturing along the main fault, many M < 2 focal mechanisms show >45° angular differences from the 2004 Mw 6.0 Parkfield earthquake. The variation of focal mechanism properties is spatially related to the variation of fault strength and geometry, and temporally correlated with the shear stress variations before, during, and after the 2004 Mw 6.0 earthquake. The observations highlight the potential of applying REFOC to monitor unprecedented details of fault zone structure and stress field, providing new insights into fault rupture physics, seismotectonic processes, and seismic hazards.

Due to the small radius and short rupture duration of small earthquakes, they are usually approximated as a double-couple point source, which has two nodal planes and divides the reference focal sphere into four quadrants with different first-motion polarities (Lay & Wallace, 1995). Note that many small earthquakes may also have considerable non-double-couple components but are difficult to constrain, like the co-rupture of multiple faults, fluid-related earthquakes, and rock-damage-related earthquakes (e.g., Cheng et al., 2021;Martínez-Garzón et al., 2017;Miller et al., 1998). Traditionally, the polarities are picked from seismic stations, mapped back to the focal sphere based on ray paths, and used to fit four quadrants and two nodal planes, such as the FPFIT method (Reasenberg, 1985). These types of methods have considerable uncertainties (Hardebeck & Shearer, 2002) because of various possible errors in polarity picking, velocity model, and event location. To further constrain the nodal planes, many studies increase observations per event using other waveform features, like S/P amplitude ratios (e.g., Hardebeck & Shearer, 2003;Julian & Foulger, 1996;Kisslinger et al., 1981), and absolute P-wave and S-wave amplitudes (e.g., Ebel & Bonjer, 1990;Kwiatek et al., 2016;Nakamura et al., 1999;Rögnvaldsson & Slunga, 1993) (Figure 1). Recently, the rapid development in machine learning algorithms also helped increase the observations by using high-accuracy polarity picker and phase picker trained using millions of labeled waveforms (Cheng & Ben-Zion, 2020;Ross et al., 2018).
Despite the above-mentioned progress, the number of observations for each earthquake is still limited because of sparsely distributed stations, the difficulties in polarity picking, the requirement of three-component waveforms for S-wave and P-wave amplitude measurements, as well as the challenges in suppressing the other factors affecting the recorded amplitudes, like earthquake magnitude, path attenuations, station site effects, and noises from the other sources. The limited station coverage and the small number of high-quality seismic records cannot provide sufficient constraints on focal mechanism estimations. Therefore, even for the areas with dense station Supervision: Richard M. Allen, Taka'aki Taira Visualization: Yifang Cheng Writing -original draft: Yifang Cheng, Richard M. Allen, Taka'aki Taira Writing -review & editing: Yifang Cheng, Richard M. Allen, Taka'aki Taira Figure 1. Input data used in different focal mechanism calculation methods. All methods use P-wave first-motion polarity as input. HASH and REFOC also incorporate S-/P-wave amplitude ratio into the calculation. The REFOC method performs first iteration focal mechanism calculation following the HASH method and further refines the focal mechanisms by performing a second iteration calculation using the relative amplitude ratios between the target event (tar) and the focal mechanisms of neighboring events (nei) obtained from the first iteration calculation.

10.1029/2022JB025006
3 of 19 coverage like southern California, only around 20% of the located events have uncertainties less than 25° (Yang et al., 2012). The limited number of high-quality focal mechanisms further leads to large uncertainties in the derived fault orientations, slip directions, and stress fields, which impedes the resolution and interpretability of the obtained results. Therefore, new methods for improving the focal mechanism resolutions and further enhancing focal mechanism quality and quantity are still in great demand.
Although each individual small earthquake has limited high-quality observations, small earthquakes occur more frequently than large earthquakes and tend to be clustered in space, like earthquake swarms and mainshock-aftershock sequences (Vidale & Shearer, 2006;Zaliapin & Ben-Zion, 2013). The highly clustered earthquakes provide a great opportunity to utilize the relative relationships between co-located earthquakes (inter-event relationships) to suppress the common path and site effects and extract reliable relative source properties. This strategy has been widely applied to constrain earthquake properties, like earthquake location, composite focal mechanism, moment tensor, and stress drop, which has made dramatic improvements in reducing solution uncertainties (e.g., Jia et al., 2022;Plourde & Bostock, 2019;Shearer, 1997;Shearer et al., 2006;Shelly et al., 2016;Waldhauser & Ellsworth, 2002) and been applied to large seismic datasets (Lin et al., 2007;Waldhauser & Schaff, 2008). Despite of these successful applications, the inter-event relationships have not been used in the focal mechanism calculations because many previous studies assumed that co-located earthquakes are located on the similar fault structure and share similar radiation patterns (Godano et al., 2014;Shelly et al., 2016). In addition to the assumption of physical processes, the implementation of inter-event relationships of hundreds of neighboring events to calculate focal mechanisms of millions of small earthquakes requires high computational cost. To overcome these limitations, we incorporate the relative radiation patterns into the traditional focal mechanism calculation method (HASH; Hardebeck & Shearer, 2002, 2003, propose a new focal mechanism calculation method (REFOC: RElative FOCal mechanism calculation), and use one NVIDIA V100 graphics processing unit to perform the method. We first present the details of methodology (Section 2) and then apply the method to earthquakes in Parkfield area (Section 3). The results are discussed in more detail in Section 4.

Methods
In this section, we introduce the algorithm and computational details underlying the REFOC method. We first explain the required input data for the focal mechanism calculation. Next, we describe the basic algorithm used by the REFOC method, which first calculates the initial focal mechanisms using polarities and S-/P-wave amplitude ratios (first iteration) and computes final focal mechanisms by adding the initial focal mechanisms as well as the inter-event relative P-/P-wave amplitude ratios and S-/S-wave amplitude ratios into the calculation (second iteration). Finally, we show the approach that REFOC uses to assess the uncertainties of focal mechanism solutions. The workflow is outlined in Figure 2.

Data Preparation
The fundamental input data for the REFOC algorithm are P-wave first-motion polarities and P-wave and S-wave amplitudes obtained from earthquake waveforms. Earthquake polarities are usually picked by data analysts. Many automated algorithms are proposed to pick polarities based on the first local maximum after P-arrivals (Chen & Holland, 2016), machine learning models trained using labeled waveforms (Ross et al., 2018), or relative polarities obtained through inter-event waveform cross-correlations (Shelly et al., 2016). For the results presented in this study, we only use the catalog polarities picked by data analysts. In addition to polarities, we also need to measure P-wave and S-wave amplitudes. We apply a 1-10 Hz band-pass-filter to the raw waveform data and choose 0.5 s before to 1.5 s after P-wave and S-wave arrivals as the signal windows and 2.5-0.5 s before P-wave arrivals as the noise windows. We take the difference between the maximum and minimum amplitude values in each time window to be the estimated signal and noise amplitudes. If the time difference between P-wave and S-wave arrivals is larger than 2 s and the SNR is larger than 3, we incorporate the P-wave and S-wave amplitudes into focal mechanism calculation. Note that the S-/P-wave amplitude ratios are only calculated using the vector summation of amplitudes over three-component waveforms.
For focal mechanism calculation, the measured polarities, P-wave amplitudes, and S-wave amplitudes need to be projected to the focal sphere around the source. To do so, we first compute the take-off angles and azimuths for all event-station pairs. Unlike azimuths, take-off angles are highly sensitive to the used velocity models and event depths, which may have errors and affect the best-fitting focal mechanism solutions. Based on the analysis in Hardebeck and Shearer (2002), different velocity models and depths can result in different focal mechanisms with angular differences as large as 45°. To consider these errors and avoid over-fitting solutions to the errors and other noise in the data, we collect solutions from calculations using different take-off angles computed from various combinations of possible source depths and 1D velocity models. The source depth is randomly perturbed within −1 to 1 km around the relocated depth based on the depth uncertainties of events and the velocity model is chosen from a set of models spanning a reasonable range of velocity structures.

Amplitude-Magnitude Scaling
Since we plan to calculate focal mechanisms using inter-event relative radiation patterns, which are also affected by the inter-event magnitude differences, we need to first determine the effect of inter-event magnitude differences on the amplitude ratios: where is a constant that describes the amplitude-magnitude scaling. Since local magnitude is defined as being proportional to the logarithm of the peak-to-peak S-wave amplitude on a Wood-Anderson seismograph, = 1 (Gutenberg & Richter, 1942;Richter, 1935). However, the systematic analysis of earthquake data in California shows that earthquake magnitude estimation depends strongly on the attenuation effects along the path, site effect, and frequency range so that the amplitude-magnitude scaling may not be exactly equal to 1 (Hutton & Boore, 1987;Uhrhammer et al., 2011). Moreover, for M < 3 earthquakes in Northern California, Northern California Seismic System (NCSS) usually used the magnitude estimated from the coda wave durations instead of the S-wave amplitudes as the preferred magnitude. To better estimate the effect of magnitude differences on the amplitude differences and better constrain the focal mechanisms, we decided to obtain the empirical amplitude-magnitude scaling from the event data. To do so, we calculate the magnitude differences and amplitude ratios between all available event-station-channel combinations with inter-event Figure 2. A flowchart of the new focal mechanism calculation method using NCSS event ID 73636545 as an example. In the first iteration, the polarities (blue dots and crossings) and S-/P-wave amplitude ratios (green squares scaled by S-/P-wave amplitude ratios) are used to obtain all acceptable focal mechanism solutions (gray lines). The second iteration adds inter-event P-/P-wave and S-/S-wave amplitude ratios (triangles colored by the number of amplitude ratios at each station) into the calculation to obtain a better solution (colored lines) with a smaller uncertainty (colored circles). hypocentral distances <1 km and the SNR of amplitude measurements larger than 3. We then use the linear fit to determine the best-fitting slopes for the P-wave amplitude ratios and S-wave amplitude ratios relative to the corresponding differential catalog magnitudes, which we assign to and and use them in the focal mechanism calculation.

The REFOC Algorithm
Following the data preparation step, the focal mechanism estimation begins. This section provides a conceptual outline of the two-step focal mechanism calculation algorithm. The algorithm first uses P-wave first-motion polarities and S-/P-wave amplitude ratios to estimate focal mechanisms. The estimated focal mechanisms are used to calculate theoretical S-wave amplitude ratios and P-wave amplitude ratios between co-located events, which are added to the second round of focal mechanism calculation to further constrain the solutions. The detailed steps are listed as follows.
Step 1. (first iteration) For each combination of event depth and 1D velocity model, we obtain the take-off angles and perform a grid search over strike, dip, and rake (in this study, we use 1 • step for the grid search) to identify solutions that minimize the averaged misfit of polarities 1 : Mechanisms with slightly more polarity misfits are also accepted considering the possible polarity errors due to the quality of waveform records and manual picking errors. In this study, we allow up to 10% misfit polarities or the minimum misfit plus 5% if this is greater.
Step 2. (first iteration) From the set of acceptable mechanisms in step 1, we then extract mechanisms that minimize the averaged L1-norm misfit of the S-/P-wave amplitude ratios: The theoretical S-/P-wave amplitude ratios ( 10 ( )) are calculated by assuming a pure double-couple source (Shearer, 2019). We allow the averaged misfit of 2 up to 0.3 (about one order of amplitude differences) or the minimum misfit plus 0.15 if this is greater.
Step 3. (first iteration) After collecting all acceptable solutions with the averaged misfits 1 and 2 less than the chosen thresholds from all trials with various event depths and 1D velocity models, we extract the preferred solution that maximizes the number of acceptable mechanisms within 45° angular differences from it and minimizes the averaged angular difference (avg_rot) between the preferred solution and the other acceptable mechanisms. The angular difference between two focal mechanisms is measured using the method in Kagan (1991). The avg_rot is used to represent the focal mechanism uncertainty. Step 4. (second iteration) Find all neighboring events with initial focal mechanism solutions in the first round and locate close to the target event. In this study, we use the nearest 500 neighboring events within 1 km radius from the target events. Note that each event in the catalog can be treated as the target event and also the neighboring event of other target events.
Step 5. (second iteration) Obtain the P-wave amplitude ratios between the target event and its neighboring events at all available station-channel combinations.
Step 6. (second iteration) From a set of acceptable mechanisms in step 4, we extract mechanisms that minimize the L1-norm misfit of the logarithm of the inter-event P-wave and S-wave amplitude ratios 3 based on the following equation: Since the target event and its neighboring events share similar path and site effects, we only consider the amplitude differences caused by the differences in magnitude and focal mechanism when calculating synthetic amplitude ratios: where ( 10 ( )) are solved using the magnitude differences according to Equation 1 and ( 10 ( )) are calculated by assuming a pure double-couple source (Shearer, 2019) and the neighboring events' focal mechanisms are the preferred solutions obtained from Step 3. We allow the averaged misfit of 3 up to 0.3 (about 1 order of amplitude differences) or the minimum misfit plus 0.15 if this is greater.
Step 7. (second iteration) Collect all solutions with the misfits 1 and 2 averaged from all stations as well as 3 averaged from all station-channel-phase-neighboring-event combinations less than the chosen thresholds (steps 1, 2, and 6) from step 6, from which obtain the preferred solution and solution uncertainty following the details in step 3.

Workflow and an Example Event
Figure 2 outlines the workflow for the focal mechanism calculation as well as an example event (NCSS event ID: 73636545). There are 16 stations used for calculating the focal mechanism with 16 polarities (Step 1) and 12 S/P amplitude ratios (step 2). Each station has at most 1 polarity (blue dots and crossings) and 1 S/P amplitude ratio (green squares). In contrast, there are 10 neighboring events within 1 km hypocentral distance. Leading to multiple inter-event S-wave and P-wave amplitude ratios from multiple event-pair/channel/phase combinations (colored triangles) at each station (steps 4 and 5), which significantly increases the number of inputs. By sequentially adding these input measurements to constrain the focal mechanism (steps 1, 2, and 6), the number of acceptable solutions decreases gradually (gray nodal planes), leading to a better focal mechanism solution (green and red nodal planes) with much smaller solution uncertainty (green and red circles at the nodal plane intersections) (steps 3 and 7).

Application to the Parkfield Area
The Parkfield section of the San Andreas fault (SAF) connects the northern creeping section and southern locked section, where there were many large historic earthquakes (Murray et al., 2001). The Parkfield area is known for hosting recurring moderate-sized (M ∼ 6) earthquakes (Bakun & McEvilly, 1984) including two M ∼ 6 earthquakes in 1966 and 2004. These two earthquakes mainly ruptured the area between Middle Mountain (MM) and Gold Hill (GH) with initiation locations near MM and GH, respectively (Figure 3; Brown, 1967;Rymer et al., 2006). The mapped fault traces and the modeling result of long-term evolution of fault orientations due to fault strength variations along the fault suggest the twist of fault structures near the MM and GH (Perrin et al., 2019). However, both M > 4 focal mechanisms and the seismicity distribution suggest highly straight and localized fault zone dominated by right-lateral strike-slip ruptures along the main fault ( Figure 3). Neither detailed imaging of complex fault geometry at depth nor its relationship to stress evolution and major earthquake ruptures have been elucidated. The fine-scale fault orientations and slip motions of small earthquakes illuminated from the focal mechanisms using REFOC method provide a great opportunity to investigate the variation of fault zone structure and its impact on seismicity behavior and stress field in the Parkfield area.

Overview of the Input Data
To demonstrate the use and efficacy of the REFOC method, we apply it to all available 38,413 relocated events (Waldhauser & Schaff, 2008, extended to later years) in the box AA′ in the Parkfield area from 1984 to 2021 ( Figure 3). Around 89% of the events are located within 1 km of the main fault trace and the 95% confidence level of depth uncertainties in the catalog is 2 km vertically. Therefore, we randomly disturb event depth from −1 to 1 km from the initial depth during focal mechanism calculation. Earthquakes in this region are well-recorded by the surrounding stations operated by the Northern California Seismic Network, with the 2004 M w 6.0 Parkfield earthquake that occurred in the study period. About 16% of the events in this region have focal mechanism solutions in the original focal mechanism catalog calculated using the FPFIT method (Reasenberg, 1985).
In this study, we re-calculate the focal mechanisms in the Parkfield area using the REFOC algorithm. We utilize the polarities, P-wave and S-wave phases manually picked by data analysts, as well as the relocated earthquake catalog archived by the Northern California Earthquake Data Center (NCEDC). The P-wave velocity models used are shown in Figure S1 in Supporting Information S1 and the S-wave velocity models are derived from P-wave velocity model by assuming that the P-/S-wave velocity ratio equals 1.732. We utilize the earthquake waveforms from NCEDC around the P-wave and S-wave arrivals to obtain S-/P-wave amplitude ratios for each individual event as well as the inter-event P-/P-wave amplitude ratios and S-/S-wave amplitude ratios. Figure 4 shows the overview of input data for focal mechanism calculation in the Parkfield area. Before 2003, there were limited stations deployed in Northern California and the majority of them only have a vertical component. As a result, during this period, there are about 70% of events with more than 8 polarities, about 25% of events with more than 15 polarities, and nearly 0% of events with S/P amplitude ratios. After 2003, many three-component seismometers were deployed around the Parkfield area, leading to over 95% of events with more than 8 polarities, and over 60% of events with at least 16 polarities and 8 S/P amplitude ratios. Because of the large number of event-pair/channel combinations for clustered events, there are many more inter-event P-wave and S-wave amplitude ratios than the polarities and S/P amplitude ratios. There are more than 43% of events with at least 24 inter-event P-wave and S-wave amplitude ratios.
We fit the linear slope between the logarithm of the obtained P-/P-wave and S-/S-wave amplitude ratios of all combinations as well as the corresponding differential catalog magnitudes ( Figure 5). The determined best-fitting slopes between magnitude differences and the logarithm of amplitude ratios are 1.05 ( ) for P-wave and 1.12 ( ) for S-wave.

Results and Comparisons With the HASH Algorithm
Since the REFOC algorithm is built on the HASH algorithm, it is instructive to check the consistency and the differences between the results obtained using REFOC and those obtained using the HASH algorithm (Hardebeck & Shearer, 2003). For the comparison, we used the identical data inputs (Section 3.3) except that the REFOC algorithm requires inter-event amplitude ratios. Figure 6 shows the misfit between the input observations and the corresponding synthetic values using the output focal mechanisms. The focal mechanisms solved using the HASH and REFOC algorithms share similar polarity misfits and S/P amplitude ratio misfits. Most events have the averaged polarity misfit ( 1 ) smaller than 0.3 and the averaged logarithm S/P amplitude ratio misfits ( 2 ) between 0.5 and 1.5 (Figures 6a and 6b). We further check the consistency between the focal mechanisms determined using these two algorithms. For each focal mechanism pair, we measure the similarity of two focal mechanisms using the Kagan angle, which is defined as the minimum rotation angle needed to make two focal mechanisms identical (Kagan, 1991). For all 27,750 events with nodal plane uncertainties less than 55° using the REFOC algorithm, over 80% of events have Kagan angles less than 40°. Therefore, the results obtained using the REFOC are highly consistent with those using the HASH. To evaluate the reliability of the results, we further compared the focal mechanisms obtained using HASH and REFOC algorithms with the NCSS moment tensor solutions obtained by waveform modeling (Dreger, 2018). The focal mechanisms obtained using REFOC algorithm are more consistent with moment tensor solutions compared with those using HASH algorithm (Figure 6d), suggesting that the REFOC method can help to obtain more accurate focal mechanism solutions. Based on the HASH algorithm, the REFOC algorithm incorporates the inter-event amplitude ratios into calculation. Overall, most events have the averaged logarithm inter-event P-/P-wave and S-/S-wave amplitude ratio misfits ( 3 ) between 0.2 and 1.0, which is relatively small compared with the misfits of S-/P-wave amplitude ratios ( 2 ) (Figures 7a and 7b). The incorporation of inter-event amplitude ratios further constrains the focal mechanisms and helps to obtain more focal mechanisms with higher resolution and less uncertainties (Figure 7c). Figures 7b-7d show the temporal histograms of focal mechanisms solved using FPFIT, HASH, and REFOC, respectively. Note that we only show the focal mechanisms with more than eight P-wave first-motion polarities. The FPFIT method only uses the P-wave first-motion polarities, resulting in 6,190 focal mechanisms (16% of all catalog events) determined by NCSS. By incorporating the S/P amplitude ratios, the HASH method can help to solve 17,484 (46% of all catalog events) and 7,478 (19% of all catalog events) focal mechanisms with uncertainties less than 55° and 25°, respectively. The REFOC algorithm can solve 27,750 (72% of all catalog events) and 16,466 (43% of all catalog events) focal mechanisms with uncertainties less than 55° and 25°, respectively. Since the REFOC algorithm can help to solve high-quality focal mechanism solutions of nearly half of the located events, most spatiotemporal analyses that have been applied to seismicity can also be applied to focal mechanisms and help to reveal more detailed earthquake and fault zone properties.
One of the important fault zone properties is fault geometry, which can be inferred using the spatial distributions of earthquakes. In contrast, each individual focal mechanism contains the fault orientation information Figure 6. The 2D histograms of the (a) polarity misfit, (b) S-/P-wave amplitude ratio misfits of the focal mechanisms solved using the HASH and REFOC methods, respectively. (c) Histogram of the Kagan angles between common focal mechanisms solved using the HASH and REFOC methods. (d) Histogram of the Kagan angles between common focal mechanisms solved using (blue) the HASH and moment tensor inversion methods and those using (green) the REFOC and moment tensor inversion methods. Vertical lines denote the 80th percentile of the histogram. of each small earthquake and can reveal the orientations of much smaller-scale faults. Figure 8 shows all focal mechanisms with uncertainties less than 55° solved using the HASH and REFOC algorithms. To compare diverse focal mechanisms in a uniform way, we further express the mechanisms' faulting styles on a continuous scale from −1 to 1, with normal faulting having a value of −1, strike-slip denoted by 0 and thrust faulting denoted by 1 . Most M > 2 earthquakes show strike-slip focal mechanisms due to the simple main fault geometry in the Parkfield area. In contrast, M < 2 earthquakes show more diverse faulting styles with many transpressional and transtensional faulting events, which may occur on the secondary faults and have not been observed from M > 2 earthquakes (Thurber et al., 2006). The increase of focal mechanism diversity with decreasing magnitude is not only observed from the M < 2 earthquakes but also from the well-constrained M > 2 earthquakes, suggesting that this increase of focal mechanism diversity is not purely caused by the increasing uncertainties of the focal mechanism solutions for small earthquakes but might be due to the increasing complexity for small scale faults. Small-magnitude earthquakes have smaller source areas and require less stress perturbations for the final rupture, which is more sensitive to small-scale spatiotemporal variation of stress field and also provide a great opportunity to track the spatiotemporal variations of fault zone properties in the Parkfield area (Figure 8).

Spatiotemporal Variations of Focal Mechanism Properties
To compare the small earthquake focal mechanisms with the main fault rupture, we calculate the angular differences between the focal mechanisms of the M w 6.0 Parkfield earthquake, which is highly similar to the other M ≥ 4 earthquakes in this area, and the others using the Kagan angle (Kagan, 1991). Figure 9 shows the locations of all focal mechanisms with nodal plane uncertainty less than 25° solved using the HASH and the REFOC algorithms. These two algorithms provide similar spatial distribution of focal mechanisms because most events are very localized around the main fault horizontally and at several narrow depth ranges vertically. The study area can be divided into three sections based on the inter-seismic slip rate estimated from surface geodetic data (Murray et al., 2001): the locked section around the M w 6.0 Parkfield earthquake with near-zero slip rate (from 13 km NW to 15 km SE of the mainshock epicenter); the creeping section with over 20 mm/yr slip rate (27-50 km NW to the mainshock epicenter); and the transition section in between (13-27 km NW to the mainshock epicenter) with ranging from 0 to 25 mm/yr. The slip rate in the shallow part (depth < 7 km) of the transition section is lower than the slip rate in the deeper part (depth > 7 km). For the events in the locked section around the M w 6.0 Parkfield earthquake, most of them are located below 4 km depth with focal mechanisms highly consistent with the mainshock (Figures 9i, 9j, 9p, and 9q). Events in the transition section are distributed in a wide depth range from 1 to 15 km with shallow focal mechanisms (depth < 7 km) consistent to the mainshock while the deeper  ones (depth > 7 km) are not (Figures 9h and 9o). In the creeping section, most events are in the top 7 km with focal mechanisms highly different from the M w 6.0 mainshock (Figures 9e, 9f, 9i, and 9m). Most of the M ≥ 4 earthquakes are located in around 10 km depth in the transition and locked sections with focal mechanisms highly similar to the M w 6.0 Parkfield earthquake.
For a better understanding of the physics of earthquakes, it is important to know if there are any changes of earthquake parameters with time and how it is related to the large earthquake ruptures or any other processes. Figure 10 shows the temporal variation of focal mechanisms along AA′. The additional focal mechanism solu-  (a, b). Large and small beachballs denote the M w 6.0 earthquake and the other M ≥ 4 earthquakes, respectively. Black dashed lines denote the occurrence time of the M w 6.0 Parkfield earthquake. Histograms of the angular difference from the M w 6.0 earthquake for the focal mechanisms before (blue), 0-1 year after (green), 1-5 years after (orange), and 5-10 years after (red) the M w 6.0 earthquake in the (c, d) Creeping, (e, f) Transition, (g, h), and Locked sections that are solved using the (c, e, g) HASH and (d, f, h) REFOC methods. tions solved using the REFOC algorithm provide more detailed temporal evolution of earthquake properties along the fault, especially for small magnitude earthquakes in the locked section and after the M w 6.0 Parkfield mainshock ( Figure 10, Table S1 in Supporting Information S1). Therefore, we mainly focus on the analysis of focal mechanisms using REFOC algorithm. To quantify the temporal variations of focal mechanisms in the REFOC catalog, we calculate the normalized histogram of the angular differences from the M w 6.0 Parkfield earthquakes for events occurring before and after the 2004 M6.0 mainshock in the creeping, transition, and locked sections, respectively (Figures 10c-10h). In the creeping section, the focal mechanisms show highly varying angular differences from the M6.0 mainshock, suggesting that the focal mechanisms are highly heterogeneous in the creeping section in all time periods. In the transition section, the focal mechanisms are more similar to the M6.0 mainshock within 0-5 years after the mainshock and gradually become more different from M6.0 mainshock afterward. On the contrary, in the locked section, the angular difference from M6.0 mainshock is very small before the mainshock, increased within 0-5 years after the mainshock. The angular differences become slightly smaller in the period 5-10 years after the mainshock but with smaller deviation (Figure 10h), which might suggest the change of tectonic stress direction. We further investigate the faulting type variations in space and time in the study area ( Figure 11). Most events are strike-slip faulting between the MM and GH while there are many reverse faulting events NW to the MM and normal faulting events SE to the GH, especially after the 2004 mainshock, corresponding to the 5° restraining bend near MM and extensional jog near GH (Figures 9a  and 9b; Lindh & Boore, 1981), respectively. These results suggest that the focal mechanism variations show strong connection to fault geometry, fault strength, and the occurrences of major ruptures. Many of these features cannot be captured using the focal mechanisms of large magnitude events because of the limited quantity and their highly consistent focal mechanisms (Figure 8).

Comparison of REFOC Method With Other Source Inversion Methods
Although other source characterization methods that utilize the advantage of clustered events Peng & Zhao, 2009;Shearer et al., 2006;Shelly et al., 2016), our method does not assume that co-located events are located on the same fault patch, share similar fault rupture patterns, or have similar recorded waveforms at the same station. The reason is that many physical mechanisms can cause different earthquake ruptures at the same location, like local fault geometry complexity, stress heterogeneity, and changing slip motions on the same fault patch with time. We first calculate the amplitude ratios using the P-wave and S-wave amplitudes that are measured from each individual earthquake waveform instead of waveform cross-correlations. Second, we considered the amplitude ratio differences caused by the differences in focal mechanisms and magnitudes of the event pair (Equation 4). These inter-event differences are not noise but the key information in constraining focal mechanisms in the REFOC algorithm. Although the algorithm can help to obtain many more focal mechanisms in the area with highly clustered earthquakes, the method does not have any preference for one type of event over another. Therefore, the focal mechanism catalog solved using the REFOC algorithm does not have an artificially high number of a certain type of focal mechanisms than others in the same location. It has the potential to improve the resolution of stress inversions without introducing aliases.
Instead of performing focal mechanism inversion using all inputs together with different weights, we sequentially add different types of inputs to reduce the number of acceptable solutions and further reduce the uncertainties of final solutions (Hardebeck & Shearer, 2002). There are several reasons for doing so. First of all, there are many uncertainties in the input data, like the errors in polarity picking and phase picking, the noise in earthquake waveforms, and the inaccuracy in the velocity models and locations. Instead of solving for a single solution, we perform multiple trials using different velocity models and locations, and extract all acceptable focal mechanisms with up to a certain fraction of misfits. These solutions can be used to estimate the final solution and the associated uncertainty. Second, different inputs are added sequentially to the focal mechanism calculation because they provide different constraints: amplitude ratios only constrain the nodal plane orientations while polarities help to determine the sign of the strain (compressional or tensional). Large polarity misfits may result in highly inaccurate or even opposite focal mechanism solutions. To obtain accurate focal mechanism solutions, it is important to use the polarity to determine the focal mechanism first and use amplitude ratios to constrain the nodal plane orientations later. Furthermore, there is no clear practical way to determine the weights of different constraints and different weights may result in highly different solutions. Therefore, the inversion method used in the REFOC algorithm can provide stable estimation of the focal mechanism solutions and their uncertainties without introducing many hyperparameters.
The focal mechanism solutions can be further constrained by performing the second round calculation for multiple iterations. The reason we only use one iteration with inter-event relative amplitude ratios is that they only play a secondary role in constraining the focal mechanism solutions compared with the polarities and S-/P-wave amplitude ratios. The synthetic inter-event amplitude ratios are not only related to focal mechanism differences but also affected by the assumed velocity models. Due to the uncertainties in the estimated take-off angles and measured amplitude ratios, iteratively resolving the focal mechanisms may result in data overfitting rather than obtaining better solutions. Therefore, we only use the inter-event amplitude ratios. In contrast, for many other source parameters with limited absolute constraints, like relative earthquake locations and stress drops, it is important to perform many iterations of calculations to better decompose the source, path, and site effects from the seismic records and better estimate the source properties Waldhauser & Ellsworth, 2002). The source inversion methods using inter-event relative information can significantly enhance the resolution of source properties in swarms and aftershock sequences and has the potential to provide important insights into many related physical processes.
The solved high-resolution focal mechanisms can further be used as inputs to constrain other source properties. Large earthquakes are usually modeled as a fault plane rupture for solving their source dimensions, rupture velocities, and directions. These unknown parameters are too many to be well-constrained for small earthquakes because of the small number of high-quality observations and limited coverage. Therefore, small earthquakes are usually modeled as a point source to estimate location and source dimension Waldhauser & Ellsworth, 2002) and as a line-source to estimate rupture velocity and direction (Abercrombie et al., 2017). The simplified source models and fewer parameters provide more robust estimations of source properties. Meanwhile, the solved properties do not have any fault plane information and are less interpretable. The missing fault plane information can be provided by the focal mechanisms with nodal plane uncertainties less than 25° (>40% of the catalog events) solved using the REFOC method. These nodal planes can be used as additional inputs for estimating rupture velocities and directions, which can both keep the 3D-plane-source assumption and provide robust estimation of source properties. Therefore, it is promising to use high-quality focal mechanisms to robustly estimate many source properties under plane-source assumption.

Implications of the Focal Mechanism Properties in the Parkfield Area
Unlike the highly consistent focal mechanisms for M > 4 earthquakes and in the locked section, focal mechanisms in the creeping section are highly heterogeneous (Figures 10c and 10d) with high angular difference from the M6.0 mainshock (Figure 9), suggesting that most seismic deformations in the creeping section occur along the faults that are not parallel to the main fault orientations and do not represent the deformation along the main fault. The heterogeneous focal mechanisms in the creeping section might be caused by the following reasons. The main fault in the creeping section has low frictional strengths (μ = ∼0.1) and velocity strengthening friction behavior (Carpenter et al., 2015) so that the main fault in the creeping section is mainly dominated by aseismic deformation and the surrounding crustal volume is more likely to host seismic deformation. Due to low fault coupling (Jolivet et al., 2015), there is limited shear stress accumulated in the creeping section, which means that the right-lateral strike-slip faulting along the main fault is not the preferred faulting type. Moreover, the large number of reverse faulting events NW to the MM (Figure 11) suggests that the slightly transpressional bend near MM may also affect the focal mechanism properties in the creeping section. As a result, most small earthquakes in the creeping section may occur on secondary faults and cannot represent the seismic behaviors of main fault.
Temporally, the focal mechanism distribution changes significantly after the 2004 M w 6.0 Parkfield mainshock, suggesting that the focal mechanism properties reflect the evolution of stress state along the fault. Before the M w 6.0 mainshock, the shear stress is highly concentrated in the locked section compared with the other areas, causing highly consistent right-lateral strike-slip focal mechanisms. After the M w 6.0 mainshock, the coseismic slip leads to stress drops in the locked section but stress concentrations in the transition section, which is consistent with the observed increasing angular differences from the M w 6.0 mainshock in the locked section but decreasing differences in the transition section. After 5 years, the focal mechanism properties in both transition and locked sections show slight recovery to the pre-mainshock period, indicating stress recovery in both sections with possible post-seismic changes of tectonic stress direction. The focal mechanism properties in the creeping section do not show significant variations related to the mainshock rupture because the stress perturbations caused by major earthquakes are quickly accommodated by the aseismic deformation in the creeping section. Therefore, it is promising to use the focal mechanisms of small earthquakes to track tiny stress changes in the crust and help the understanding of seismic behavior in different stages of earthquake cycle.
The variations of earthquake faulting style also have important implications to the fine-scale fault structures and the major rupture slip behaviors. The correlation between surface fault traces and focal mechanisms at depth suggests that the observed restraining bend and extensional jog at the surface may extend deep into the crust and cause the observed changes of faulting type. However, fault geometry complexity can only explain the observed spatial variations. The increasing number of normal and revere faulting events after the mainshock may suggest the effect of rupture directivity of the 2004 M w 6.0 mainshock, which affect the mass movement on either side of the fault (Titus et al., 2011) and hence affect the earthquake faulting type. The 2004 M w 6.0 mainshock is dominantly unilateral to the NW with a hypocenter near the SE end of the rupture (Kim & Dreger, 2008) which can cause mass and stress concentration to the NW and mass deficit to the SE. The effect of mainshock rupture directivity on the mass movement can enhance the compressional stress to the NW and extensional stress to the SE and cause the variations of faulting type with time. Therefore, faulting type variations have the potential to provide information on fine-scale fault structure at depth and the directivity of main fault slip.
Based on the above-mentioned observations, earthquake focal mechanisms solved using the REFOC algorithm delineate unprecedented details about subtle changes of focal mechanism properties in the Parkfield area, where the events were always assumed to be on the main fault (Thurber et al., 2006). The spatiotemporal variations of focal mechanism properties provide useful constraints on many important scientific questions, like the relationship between seismic/aseismic deformation, fault strength, small-scale fault structures at depth, stress evolution along the fault, as well as the directivity effect of fault slips. The observed focal mechanism properties can also be combined with other observations to better constrain and understand the fault and crustal properties. The analysis done in this study in the Parkfield area is particularly valuable for the fault monitoring in northern and central SAF. This region has many large-scale localized faults, like the Calaveras Fault, Hayward Fault, San Andreas Fault, and San Gregorio Fault, with large historic earthquakes and many small earthquakes occurring frequently around them (Chaussard et al., 2015;Harris et al., 2021). The small earthquakes' focal mechanism properties provide a means to image and monitor these major faults' structure, strength, and stress in high resolution, providing important clues for the fault rupture physics, seismotectonic activity, as well as the time-dependent seismic hazard assessment.