## Abstract

Large earthquakes often lead to transient deformation and enhanced seismic activity, with their fastest evolution occurring at the early, ephemeral post-rupture period. Here, we investigate this elusive phase using geophysical observations from the 2004 moment magnitude 6.0 Parkfield, California, earthquake. We image continuously evolving afterslip, along with aftershocks, on the San Andreas fault over a minutes-to-days postseismic time span. Our results reveal a multistage scenario, including immediate onset of afterslip following tens-of-seconds-long coseismic shaking, short-lived slip reversals within minutes, expanding afterslip within hours, and slip migration between subparallel fault strands within days. The early afterslip and associated stress changes appear synchronized with local aftershock rates, with increasing afterslip often preceding larger aftershocks, suggesting the control of afterslip on fine-scale aftershock behavior. We interpret complex shallow processes as dynamic signatures of a three-dimensional fault-zone structure. These findings highlight important roles of aseismic source processes and structural factors in seismicity evolution, offering potential prospects for improving aftershock forecasts.

## INTRODUCTION

Observations of crustal deformation and seismicity associated with major earthquakes are increasingly common, often with unprecedented details, e.g., (*1*–*6*). It is well known that many large earthquakes are followed by vibrant aftershock sequences and transient, longer-term deformation, both of which release and redistribute stresses over a rapidly evolving post-rupture period. As our experience of characterizing the aftermaths of major seismic events continue to grow, our knowledge of how earthquake ruptures transition to post-rupture processes remains limited. Long-standing challenges include the scarcity of on-scale broadband observations and the high instrumental noise and confounding signals of multiple source processes during this ephemeral phase (*7*). Nonetheless, this transition affects estimates of slip budget and moment deficit in the long-term seismic hazard assessment and remains a critical phase for us to understand what determines the eventual size of an earthquake rupture and the style of faulting under a wide range of strain rates.

The Parkfield section of the San Andreas fault (SAF) in central California has long provided a natural laboratory for studies of seismic and aseismic fault zone processes, owing to the dense, continuous regional monitoring network. This fault segment, known for hosting recurring magnitude ~6 earthquakes (*8*), connects the northern creeping section and southern locked section that historically produced large, destructive earthquakes over the regional scale (*9*). Parkfield and nearby fault sections also host a diverse spectrum of smaller-scale fault slip phenomena, including frequent occurrences of characteristically repeating microseismicity (*10*), along with tectonic tremor and slow slip events below the conventional seismogenic layer (*11*, *12*). The most recent moment magnitude (*M*_{w}) 6.0 earthquake at Parkfield in 2004 was documented by a diverse range of geophysical instruments (Fig. 1). The rupture evolution during the earthquake and aftershock sequences explored in seismological studies (*8*, *13*, *14*) are complemented by longer-term interseismic, transient, and postseismic processes investigated with geodetic observations (*15*–*17*). The complexities in these seismic-cycle processes shed light on the mechanisms of crustal faulting and underpin the development of physics-based models of earthquakes with predictive power (*18*). The conditions and behavior of the Parkfield segment are also important for mediating the interaction between large neighboring SAF segments that can influence regional and potentially system-level seismic hazard in California (*19*, *20*), with implications for other major fault systems.

Here, we take advantage of the near-field observations from the 2004 Parkfield earthquake to study the early post-rupture dynamics of the SAF, focusing on two distinct processes, afterslip (aseismic fault slip) and aftershocks. We use a combination of high-rate (1 Hz) and daily GPS data, aided with seismological data, to retrieve the surface displacement records that span the complete coseismic-to-postseismic periods since the earthquake initiation. We use GPS-derived coseismic and early postseismic fault slip evolution to characterize the relation between afterslip and aftershock activity and investigate the rheological and structural conditions of this fault section, as well as potential mechanisms of aftershock generation.

## RESULTS

### Distinguishing coseismic and postseismic phases

We use high-rate GPS displacement time series, in tandem with trigger-mode strong-motion accelerations, to first identify the coseismic period when near-source ground shaking persists. The time-varying ground displacements settle into a steadier level over a period of 1 to 2 min after the earthquake initiation time [17:15:24 universal time coordinated (UTC), 28 September 2004] (Fig. 2A). Visual examination of the accelerometer waveforms suggests a shorter duration of ground shaking for up to 30 s at each seismic station, whereas a spectrum analysis indicates that high-frequency GPS signals fade at ~60 s (Materials and Methods; figs. S1 and S3). Taking into account the dynamic ground displacement and additional time required for the passage of seismic waves at all stations, we set the time of 90 s after the initiation of this earthquake as the division between the coseismic and postseismic phases throughout this paper. The assumption is that the coseismic phase is somewhat longer than the actual duration of dynamic rupture of the fault zone. Next, we merge processed high-rate and daily GPS solutions in a consistent manner to form continuous displacement records for the coseismic-to-postseismic periods from seconds to years, spanning time scales of more than seven orders of magnitude (Materials and Methods; figs. S7 to S9).

The overall postseismic signals are characterized by rates that largely decay logarithmically with time, which motivates us to focus on postseismic displacement history on a logarithmic scale (Fig. 2, B and C). We project the east and north components of displacement data to fault-parallel and fault-perpendicular directions, using a strike angle of 319° (*14*), to help us distinguish GPS signals due to physical processes from noise. The fault-parallel components exhibit a largely monotonic trend, except for station CARH, which records near-source complexity as we will describe later (Fig. 2B). The fault-perpendicular components show only small displacements over the long term, albeit with notable variability at various periods, e.g., 2 to 10 min and 3 to 5 hours. After assessing the site- and satellite-related conditions, we conclude that stations PKDB and MIDA have higher noise levels than other stations, with issues at the very early postseismic time, and that ~3 to 5 hours after the earthquake is an anomalous period for the entire network, as we will describe below (Materials and Methods).

### Continuous evolution of surface displacement field

The surface displacement field shows distinct patterns over the coseismic and several postseismic intervals (10 min, 100 min, 17 hours, and 7 days) (Fig. 3A). Comparisons of fault-parallel GPS displacements reveal characteristic profiles of fault motion at different time steps, including the changing location of the primary slip zone (Fig. 3B). While the coseismic rupture clearly occurs between stations CARH and POMM, the location of the fault plane hosting postseismic slip varies with time, moving from the southwestern side of CARH to its northeastern side within days after the earthquake. The different locations are associated with the surface traces of the Southwestern Fracture Zone (SWFZ) and the main SAF, which are suggested to merge into a single fault at a depth of ~6 km based on aftershock clusters (Fig. 3C) (*8*). The switching of slip plane from SWFZ to SAF is suggested to occur at around 3 to 5 days based on GPS data (Fig. 3C) (*16*).

The use of high-rate, as opposed to daily, GPS displacements directly affects the estimates on the partitioning of coseismic and postseismic signals (*21*, *22*). The high-rate solutions record large aseismic signals that would be indiscernible from coseismic signals in the daily solutions (Fig. 4). For the Parkfield event, the coseismic displacements estimated from the high-rate solutions (0 to 90 s) are around 30 to 70% of the coseismic offsets based on daily solutions before and after the event (Fig. 4). This ratio is generally smaller for stations closer to the fault, indicating more substantial contributions from early postseismic processes. A notable exception to this trend is station CARH, located between SWFZ and SAF, due to the complex changing directions at this site in response to processes below. Using postseismic displacement records, we also estimate characteristic times in a logarithmic or exponential function (*23*) with different observational windows (from 2 min to a maximum of 14 years) (fig. S10). The derived effective characteristic times are not only variable among stations but also dependent on the chosen time windows, demonstrating that the postseismic data cannot be reduced to a simple logarithmic or exponential function with fixed model parameters, as noted in (*24*).

These near-source features—including the switch of primary slipping fault, spatiotemporal variations in coseismic-to-postseismic signals, and location- and time-dependent characteristic time scales—indicate the complexity of source processes captured by the near-field GPS network. Using finite-fault inversion, we explore to what extent subsurface fault slip processes can reproduce these observations. To account for the change in involved fault planes, we construct fault geometries for both the SWFZ and SAF and use them in fault slip modeling for the early postseismic phase before and after ~4 days, respectively (Materials and Methods; figs. S11 and S12). We also incorporate measurements of local surface slip from one creepmeter installed on SWFZ and four creepmeters on SAF (Materials and Methods; fig. S1).

### Coseismic rupture and postseismic slip evolution

We adopt a Bayesian formulation of the inverse problem to obtain closed-form solutions for the posterior probability distribution of fault slip models, subject to relatively weak a priori constraints on model and data misfit metrics that account for both observational and modeling uncertainties (Materials and Methods). We first perform inversions on synthetic scenarios of fault slip to test the resolution of our imaging approach and guide our choices of inversion parameters and interpretation of source models (fig. S13). As expected, the resolution of peak slip decreases with depth, especially for more localized slip sources, whereas the potency (the product of the subfault area and the slip amount) can be better resolved (Materials and Methods).

We infer the spatial distribution of coseismic and cumulative postseismic slip on the SWFZ using the posterior mean models constrained by surface displacement data (Fig. 5 and fig. S14). The coseismic model features two main slip patches near the hypocenter and further to the northwest, with the most potency release around the depths of ~7 km (Fig. 5A). The upper limit of the largest slip patch locally extends to shallow depths above 4 km, whereas the streaks of aftershocks later appear at the depths of ~5 km (Fig. 5D). The surface displacements are well reproduced by the source model (Fig. 5F). A few subfault patches near the surface may contain localized slip with varying slip directions, observed similarly in postseismic models. This is likely due to the need to fit near-fault stations, where our fault geometry may oversimplify the complex shallow fault zone indicated by the discontinuous surface traces. The actual rupture processes may contain more localized slip patches, as suggested by strong motion–based source models (*13*, *25*), but the overall potency should be better constrained by geodetic displacements.

Acknowledging that geodetic data primarily constrain seismic potency rather than moment, we estimate that our coseismic slip model has a moment of ~1.7 × 10^{18} N ∙ m, equivalent to a moment magnitude of 6.09. Our postseismic models at 34 min, 24 hours, and 15.4 days have a cumulative moment that is about 7, 34, and 100% of the coseismic moment release, respectively (fig. S15).

The surface displacements that rapidly vary over the early postseismic period can be explained by evolving afterslip on the fault for the most part, from which we select several representative epochs (Fig. 5 and fig. S14). As early as ~2 min after the coseismic rupture (3- to 4-min postseismic time), near-field signals start to appear at several stations (most evidently at HOGS and CARH and probably at HUNT and POMM), reaching a peak value at a postseismic time of ~5 min (fig. S14A). The predominantly fault-perpendicular motion up to ~5 min can be reproduced by a localized shallow slip patch of ~3 km deep beneath the network (Fig. 5B). This displacement pattern then decays, disappears before ~10 min (fig. S14B), and reemerges afterward with growing fault-parallel movements (Fig. 5C and fig. S14C). Postseismic signals at several later time steps, 34 min, 8.5, and 24 hours, can be reproduced by afterslip on the SWFZ that increases in amplitude and expands in space, both laterally and downdip (Fig. 5, C to E). After ~4 days, the reversal of displacement vectors at station CARH requires that the ensuing afterslip takes place on the main SAF (fig. S14, F and G). Despite larger uncertainty in local slip amplitude, the depth profiles of coseismic and postseismic slip at different stages (Fig. 5, H to L) suggest that coseismic slip extends above the shallow seismicity streaks, and the centroid of cumulative postseismic slip moves downdip with time.

For these snapshots of fault slip distribution, we also show the spatial distribution of background seismicity and aftershocks that occur in corresponding postseismic periods. In this and later comparisons, we consider two high-resolution earthquake catalogs, a relocated catalog from the Northern California Seismic Network (NCSN) (*26*), and a 3-day aftershock catalog that is more complete for small events (hereinafter referred to as “PZ09”) (*27*). The aftershocks in PZ09 catalog are relatively scarce before 34 min and form more pronounced streaks in later times of the first day.

The postseismic time of ~3 to 5 hours appear to be an “anomalous period” for GPS observations, where surface displacements exhibit back-and-forth fault-perpendicular movements that, at times, even exceed fault-parallel movements at several stations (Fig. 2). Although some slip models can reasonably reproduce the displacement data, these models appear to contain notable backward slip and large variation in local slip vectors (fig. S14C). While we cannot exclude that some variability in the data is due to anomalously higher GPS noise levels, the large nonmonotonic signals observed in the fault-perpendicular component may suggest other fault zone processes or structure that our simplified dislocation models cannot capture (Materials and Methods). In the movie S1, we show our time-dependent models of afterslip at finer time steps over an extended postseismic period, along with the metrics for data misfit, which also appears anomalously larger for this period.

### Relationship between early afterslip and aftershocks

The early afterslip process constrained by high-rate geodesy complements seismological observations of early aftershocks immediately following the Parkfield earthquake (*14*, *27*). These independent observations and inference of aseismic and seismic moment release provide a unique opportunity to examine the relationship between the two physical processes over early postseismic time scales. We compare the temporal behavior of afterslip and aftershocks, along with stress change (described later), within three larger fault areas at shallower (1.5 to 6 km) and greater depths (8.25 to 12.75 km) (Fig. 6). The selected regions are either directly above (R1) or below (R2) the location of coseismic peak slip, with implications for the depth-dependent properties of a seismogenic fault, or in adjacent areas to the northwest (R3) where the occurrence of two magnitude 5 (M5) aftershocks on 29 and 30 September is of great interest.

We first focus on the temporal behavior of afterslip and aftershocks in the three regions (Fig. 7). We consider the history of area-averaged slip (proportional to potency) within these regions to reduce variability and uncertainty of inferred slip history on individual subfault patches, along with the estimates of slip rates shown on a logarithmic scale (Fig. 7, A to C). For aftershocks, we compute time-dependent moment rate *M*(*t*) and seismicity rate *R*(*t*) for “near-fault” aftershocks (within a distance of 3 km of the fault plane) with a cutoff magnitude of 1.5 in the corresponding periods (Materials and Methods). We compare *M*(*t*) (Fig. 7, D to F) and *R*(*t*) (Fig. 7, G to I) with the amplitude and rate of cumulative afterslip on a logarithmic time scale. Since the NCSN catalog clearly misses smaller events throughout early postseismic period compared to the PZ09 catalog, we perform further analysis of seismicity using only the latter.

Over shallower fault areas (R1), the cumulative fault slip largely follows a logarithmic-like postseismic response, starting a few minutes after the earthquake. However, a notable deviation from the monotonic increasing trend occurs between 5 and 15 min, where the increased slip of ~30 mm reverses its direction before slipping in the right-lateral direction again at ~15 min (Fig. 7A). We consider this process a potentially short-lived reversal of early afterslip. The total moment rate of aftershocks is more variable, whereas the seismicity rate may be expected to follow Omori-like laws. We fit two functions, an Omori-type law *R*(*t*) = *a* ∙ *t*^{−p} (*28*, *29*) and a logarithmic function *R*(*t*) = *a* ∙ log (*t*) + *b*, for the monotonic decreasing trend of seismicity rates over the period from 7000 to 197,000 s (the later time span of PZ09 catalog) and extrapolate the two functions to the earlier postseismic period where they diverge. Aftershocks within the first 2-min time window can be plausibly attributed to dynamic stress perturbations from the mainshock rupture (Fig. 7D). The productivity rate of early aftershocks after 2 min is deficient compared to the Omori law prediction and fits the logarithmic pattern more closely, with the temporal variability probably due to noise and small total event numbers. It is worth noting that the initial moment release rate of afterslip (average slip rates of 10^{−5} to 10^{−4} m/s over an area of 10 km × 4.5 km) is on the order of 10^{13} to 10^{14} N ∙ m/s, much higher than the seismic moment release rate of aftershocks on the order of 10^{11} N ∙ m/s within the same area.

At greater depths below the coseismically ruptured patch (R2), the cumulative fault slip remains insignificant up until at least ~30 min after the earthquake rupture (Fig. 7B). The inferred slip appears to first increase and then decrease from ~1.5 to 6 hours, corresponding to the “anomalous period” of surface displacements that we identified earlier. It seems unusual and physically implausible that substantial left-lateral slip can occur at depth during that period. We speculate that if the higher GPS noise and oscillating fault-perpendicular movement can be accounted for, the underlying afterslip process, which is presumably continuous in time, could follow a simpler trend. Regardless, our results suggest delayed occurrence of deeper afterslip within hours. The rates of deeper aftershocks show a delayed increase at ~20 min and a decrease within the first hour.

The deeper fault areas further to the northwest (R3) show slightly different behavior compared to R2 (Fig. 7C). While more aftershocks occur in the period with seismic waves, the most intense aftershock activities occur after 10 min and are often associated with larger aftershocks, including two M5 aftershocks. In addition, the enhanced moment release rates of aftershocks are preceded by the initiating or ongoing afterslip. The inferred afterslip history similarly contains complex evolution during the “anomalous period” of ~3 to 5 hours.

We further compare local stress change due to afterslip with aftershock behavior. Using posterior models of fault slip, we compute the Coulomb stress change history within these regions, taking into account the full uncertainty of slip models (Fig. 7, J to O; Materials and Methods). The coseismic rupture increases stress near the surface and decreases stress over deeper ruptured areas (Fig. 5P); then the evolving afterslip continuously redistributes stress, resulting in a large spatial variability in local stress change history. We note that stress change during the coseismic and postseismic stages have large uncertainties, with their area-averaged amplitudes depending on the location and size of selected regions.

Focusing on the shallow fault areas that contain the streaks of aftershocks, we find that the very early temporal trend in the stress change history (i.e., the sign of stressing rates) is the most robust feature for different selections of region R1 that enclose the aftershock streaks (Fig. 6 and S18). The initial negative stress change at 3 to 5 min and later positive stress change between 5 and 10 min (Fig. 6, B and C), with stressing rates on the order of ~100 Pa/s, can be directly attributed to the initial afterslip episode and ensuing slip reversal. These stress changes due to afterslip evolution coincide with the notable decrease and later increase of seismicity rates at the same time. However, stressing rate variation after 10 min are more sensitive to the selected regions and much reduced in amplitude due to rapidly decaying slip rates, which preclude us from making further robust, meaningful comparisons. Although stressing rates in R2 and R3 are less constrained, the dearth of aftershocks before 30 min in R2 and 15 min in R3 suggests different aftershock-triggering conditions of the deeper fault zones.

## DISCUSSION

A diverse range of geophysical studies over the years have contributed to illuminating source processes of the Parkfield earthquake sequence, e.g., (*13*, *15*, *16*, *18*, *24*, *25*, *30*). Our coseismic slip model is largely consistent with models derived from high-rate GPS data (*30*, *31*) and seismological data (*13*, *25*) in terms of the lateral locations of large-slip areas and overall moment, although different resolutions between seismic and geodetic source models and different spatial smoothing applied in geodetic models preclude a detailed comparison. Our longer-term postseismic slip models are also consistent with geodetic studies that focus on postseismic periods of days or longer (*15*, *16*, *24*, *32*). To our knowledge, the evolution of coseismic-to-postseismic fault slip processes from minutes to days has not been explored so far, and our study fills an important observation gap.

The coseismic surface displacement field indicates a dynamic rupture process that involves the SWFZ near the surface, rather than the shallow portion of the SAF, as corroborated by field surveys of immediate surface slip and creepmeter measurements (*33*). It is likely that the earthquake at least partially ruptured updip beyond the shallow microseismicity streaks at the depth of ~5 km, indicators of subsurface geometric connectivity of SWFZ and SAF, due to large dynamic stress changes on the shallow low-strength fault areas. A coseismic rupture through geometric complexity or heterogeneous stress fields may generate the observed high-frequency seismic radiators during the event (*34*), which appears to coincide with the location of the highest spatial gradient of coseismic slip, and higher stress change in our model. Penetration of coseismic slip into nominally creeping regions is also consistent with the low activity of nearby aftershocks immediately after the mainshock (*35*, *36*).

The postseismic displacement field exhibits coherent, small-scale signals that occur shortly after the earthquake arrest and are geodetically resolvable as early as ~3 min after the earthquake initiation time. We have explained these very early signals by a localized shallow source of afterslip on the same fault plane as earthquake rupture, which is consistent with later measurements from the creepmeter on SWFZ (*24*, *33*). During the early postseismic period, the adjustment of pore fluid flow after coseismic disturbance to the pore pressure may be a possible cause for ground displacements, often over time scales of weeks to years for some major earthquakes, e.g., (*37*, *38*). To reproduce early postseismic signal over several kilometers within minutes at Parkfield would require much higher hydraulic diffusivity than in the previously documented scenarios. We estimate the expected displacement field of such poroelastic rebound and find that the predicted displacement pattern does not explain the directions and localization of the observed signals (Material and Methods; fig. S17). Viscoelastic relaxation of the lower crust is also shown to be negligible within a short postseismic time window (*39*). Therefore, afterslip is likely the predominant and sufficient mechanism for surface displacements (*40*) over all early postseismic time scales.

The early afterslip after the Parkfield earthquake shares a considerable fraction of the total moment release within the first day, a phenomenon noted by previous studies in this region (*16*, *24*) and similar to recent observations from several subduction zone earthquakes, e.g., (*22*, *41*). With the high-rate GPS data at Parkfield, we also find possible occurrences of short-lived left-lateral slip that is geodetically resolvable around 5 to 10 min. This inference is ultimately based on the cumulative movements in the fault-perpendicular direction most clearly seen at HOGS (~8 mm), CARH (~5 mm), and possibly POMM at 5 min and vanishing of these cumulative displacements between 5 and 10 min. The distinct motion in the directions parallel versus perpendicular to the fault plane and good signal-to-noise ratio at the involved stations lead us to conclude that such subtle localized signals are related to real fault-zone processes, rather than random noise or data artifacts. Our models of a shallow afterslip source with primarily right-lateral slip and later left-lateral slip above the coseismic rupture area can explain these fault-perpendicular movements at the earliest postseismic period.

Reversed fault slip is physically plausible for low-strength faults undergoing major adjustment of stress state, with some observational evidence from reversed surface creep under earthquake-induced stress perturbations (*42*) or reversed aftershock mechanisms (*43*, *44*). On the shallow SWFZ, transient left-lateral slip may occur near the boundary of coseismically ruptured areas, possibly as a post-rupture response to dynamic stress overshoot in a fault zone with subsurface geometrical and material heterogeneity. The immediate afterslip is likely driven by rupture-induced stress elevation at shallow depths. Since the resolution of our models only allows us to infer area-averaged potency release over spatial scales of ~4 km, our inferred afterslip and slip reversal episodes may plausibly result from a combination of two processes that occur over different fault areas with different time evolution.

Over later periods, the shallow afterslip continues to grow in amplitude and expand in space well into 1 month and longer after the earthquake, consistent with observations from longer periods (*39*, *40*, *45*). In the framework of the rate-and-state friction laws, afterslip can be explained by velocity-strengthening frictional behavior (*46*, *47*), with the initial slip rates inversely proportional to the product of frictional parameters and effective normal stress, (*a* − *b*)σ, for a given stress change. Considering the overall power law–like decay of shallow afterslip rates, with the only exception for the first ~10 min, we expect that the effective (*a* − *b*)σ derived from earlier postseismic periods to be largely consistent with estimates from longer-term postseismic records, e.g., (*31*, *32*). We do not quantitatively estimate these parameters due to large uncertainties in coseismic stress perturbations and our focus on the qualitative differences between postseismic periods and between fault areas. In contrast to the shallow region, the deeper afterslip in both R2 and R3 shows a delayed or weak start, with the cumulative slip exhibiting a non-monotonic trend around the anomalous period (~3–5 hr) before resuming a steady logarithmic increase. The early pattern suggests an apparent downdip migration of afterslip and can be attributed to larger (*a* − *b*)σ and its potential change with time.

The concurrent evolution of afterslip and aftershocks have been observed in some subduction zones and crustal faults, typically over time scales longer than hours or spatial scales larger than tens of kilometers (*1*, *41*, *48*, *49*). Questions remain on whether the two processes reflect similar or common statistical fault properties, physical mechanisms, or both. Our results show a robust correlation between the decrease and increase of earliest shallow aftershock rates and of stress change history that can be directly attributed to the initial episodes of afterslip and slip reversal. These independent observations lead us to conclude that early afterslip, with a larger moment release rate than that of aftershocks by two to three orders of magnitude, likely modulates aftershock activity at fine spatial scales and across early postseismic time scales.

The rise and fall of seismicity rates within a period of overall reduced early aftershock activity can be explained by negative stress changes instantaneously incurred by coseismic rupture and time-dependent stressing from nearby afterslip episodes, providing support to the interpretation of migrating early aftershocks driven by afterslip (*27*). We normally do not expect the time evolution of afterslip-induced stress change and seismicity rates to be a simple correlating function due to the mainshock-induced stress perturbation and nonlinear dynamics of rock friction (*46*). In our case, the apparent correlation and physical connection between stress change due to early afterslip episode and seismicity rates seem more plausible because of the reversal of slip direction that result in a flipped sign of stress change (Fig. 8). After 1 hour, both seismicity rates and afterslip rates follow monotonic power law decay, consistent with theoretical models (*50*, *51*).

Time-varying afterslip may also play a role in the generation of larger aftershocks following the Parkfield earthquake. In R1 and R2, periods of higher moment release rates relative to the seismicity rate usually follow the onset of afterslip (Fig. 8, A and B). In the deeper part closer to the creeping SAF faults (R1), ongoing afterslip, with minor stress increase, appears to precede episodes of higher moment release, including two M5 aftershocks on 29 and 30 September. These observations are consistent with a scenario that afterslip affects not only small aftershocks but also some larger aftershocks in the region.

We interpret the complexity in the shallow seismic-aseismic processes over the earliest postseismic period as the dynamic post-rupture behavior of three-dimensional (3D) fault zone structure at Parkfield (Fig. 8). Despite the recurrence of similar M6 earthquakes, the Parkfield section of the SAF features more complex fault geometry than nearby segments, with evident discontinuous fault strands at the surface and likely intersecting subsurface fault planes, reflecting persistent structural and mechanical heterogeneity of this transitional segment subject to the differential motion of the adjacent locked and creeping SAF faults (*52*). Although the coseismic rupture had preferentially ruptured the SWFZ, the longer-term plate-boundary motion is accommodated by the main SAF fault, requiring a transfer of motion between the subparallel fault strands. We suggest that the shallow afterslip initiated immediately above the largest coseismic slip patch (Fig. 8C), followed or accompanied by slip reversals over nearby fault areas within minutes, before expanding within hours and gradually migrating from SWFZ to SAF over days. After the coseismic rupture, the evolving afterslip redistributes stress and influences the behavior of shallow aftershocks that are indicators of structural heterogeneity of the shallow fault zone. Therefore, observations of early post-rupture phase can provide insights about the roles of aseismic source processes and structural factors in the production and evolution of seismicity.

Our source imaging approach incorporates relatively weak constraints on the amplitude and rake angles of fault slip, thereby favoring the weight of observations in the inference of qualitative source features. Because of large uncertainties in our source models, we have not yet established a quantitative relationship between aseismic slip, stress change, and seismicity. Improving source imaging with additional constraints could help us test and develop physical models of seismicity. Furthermore, the unclear origin of additional complexity in our observations, such as the fault-perpendicular surface displacements over the identified anomalous period (3 to 5 hours), requires a future investigation of high-rate GPS noise mitigation and other fault-zone processes.

In summary, we use high-resolution geophysical observations to provide a more complete picture of the early post-rupture processes on the SAF following the 2004 M6 Parkfield earthquake. Our results reveal closely connected afterslip and aftershock processes over early postseismic time scales and synchronous complexity in both processes that may point to dynamic signatures of 3D fault zone structure at Parkfield. While the combination of certain source and structural factors may be unique for this SAF fault segment, our study demonstrates the values of expanding our observational window, aided with robust source imaging approaches and uncertainty characterization, to better resolve detailed aseismic and seismic processes under varying strain rates. These efforts could have important implications for understanding earthquake triggering; improving aftershock forecasts; and estimating the slip budget, timing, and patterns of major earthquakes. Other major fault zones around the world would benefit from long-term seismological and geodetic monitoring networks with near-source coverage, which offer promising opportunities to probe the dynamics of fault zones, major earthquakes, microseismicity, and aseismic transients over multiple time scales.

## MATERIALS AND METHODS

### Strong-motion records

Strong-motion records from 60 stations during the earthquake were retrieved from data archives of the National Strong-Motion Project and the California Geological Survey (*53*). Recording of these acceleration seismograms was triggered when the acceleration amplitude reached a threshold value of 0.005 g at most stations. The horizontal components of accelerographs show that ground shaking lasted for up to 30 s (fig. S1, A and B), which is shorter than the dynamic displacement period of ~1 to 2 min that we estimate from the high-rate GPS data (Fig. 2). Because of potential errors in integrating accelerographs without colocated permanent offset measurements, we only use the seismic records to assess ground shaking and wavefield complexity, rather than for fault slip modeling.

### High-rate GPS solutions

A high-rate real-time GPS network was installed around the Parkfield area in July 2003, consisting of 12 near-fault GPS stations and one reference station (CRBT), about 40 km southwest of the SAF (Fig. 1) (*54*). The 1-Hz GPS dual-frequency phase and pseudorange observations were processed over 2 weeks from 22 September to 5 October 2004 to span the *M*_{w} 6.0 Parkfield earthquake on 28 September (17:15:24 UTC). We used the method of instantaneous positioning (*55*) to estimate displacements relative to station CRBT (fig. S2) using precise satellite ephemerides from the International GNSS (Global Navigation Satellite System) Service (IGS) while estimating single-epoch troposphere delay parameters. Examination of high-rate time series of this type exhibit a flat power spectral density from 1 to 10 Hz (*21*). Therefore, we use the typical root mean square (RMS) value from a large set of 1-Hz displacement time series as a measure of the precision of a single 1-Hz displacement, ~1 cm in the horizontal components and ~5 cm in the vertical component. To avoid a local reference station, we also estimated “absolute” displacements with respect to a global reference frame using the method of precise point positioning (*56*) but found them to be much noisier.

We applied sidereal filtering (*57*) to reduce multipath effects in the GPS displacement time series and increase the signal-to-noise ratio for estimating the permanent (static) coseismic and short-term postseismic displacements (fig. S2). To do so, we estimated the errors associated with multipath by averaging and correcting the displacement records using shifting 3-day windows offset by multiples of a sidereal day (23 hours 56 min 4 s) up to and after the event. Multipath effects are minimized during seismic shaking (*21*). Therefore, the sidereal filtering is only intended to reduce multipath effects over static periods in the frequency range of about 1 hour^{−1} to 1 min^{−1}, which removes a random walk–type noise contribution to uncover a lower-energy flicker noise signature (*21*), the latter characteristic of low-frequency GPS static positions estimated from 1–24 hours of GPS phase and pseudorange observations (next section). Sidereal filtering increases the signal-to-noise ratio by about a factor of two to three (*21*). As expected, we found no pronounced motion in the vertical direction for this strike-slip event. Therefore, only the horizontal precision is of concern.

### Daily GPS solutions

The GPS data from the 14-station Parkfield network (Fig. 1) have been processed since 2003 as part of the operational analysis of daily displacements at the Scripps Orbit and Permanent Array Center (SOPAC). The daily displacement time series with respect to the International Reference Frame 2014 (ITRF2014) (*58*) are the result of combining independent solutions by SOPAC using the GAMIT software and the Jet Propulsion Laboratory using the GipsyX software, as well as a common source of metadata (*59*). The result is a single set of 3D displacements using 24 hours of GPS phase and pseudorange data sampled at 15 s and precise satellite ephemerides from the IGS.

A principal component (PC) analysis method is applied to estimate and eliminate network-wide common mode noise from the daily displacement residuals after outliers, coseismic and noncoseismic offsets, and linear trends have been parametrically modeled (*60*). This procedure improves the signal-to-noise ratio but does not affect the signals of interest, here, the coseismic and postseismic motions. We removed the interseismic, coseismic, and postseismic signals associated with the 22 December 2003 M6.6 San Simeon earthquake from the time series *d*(*t*) according to *d*(*t*) = *Vt* + *a* + *b* log (*t*/*c*) with fitting parameters *V*, *a*, *b*, and *c*, where *V* is the interseismic linear trend and *a* is the coseismic offset (fig. S5). The postseismic motion is fit by a logarithmic function with coefficient *b* and decay time *c*. The negligible postseismic signals ~1 year after the San Simeon event justifies our treatment of the remaining signals as those only associated with the 2004 Parkfield earthquake (*45*).

Power spectrum analysis of the uncertainties of 20-year-long daily GPS displacement time series at hundreds of stations in the western United States indicates a combination of white, flicker, and random walk noise processes, e.g., (*61*). Typically, the “one-sigma” uncertainty for a single daily displacement is about ~1 to 2 mm for the horizontal components and ~5 mm for the vertical component (*59*).

### Combining high-rate and daily GPS solutions

To retrieve continuous coseismic-to-postseismic ground displacement time series, we need to merge high-rate and daily GPS solutions over minutes, hours, and days since the earthquake initiation time. First, we address the issue of different reference frames and examine the consistency of the two datasets. Since high-rate positions are relative to the far-field station CRBT and daily positions are with respect to the ITRF2014, we use the following steps to obtain pseudo “absolute” high-rate positions for the Parkfield stations (Fig. 1). We first subtract from each component of the high-rate and daily solutions the pre-event displacement as the mean value of respective solutions over 6 days before the earthquake (22 to 27 September). The daily solutions at the station CRBT also exhibit minor postseismic signals (figs. S5 and S8). To take that into account, we interpolate the daily postseismic signal at CRBT to 1 Hz and add it back to the high-rate solutions of the other Parkfield stations. After these corrections, the high-rate solutions show positions in the ITRF2014, at least for longer time scales and localized signals (e.g., associated with afterslip). The good agreement between high-rate and daily solutions (fig. S7) supports the consistency of our approach applied to the two datasets. However, the daily solution for near-fault stations clearly cannot resolve complex displacement history during the day of the event (28 September), as seen in high-rate time series. Some minor discrepancies between the two datasets still remain because of the ambiguity in estimating the pre-event displacements.

Second, we use daily GPS solutions to fill in the gap after the end of high-rate GPS solutions while ensuring a close agreement of the two datasets in the overlapping postseismic time windows (1 to 7 days after 28 September). To do that, we first apply a temporal median filter with an adaptive filtering window to both high-rate (after sidereal filtering) and daily time series to reduce overall noise levels. The filtering time steps are equally spaced on a logarithmic scale and fine enough (finer than the downsampling time steps described below) to retain prominent features in the time series, with adaptive filtering windows bracketed by logarithmic midpoints of adjacent time intervals. The original daily solutions are adjusted by a constant offset so that filtered daily and high-rate solutions give the same displacement at the first day after the earthquake, thus ensuring a better alignment in the overlapping period (fig. S8). Note that this temporal filtering only serves to facilitate this alignment and aid visualization in Fig. 2 and figs. S7 and S8. The final merged GPS time series consist of the original high-rate GPS solutions, with later periods filled in with adjusted daily GPS solutions.

### Estimating coseismic and postseismic displacements

We obtain the static (permanent) coseismic displacements as the difference between the mean high-rate GPS displacements over a 10-s pre-event period and 60- to 120-s post-event period. Data errors are estimated assuming a normal distribution. The static offset uncertainty is based on the variance of the pre-earthquake data points, which is also assumed for the postearthquake period, followed by error propagation of the variances. To obtain a smoothly varying postseismic displacement history that is suitable for time-dependent fault slip modeling, we estimate the downsampled cumulative displacements, *D*(*t*), and associated RMS errors at representative postseismic times, *t*(*k*), that are equally spaced on a logarithmic scale relative to the end of the coseismic phase (~90 s). Specifically, we use *t*(*k*) = 120 + 60 × 2^{k/2} s, ∆*t*(*k*) = 60 × 2^{k/2 − 1} s, *k* = 0,1,2, …,43, and *n* is the number of data points in 2∆*t*. The downsampled GPS data are shown in fig. S9 and provided in Supplementary Data.

### Quality control and data editing

Several GPS stations have higher noise levels than others in the network and potential instrument-related issues for certain periods. For example, MIDA is installed on a rooftop, collocated with an electronic distance measurement instrument and, hence, susceptible to more severe multipath effects. Some offsets and outliers around ~5 min in MIDA time series are removed in data processing. PKDB is the only GPS station in the network with a colocated seismometer. However, its siting is suboptimal on the lower slope of a hill obstructing the satellite view to the east. Furthermore, its sampling rate is 30 s and postseismic data within hours are noisy (Fig. 2 and fig. S2), so we exclude the first 2-hour postseismic data from PKDB in our modeling.

Postseismic displacements at almost all stations exhibit noisy oscillations over the period of ~3 to 5 hours after the earthquake (Fig. 2). To assess whether these variations are network-wide common noise or actual signals, we performed a PC analysis to identify features in the orthogonal space that explain most of the data variance (fig. S6). While the first PC captures the monotonic postseismic trend, the second PC indicates that this noisier feature around ~3 hours appears to occur throughout the network, and other PCs do not exhibit distinct time-varying signals. On one hand, the repeat of higher noise levels after 24 and 48 hours points to a possible cause of common mode errors associated with GPS satellite configurations, which were less optimal at the time of the 2004 Parkfield earthquake than those available today from multiple satellite GNSS constellations. On the other hand, the more pronounced variability in the fault-perpendicular directions suggests an influence from fault-zone processes. Correcting the second PC from data does not improve displacement data at all stations, suggesting the spatial variability of this feature within the network. Since we cannot determine the cause of this data feature, we do not apply principal components analysis–based noise mitigation to the high-rate GPS data and note the caveat that the ~3- to ~5-hour postseismic period is associated with higher noise than other periods.

### Fault geometry

The fault geometries for the SWFZ and main SAF are constrained by seismicity (mostly aftershocks) and surface fault traces (*14*, *26*, *27*), similar to a previous study that used near-field data (*16*). The SWFZ is approximated by a planar near-vertical fault, with its strike (322°) and dip (87°) angles chosen to delineate discontinuous surface traces and closely track the concentrated seismicity clusters at depths of ~6 km (fig. S11). The SAF is modeled by a nonplanar fault geometry, which closely fits the fault trace at the surface, curves downward to meet with SWFZ at depths of ~6 km, and dips vertically at greater depths (figs. S11 and S12). The SWFZ geometry (~32 km by 15 km) consists of 45 × 20 rectangular patches (0.71 km by 0.75 km) along the strike and dip directions, respectively; the nonplanar SAF geometry (~45 km by 18 km) consists of 59 × 24 rectangular patches (0.66 to 0.9 km by 0.75 km). We provide detailed fault geometry in Supplementary Data. Because of the overlap of the two geometries, we consider only one fault in the inversions of displacement data.

### Creepmeter records

Several creepmeters are installed across the SAF, and only one creepmeter (XRSW) is installed across the SWFZ (fig. S1D). These creepmeters have been operating for several years with a data sampling interval of ~10 min. Interseismic linear trends before the Parkfield event have been removed from the data (fig. S1C). Creepmeters XPK1, XVA1, and XMD1 were disrupted by coseismic shaking and contain missing records for the first day after the earthquake. Creepmeters CRR1 and TABC on the SAF show immediate signals that are followed by a more gradual increase, suggesting a combination of early instantaneously triggered surface slip due to passing seismic waves and later fault creep. The creepmeter XMBC, located farther north on the creeping section, shows several jumps after the earthquake, probably due to slip triggering or instrumental issues associated with passing seismic waves. We use the measurements of cumulative creep, when available, as a priori constrains on the slip amplitude of the nearest subfault patch. We expect that creepmeter data only constrain very localized slip at the surface, with little sensitivity to most subsurface fault slip.

### Deformation and stress change modeling

Two sets of 1D subsurface structures (multiple layers atop a half-space) are used to model crustal deformation on two sides of the SAF due to slip on the fault (table S1). The velocity structure is slightly modified from the one used in a seismological source study (*13*). The Green’s functions that map fault slip to surface displacements are computed using a propagator matrix method (*62*). We use this velocity structure to provide better constraints and reduce bias in the inferred depth distribution of fault slip.

To estimate postseismic rebound due to pore fluid flow in the upper crust, we use a coseismic slip model to compute ground deformation for the undrained scenario (immediately after coseismic disturbances to pore pressure) and the drained scenario (after pore pressure gradient dissipates through fluid flow) in a half-space with a Poisson’s ratio of 0.3 and 0.25, respectively (*37*). The difference in ground displacement of the two scenarios is the expected postseismic poroelastic rebound in a homogeneous structure (fig. S16).

To estimate stress changes that result from any fault slip model, we use the half-space solutions to compute for components of the stress tensor **Γ*** _{ij}* on the receiver fault, where

*i*,

*j*=

*x*,

*y*,

*z*(axes

*x*,

*y*, and

*z*are the strike-slip, fault-perpendicular, and dip-slip directions, respectively). We calculate the Coulomb stress

**Γ**

*=*

_{cs}**Γ**

*− μ*

_{xy}**Γ**

*, with μ = 0.6. In both cases of poroelastic rebound and stress changes, we use half-space solutions (*

_{yy}*63*) for numerical efficiency. The half-space approximation is appropriate, since we are only concerned with the qualitative effect in the former case and stress changes at seismogenic depths in the latter case.

### Bayesian inference of slip and stress change

We perform inversions of surface GPS displacements to retrieve the spatial distribution of subsurface fault slip and associated uncertainties for the coseismic or postseismic periods. The model for each epoch, of either static coseismic slip or cumulative postseismic slip, is independently derived from downsampled data. Each set of GPS displacement vector **d** is used to infer two-component slip vector **m=** [**m**_{⊥}, **m**_{∥}], (**m**_{⊥} for dip-slip and **m**_{∥} for strike-slip directions) on each subfault patch, subject to our choices on prior model covariance matrix **C**_{m} and prior mean model vector **m**_{0}. In a Bayesian framework, we formulate the prior probability distribution for model parameters, *P*(**m**), and data likelihood function, *P*(**d**∣**m**), in Gaussian forms**G** is the Green’s functions matrix and **C**_{χ} is the data misfit covariance matrix. Gaussian priors and data likelihood allow closed-form solutions for the posterior probability distribution, *P*(**m**∣**d**), also a Gaussian function (*52*)

The use of independent Gaussian functions for model priors and data likelihood can be justified by the principle of maximum entropy (*64*). We discuss how this choice affects our interpretation of models at the end of this section. The posterior mean model **S** that maps a slip model to stress change, we compute the mean stress change and covariance matrix **Γ*** _{cs}*, with 1σ credible intervals obtained from

**C**

_{Γ}.

For prior information, we choose zero mean values for both dip- and strike-slip components, i.e., **m**_{0} ** = 0**, and uniform prior model covariances based on a spatial Laplacian function

_{0}is the SD of prior models,

*D*is the distance between patches

_{ij}*i*and

*j*, and

*L*is a correlation length scale. The use of creepmeter data slightly modifies the prior constraints on individual patches, as described in the next section.

The error structure of the problem, characterized by the misfit covariance matrix **C**_{χ}** = C**_{d}*+*C_{p}, includes contributions from the observational error due to inaccuracy in the data (**C**_{d}) and the model prediction error that is associated with approximations in modeling (**C**_{p}) (*65*, *66*). We choose both **C**_{d} and **C**_{p} to be diagonal matrices with **C**_{p} = diag (α^{2}*d*^{2}), where σ* _{d}* is the error for each component of displacement vector,

*d*, and α is a constant prefactor.

**C**

_{p}is more important than

**C**

_{d}when surface displacements are larger.

Since we do not impose explicit prior constrains on the positivity of slip along the strike direction, the posterior model distribution at each epoch most likely include models with cumulative backward slip on certain subfault patches. These individual models are less physically plausible than other models, reflecting a larger inconsistency between the data and simplified fault geometry and forward models. With the inclusion of these models, the uncertainty estimates of our models will be larger than those in inversions adopting stronger constraints (typically lacking closed-form solutions). Note that these inversions with different prior constraints will have the same best-fitting MAP (posterior mode) model, as long as our MAP model satisfies the constraint (i.e., no backward slip). However, the posterior mean model will change if model uncertainty estimates are different. Overall, we view that our approach is appropriate to derive best-fitting models unaffected by our prior assumptions and conservative for estimating model uncertainties in the sense that model errors can clearly be further reduced.

### Synthetic source inversions

We apply our inversion approach to synthetic scenarios of fault slip to quantify the resolution of our inversion and guide our choices of inversion parameters for real data (fig. S13). Each synthetic scenario features a localized 2D Gaussian strike-slip source on the SWFZ, which is parameterized as *z*_{c} and *L*_{s} are the centroid depth and spatial length scale of the source, respectively, and *S*_{max} is chosen as 0.5 m to resemble the coseismic scenario of the Parkfield event [e.g., (*16*)]. We only focus on changing the source depth, rather than lateral locations, since the resolution of source inversion varies more rapidly along depth. The surface GPS displacements that result from the given source are the synthetic data that we use in the inversion (fig. S13A). We do not include synthetic noise in these tests. For linear inversions with Gaussian assumptions on the error structure, we can interpret each inverted model as the mean of a large ensemble of inversions using synthetic data plus realizations of random noise. We perform a series of inversions using different synthetic source centroid depths *z*_{c} (1.5, 3, 4.5, 6, 7.5, and 9 km) and sizes *L _{s}* (2, 3, 4, and 5 km) and different parameters in the inversion scheme. These parameters include those for

**C**

_{m}: the spatial correlation scale

*L*(2, 3, and 4 km) and prior model uncertainty σ

_{0}(0.25, 0.5, and 0.75 m); and for

**C**

_{p}: model prediction error prefactor α (0.05 and 0.1). The inverted source models are compared with the input source models in terms of total potency, peak slip, and length scale (fig. S13, B and C).

These inversion results demonstrate depth-dependent resolution given the coverage of GPS stations (fig. S13, D to F). As expected, the inferred average fault slip (in areas with slip larger than 5% of the maximum slip amplitude) decreases quickly with depth, corresponding to the increase in resolution scales. The ability to resolve peak/average fault slip is poorer for input sources with a smaller dimension, suggesting a maximum spatial resolution of ~3 km even at shallow depths. Despite the limited resolution on local slip amplitude, the total potency is well recovered for sources at various depths.

### Inversion parameter choice

Guided by results of the synthetic tests, we adopt the following parameters in the inversions of actual data. We choose *L* = 4 km, to be consistent with the estimated maximum source resolution. We choose α = 0.1, which represents 10% of model prediction uncertainty as a means to avoid overfitting data. We choose σ_{0} = 0.5 m, comparable to the maximum amplitude of coseismic fault slip. Note that both *L* and σ_{0} control the prior model distribution, with a smaller *L* and a smaller σ_{0} generally leading to rougher slip models and smaller posterior model uncertainty, respectively. While the choice of *L* is supported by the consideration of model resolution, the choice of σ_{0} is more arbitrary and subject to how we determine our prior knowledge before analyzing the data. To model the cumulative GPS displacements at different epochs, we scale prior model uncertainty σ_{0} with the maximum amplitude of data *d*_{max} for the epoch in consideration, relative to the coseismic phase with the maximum displacement *d*_{c} at each epoch, we modify the prior constraints on an individual subfault patch *i* that is closest to each creepmeter (**m**_{0})_{i} *=**d*_{c} and (**C**_{m})_{ii} ** =** (0.1

*d*

_{c})

^{2},

### Seismicity catalog and analysis

We use the double-difference relocated catalog of seismicity from the NCSN (over 1984–2018) (*26*) and template matching–based catalog (PZ09) of early aftershocks (over 3 days after the Parkfield earthquake) (*27*). During the early postseismic period, NCSN catalog is complete for magnitudes larger than 1.5, whereas the complete magnitude of PZ09 catalog decreases from 1.5 at ~100 s to −0.5 at ~10^{5} s (*7*, *27*). For the PZ09 catalog, we find that the use of different cutoff magnitudes affects the estimated seismicity rate quantitatively while preserving qualitative features during the earliest postseismic period and having negligible effect on the moment rate (fig. S17). We, hence, use a conservative threshold for a completeness magnitude of 1.5 to calculate seismicity rate and moment rate.

We examine the characteristics of near-fault seismicity that occurs within 3 km from the fault plane within three regions of the same size on the fault plane (R1, R2, and R3 shown in Figs. 4 and 5). The along-strike spatial range is (−18.8, −8.8) km for R1 and R2 and (−29, −19.2) km for R3, with the distances relative to the hypocenter. The along-dip spatial range is (1.5, 6) km for R1 and (8.25, 12.75) km for R2 and R3. We calculate the seismicity rate as a function of time, *R*(*t*) = *N*(*t*)/∆*t*, where *N*(*t*) is the total number of aftershocks that occur over a time interval ∆*t* centered on *t*, and the moment release rate, *m _{i}* for each event

*i*and summed over a time interval ∆

*t*. We choose finer time steps than those in GPS data

*t*(

*k*) = 10

^{1.5 + k/20}s, where

*k*= 0, 1, 2, …, 140. For both

*R*(

*t*) and

*M*(

*t*), we average the rates over three time intervals to reduce noise in the estimated quantities and facilitate a comparison with afterslip behavior at similar temporal sampling, noting the caveat that the precise timing of seismicity rate could be affected by one time interval due to this averaging.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/15/eabc1606/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank J.-P. Avouac, Y. Fialko, P. Fulton, G. Barcheck, R. Lohman, and G. Abers for the discussions about this project. We thank J. Freymueller, two anonymous reviewers, and M. Ritzwoller for comments that help improve the manuscript.

**Funding:**This study is supported by NASA Earth Sciences and Interior grant NNX17AD99G and NASA MEaSUREs grant NNH17ZDA001N awarded to Y.B. J.J. was partly supported by the Green Foundation Postdoctoral Fellowship at the Scripps Institution of Oceanography.

**Author contributions:**J.J. and Y.B. designed the study. Y.B. processed GPS data. J.J. processed all datasets and performed modeling and analyses. J.J., Y.B., and E.K. all contribute to discussions. J.J. lead the writing of the manuscript with contributions from Y.B. and E.K.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials and following data repositories. Strong motion data are obtained from the Center for Engineering Strong Motion Data (https://strongmotioncenter.org/). Creepmeter data are obtained from the U.S. Geological Survey (https://earthquake.usgs.gov/monitoring/deformation/data/). The daily and high-rate GPS solutions are processed at the SOPAC (http://sopac-csrc.ucsd.edu/). The original high-rate solutions are accessible from http://garner.ucsd.edu/pub/measures_earthquake_products/Parkfield_20040928/GNSS/JiangBockKlein2021/ (username: “anonymous”; password: your email). Other datasets and models in this study are available in a public repository (https://doi.org/10.5281/zenodo.4278477).

- Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).