Scattered M3–4 Slip Bursts Within Creep Events on the San Andreas Fault

Scientists have observed the surface expression of creep events along the San Andreas Fault since the 1960s. However, the evolution of slip at depth has been examined relatively little. So here we probe that deep slip by analyzing strain observations just before and during hours‐ to day‐long creep events at the northern end of the creeping section of the San Andreas Fault. We identify 71 strain offsets that are likely produced by few‐hour bursts of slip at depth. Then, we grid search to determine the location, depth, and magnitude of these slip bursts. We find that the slip bursts occur at a range of along‐strike locations, from 0 to 7 km away from the surface slip observations. Slip occurs at depths from 0 to 10 km; 42%–55% of the bursts are likely below 4 km depth. The bursts typically have moments equivalent to Mw 3.2–4.1 earthquakes. These findings suggest that creep events are not just small shallow events; they are relatively large events that nucleate at significant depths and could play a prominent role in the slip dynamics of the creeping section.

Correlation between the timing of surface rupturing creep events, as recorded by creepmeters at different locations along the San Andreas Fault, show that creep events can span large distances, despite their small displacements (Gittins & Hawthorne, 2022).Similar observations of kilometer-long creep event ruptures have been identified using tiltmeters (Mchugh & Johnston, 1976) and strainmeters (Mortensen et al., 1977).Models to explain the physical processes controlling creep events have also invoked large along-strike lengths (King, 2019;King et al., 1973;Nason & Weertman, 1973).However, those along-strike extents are based primarily on observations of surface slip.What happens at depth remains more mysterious.For instance, we do not know whether creep events rupture broad regions of approximately 4 km along-strike, from the surface to 2 km deep, as estimated using a combination of surface strainmeters (i.e., creepmeters), tiltmeters, and water level fluctuations (e.g., Evans et al., 1981;Mortensen et al., 1977) or result from slip at a few isolated locations rupturing a few hundred meters along-strike, to a depth similar to their lengths, as estimated from surface strainmeters (i.e., creepmeters), borehole strainmeters and tiltmeters (e.g., Gladwin et al., 1994;Goulty & Gilman, 1978;Mchugh & Johnston, 1976).We do not know why some creep events are longer than others or why some events contain more than one burst of slip (e.g., Gittins & Hawthorne, 2022;Schulz, 1989) as seen in Figures 2b-2d.Are there isolated patches of the fault that particularly encourage creep events and repeatedly rupture, or are we observing the triggering of different patches of the fault that vary from event to event?
At the surface, creep events accommodate >50% of the slip along the San Andreas Fault (Gittins & Hawthorne, 2022).However, we do not know how much slip occurs at depth during creep events.The majority of estimates of slip during creep events are made using surface instruments such as creepmeters that have a limited resolution to slip at depth.Small strain offsets preceding surface creep events have been observed at San Juan Bautista using a borehole strainmeter (Gladwin et al., 1994).This observation illustrates that there are still some small slip events associated with creep events that are missing in both the observations and calculations of the amount of slip creep events can accommodate.Determining how much of this slip is missing from such a calculation requires knowledge of the depth at which slip occurs during creep events and the amount of slip.
If creep events accommodate significant slip at depth, with the moment in each event comparable to M w > 4 earthquakes, these ruptures could accommodate much of the slip budget at a range of depths.Indeed, if the events are large, bursts of creep events (e.g., Khoshmanesh & Shirzaei, 2018) might accommodate a large fraction of the moment deficit on the San Andreas Fault, which is currently estimated to be equivalent to a M w 5.2-7.2earthquake every 150 years (Maurer & Johnson, 2014;Michel et al., 2018;Ryder & Bürgmann, 2008).In this study, we use strainmeter and creepmeter observations to probe several features of creep events: their depths, along-strike locations, and moments.Specifically, we use data from two strainmeters and three creepmeters to probe creep events at the northern end of the creeping section of the San Andreas Fault.We begin in Section 2 by processing the strainmeter data.In Section 3, we identify 71 offsets in the strain time series created by short bursts of slip within longer creep events.Then, we determine the location and magnitude of these slip bursts in Section 4. We find that the slip bursts originate from a range of along-strike locations.Most bursts occur below 4 km and have moments equivalent to M w 3-4 earthquakes.We discuss the implications of these locations and sizes in Section 6.This study provides an in-depth analysis of small slip bursts that create strain offsets associated with surface rupturing creep events, which previous studies have been unable to reach or observe.

Data
We analyze creep events recorded between 2009 and 2019 on the USGS creepmeters at Harris Ranch (XHR) and Cienega Winery (CWN), located at the northern end of the creeping section of the San Andreas Fault (Figure 1).We combine the creep event observations from XHR and CWN with data from two borehole strainmeters.B065, a nearby PBO instrument, provides the primary data set.As a supplement, we use data from SJT, a USGS instrument 11 km north, near San Juan Bautista.We also use data from a third creepmeter near San Juan Bautista (XSJ) to validate our grid search results (see Section 5.1.5for details).
The creepmeters XSJ, XHR, and CWN all have the same sampling rate: one sample per 10 min.Each creepmeter has a different operating duration.The most recent iteration of XSJ operated between 2007 and 2016, while the most recent iterations of XHR and CWN have operated since 2009 and 1991, respectively.These creepmeters allow us to measure firsthand the surface rupturing creep events and are one of the few geodetic instruments that can do so.
The two strainmeters have different sampling intervals.SJT, which operated from 1983 to 2016, has a sampling frequency of 18 min.B065, operating since 2007, records 20 samples per second, but we use downsampled 10min data.Both instruments can detect sub-nanostrain variations in strain, but B065 has a slightly better resolution with fewer data gaps (Gladwin, 1984;Gladwin et al., 1987;Hodgkinson et al., 2013;Roeloffs, 2010).
Figure 2 shows data from three creep events as recorded on the creepmeters (a)-(d) and the strainmeters (e)-(h).Figures 2e-2h show strain offsets that are observed to occur before the surface recordings of creep events at XHR and CWN (See Section 3 as recorded by nearby borehole strainmeter B065, buried at ∼250 m depth at ∼5 km from the fault).Figures 2a-2d show the surface rupturing creep events that are preceded by the strain offsets shown in the associated bottom panel, with the duration of the strain offset shown as the vertical orange bar.The creepmeter and strainmeter data do not show consistent features with the strainmeters responding to slip at depth over a short period of time (∼1-2 hr), whereas the creepmeters are recording the surface response to the burst of slip at depth as a form of a relaxation process after receiving an initial kick of stress from depth which occurs over several hours to days (Bilham & Behr, 1992;Scholz, 1998;Wei et al., 2013).
Using both creepmeters and strainmeters allows us to investigate the depth, horizontal location, and magnitude of slip during creep events.

Processing and Correcting the Strainmeter Data
We begin our strainmeter analysis by calibrating the data to horizontal strain.The strainmeter time series are recorded on four horizontal gauges with four different orientations in the borehole, and we map the gauge data into the three components of horizontal strain.We use the calibration of Hodgkinson et al. (2013) for B065 and the calibration of J. Langbein and others at the USGS for the SJT.We thereby obtain three strain time series at each strainmeter: areal strain ϵ E+N , differential extension ϵ E N , and engineering shear ϵ 2EN .We discard the areal strain ϵ E+N record because it is significantly noisier.
The sub-nanostrain resolution of the borehole strainmeters means they often record "noise" from non-tectonic sources, such as borehole settling, tidal loading, and atmospheric pressure (Hodgkinson et al., 2013;Roeloffs, 2010).These signals dominate the record and make it difficult to observe the typically very small offsets (≤10 nanostrain, Figure 2) associated with creep events.We, therefore, correct for these non-tectonic signals.
For strainmeter SJT, we use the corrections obtained by Hawthorne, Simons, and Ampuero (2016).They examined the ϵ E+N , ϵ E N , and ϵ 2EN strains, discarded times with obviously large noise, and removed linear mappings to tidal loading, atmospheric pressure, and a spurious 3-hr instrument signal.
For strainmeter B065, we use slightly different components and process the data following the approach of Hawthorne, Bostock, et al. (2016).We first estimate and remove two long-term exponential trends associated with the settling of the borehole (Gladwin, 1984).
Next, we identify the strain components at B065 that respond minimally to loading by atmospheric pressure.These "non-atmospheric" components tend to have smaller hydrological noise, presumably because atmospheric pressure and hydrologic loads deform the Earth in a similar way (Hawthorne, Bostock, et al., 2016).We choose two components close in orientation to the original two shear strains.These new non-atmospheric components (denoted by the -na suffix) are a linear combination of ϵ E+N , ϵ E N , and ϵ 2EN .The values by which ϵ E+N , ϵ E N , and ϵ 2EN are multiplied are empirically calculated.These multipliers are balanced such that if we are looking at the ϵ E N-na component of strain, the majority of this component comes from the original ϵ E N component.Here, we use the following linear combinations of ϵ E+N , ϵ E N , and ϵ 2EN to calculate the non-atmospheric components, such that We construct the non-atmospheric components from the original strain data before any correction, so we must now reprocess them.As before, we remove exponential trends that represent borehole settling.Finally, we remove the last non-tectonic signal: tidal loading.We identify a set of tidal periods that are likely to have signals with amplitudes of at least 0.5 times the noise level (following estimates by Hawthorne, Bostock, et al. (2016)).
We fit a set of sinusoids with these periods to the data and then remove them.

Identifying Strain Offsets Associated With Creep Events
Now that we have processed the strainmeter data to remove non-tectonic signals, we visually inspect it.We find that there are often abrupt offsets in the strain record at B065 during and just before creep events (Figure 2).The strain offsets are typically small (a few nanostrain) and accumulate over a few hours (bottom panels of Figure 2).They are likely caused by few-hour bursts of slip somewhere on the fault.
We systematically search for strain offsets during and just before and after creep events recorded at creepmeter XHR, as identified by (Gittins & Hawthorne, 2022).For each event, we estimate a long-term trend in the strain data from the 7 days before the creep event.We remove that long-term trend from the strain time series, as it is largely non-tectonic noise.Then, we visually search for offsets in the strain data, beginning 10 hr before the surface slip onset and ending when the surface slip is finished.The time window in which we search for the strain offsets may affect their amplitude.However, when we initially considered a longer window of 48-72 hr, we observed that the strain offsets typically appeared as sharp, short-term changes in slope (overprinting a long-term trend) that were typically within 1-2 hr of the creep event onset time, with the long-term trend taking over once more after this short duration offset.Due to this, we chose to look at a shorter window of 10 hr.This time window is long enough to observe and pick the strain offsets.A 10 hr window would have sufficiently captured the majority of previously reported strain offsets associated with surface rupturing creep events (Gladwin et al., 1994).
We identify 92 potential strain offsets on the B065 strainmeter associated with creep events.However, some offsets are only marginally above the noise.We will focus on the 71 clearly identifiable offsets.Sometimes, two or more strain offsets are associated with a common creep event, so the 71 strain offsets are associated with 44 surface creep events recorded at either the XHR or CWN creepmeter.Some of these 44 creep events are recorded at both creepmeters, with approximately 20%-25% of creep events rupturing both creepmeter locations (Gittins & Hawthorne, 2022).
Figure 3 summarizes some properties of the strain offsets identified at B065.The offsets are typically ± a few nanostrain (Figures 3a and 3b).They typically last 1-3 hr (Figure 3c) and occur predominantly in the 3 hr before surface slip begins (Figure 3d).
It is interesting to note that the strain offsets vary significantly in amplitude and sign.The amplitude may be controlled by several parameters: the dimensions of the patch of fault that slipped, that patch's depth, and the slip burst's moment.However, the sign of the strain offset does not depend on the magnitude of slip; it is a function of slip location and primarily a function of the along-strike location.The variation in the sign of the strain offset is thus an immediate indication that the location of slip varies from burst to burst.
This variation in the sign of the strain offset contrasts with observations of strain associated with creep events slightly further north.Gladwin et al. (1994) inferred that a single location produced all of the strain offsets near San Juan Bautista, as the signs and relative amplitudes of the strain offsets there were constant from event to event.

Identifying the Location and Magnitude of Slip
Our next task is to interpret the strain offsets quantitatively so that we can identify the slip bursts' locations and magnitudes.To do so, we first compute a forward model; we compute the strain that would be produced by slip at a range of depths, along-strike locations, and moments (Section 4.1).Then, we assess which locations and offsets match the observations (Section 4.2).

Forward Model
We calculate the strain offsets that would be produced by 500-m-wide square patches of horizontal slip on the San Andreas Fault (We obtain very similar results with 1,000-m-wide patches).The patches are embedded within an elastic homogeneous half-space (Okada, 1985(Okada, , 1992;;Thompson, 2015).They are located at depths between 250 and 9,750 m and at along-strike locations between 13 km northwest and 17 km southeast of creepmeter XHR.In the strain calculations, we change the fault's strike according to the along-strike location, as the San Andreas Fault has a subtle change in the strike around XHR (Figure 1).We use a strike of 126°between 13 km north and 7 km south of XHR, 133°between 7 and 10.3 km south of XHR, and 137°from 10.3 km south of XHR.As the final parameter of our grid search, we vary the moment of slip.We consider a wide range of magnitudes, between M w 1-6, at 0.1 M w intervals.As one would expect, the predicted strain increases linearly with the assumed moment (Figure S1 in Supporting Information S1).
Figure 4 shows the predicted strain varies with slip location by illustrating the (a) ϵ E N na and (d) ϵ 2EN na strains produced by M w 3.6 slip bursts at each location.For visual interpretation, we have shaded the regions by the sign of the predicted strain, with the red regions being positive and the blue negative.The 5 nanostrain contours illustrate how the Green's functions evolve within each region.Sections B-B′ and C-C′ in (a) illustrate how the predicted ϵ E N na varies with depth and along strike and are shown in (b) and (c), respectively.Sections E-E′ and F-F′ in (d) illustrate how the predicted ϵ 2EN na varies with depth and along strike and are shown in (e) and (f), respectively.
Figures 4b and 4e, which represent the depth sections B-B′ and E-E′ in (a) and (d), respectively, show how the surface strain offsets vary as the slip patch is moved deeper.As the depth of the patch increases, the amplitudes of the strain offsets decrease, but at different rates on the two components.Figures 4c and 4f, which represent the horizontal sections C-C′ and F-F′ in (a) and (d), respectively, show how the strain offsets observed vary as the slip patch moves along strike.As the patch moves from NE (left) to SW (right), the amplitude and sign of the strain offsets change.
As noted above, we can use the strain offsets' signs to quickly identify regions of the fault where the slip patch may be located Figure 5.For example, the histograms in Figures 4b, 4c, 4e, and 4f illustrate estimates of the strain offset in an event on 02 February 2012 (uncertainty histogram location is shifted by the observed strain offset).The strain offset associated with this event has a negative offset on both components.Those negative offsets suggest that the slip burst should be within the dark purple regions in Figure 5.There are only two regions of the fault where both strain components are negative: just to the south of XHR or more than 7 km north of XHR.We will assess these locations more quantitatively in the next section.

Matching the Observed Strain Data
To see which locations and moments match the observed data the best, we need to estimate the strain offset in each burst and then assign each scenario in our grid search, which represents several magnitude scenarios for each spatial grid location, a probability based on how well it can reproduce the observed strain offset.
Our best estimate of the observed strain offset is simple; we subtract the strain at the end of the identified offset from the strain at the beginning, as is appropriate for offset estimates in random walk noise (Langbein, 2004;Langbein & Johnson, 1997).To calculate the uncertainties in the strain offsets, we assess how much strain changes in similar-length windows at random times.We compute the change in strain during 50,000 2-hr-long windows.The uncertainties are relatively small.Of the 50,000 random strain offsets taken at B065, 70% lie between 0.63 and 0.59 nanostrain for the ϵ E N-na component and 0.75 and 0.49 nanostrain for the ϵ 2EN-na component.At SJT, 70% of the random strain offsets lie between 2.47 and 0.67 nanostrain for the ϵ E N component and 1.88 and 1.23 nanostrain for the ϵ 2EN component.The uncertainty distribution of 50,000 2-hrlong windows can then be used to determine the probability of each grid search scenario by allowing us to estimate how well each scenario fits the observed strain offset.
To determine the probability that slip with a given location and moment produced an observed strain offset, we extract the relevant probability from the observed probability distribution (i.e., a normalized version of the uncertainties distribution estimated above).We calculate the difference between the observed and modeled strain offset to determine its horizontal position along the probability distribution.We then extract the probability of this grid search scenario, thus a probability for each combination of location and moment.For instance, the probability that M w 3.6 slip at 2 km depth below XHR produced the ϵ E N-na observed in the event in Figure 4 is proportional to the value on the histogram in panel b at ϵ E N-na = 1.5 nanostrain.We do this calculation for both strain components at B065 for all events and the two strain components at SJT when it has data available.This in essence is testing whether the difference between the observed and modeled strain offsets can be solely explained by the noise of the strainmeters.We then multiply the probabilities for each component of B065 (and SJT when available) to obtain the probability for each spatial grid given a certain moment.Finally, we sum the probability of all of the moments for each spatial grid, weighting each moment equally, before normalizing so that the probability of the whole grid search region is equal to 1. Journal of Geophysical Research: Solid Earth 10.1029/2023JB028187

Slip Burst Locations
Slip bursts appear to occur at various locations along the fault.We will discuss four groups of locations, starting with the least common and progressing to the most common.

Location 1: ~3 km North of XHR
Approximately 4% of slip bursts appear to occur approximately 3 km north of XHR.They produce strain offsets that are positive on both components.Figure 2e shows an example of one of these strain offsets that occurred on 02 November 2014 associated with a creep event at XHR (Figure 2a).
Figure 6a shows a normalized stack of the probability of all of the events that have both strain offsets that are positive on both components.The dark blue contour represents 20% of the maximum probability.For the slip burst's probability summed and shown in Figure 6a, we can see a distinctive green-yellow area of the fault north of XHR that has a higher probability than elsewhere on the fault.This area highlights the best-fitting location.

Location 2: 5 km North of XHR or 5 km South of CWN
Approximately 10% of slip bursts appear to occur either ∼5 km north of XHR or >9 km south of XHR (>5 km south of CWN).They produce strain offsets that are positive on the ϵ 2EN-na component and negative on the ϵ E N-na component.Figure 2f shows an example of one of these strain offsets that occurred on 08 November 2012, associated with a creep event at XHR (Figure 2b).
Figure 6b shows a normalized stack of the probability of all of the events that have a positive strain offset on the ϵ 2EN-na component and negative offset on the ϵ E N-na component.The dark blue contour represents 20% of the maximum probability.For the slip burst's probability summed and shown in Figure 6b, a region north of XHR and a region south of CWN appear equally likely, as indicated by the presence of yellow-green patches in both locations.Which location is preferable depends on whether the observed strain offset is associated with surface slip at either the XHR or CWN creepmeter.

Location 3: Just North of XHR
39% of slip bursts appear to be located just north of XHR.They produce strain offsets that are positive on the ϵ E N-na component and negative on the ϵ 2EN-na component.Figure 2g shows an example of one of these strain offsets that occurred on 26 June 2012, associated with a creep event at XHR (Figure 2c).
Figure 6c shows a normalized stack of the probability of all of the events that have a negative strain offset on the ϵ 2EN-na component and positive offset on the ϵ E N-na component.The dark blue contour represents 20% of the maximum probability.For the slip burst's probability summed and shown in Figure 6c, a region just north of XHR appears the most likely, with the best-fitting locations following a SE-trending yellow-green streak.

Location 4: Between XHR and CWN
The largest number of slip bursts appear to be located between XHR and CWN.They produce strain offsets that are negative on both components.Slip bursts in this region produce ≈46% of the strain offsets.Figure 2h shows an example of one of these strain offsets that occurred on 15 September 2015, associated with a creep event at CWN (Figure 2d).
Figure 6d shows a normalized stack of the probability of all of the events that have a negative strain offset on both components.The dark blue contour represents 20% of the maximum probability.For the slip burst's probability summed and shown in Figure 6d, we see two regions that appear equally likely, with a yellow-green highprobability patches of fault between XHR-CWN and beneath XSJ.However, the northern location, beneath XSJ, seems unlikely given further consideration in Section 5.1.5.It can be excluded when data from SJT is used.
One double negative strain offset is slightly different from the others, however: that on 12 July 2012 (Figure S2 in Supporting Information S1).This event has a large negative offset on the ϵ E N-na component and is located to the south of XHR, directly below CWN.The slip burst on 12 July 2012 is the only event located here.

Excluding Scenarios Under XSJ
As noted above, a number of the observed strain offsets are consistent with slip at two locations: one just north of CWN and one under XSJ (e.g., Figure 6d).We, therefore, must determine whether both are plausible.
To do this, we consider the best-fitting locations of (a) all events and (b) events that occurred when we have data from the northern strainmeter SJT; between 1983 and 2012.Figure 7b shows the five most probable slip locations for all slip bursts between 2009 and 2019.Some locations plot to the north, under XSJ.However, Figure 7a shows the five most probable slip locations for slip bursts that occurred when SJT was operational.Very few are located near XSJ.There is no reason to expect a physical change in the fault after SJT stopped, so it seems likely that all of the slip bursts occurred farther south, away from XSJ.
Northern slip bursts also seem unlikely because no significant slip occurs at creepmeter XSJ at the times of the examined creep events (e.g., Figure 2d).In contrast, when slip bursts occur farther south, near XHR for example, the XHR creepmeter records slip (e.g., Figure 2c).

The Locations of Aseismic Slip Bursts and Background Seismicity
Figure 7a shows the five best-fitting scenarios for strain offsets when SJT is operational.If we compare the location of these scenarios to the locations of background seismicity (gray dots in Figure 7a) we observe that these bestfitting scenarios appear to lie within an along-strike gap in the background seismicity that exists beneath XHR. Figure 7b also shows the relationship between the slip burst locations and the local earthquakes but the picture is not as clear.The majority of the five best-fitting scenarios occur somewhat away from the earthquakes, between patches of higher seismicity to the north and south.Some locations lie within the seismicity that exists beneath XSJ, however as discussed in Section 5.1.5,these scenarios are probably located here rather than between XHR and CWN (i.e., in the horizontal gap in seismicity) due to a lack of data available from SJT to reduce the probability of these locations.Yet it is clear from both Figures 7a and 7b that the slip bursts creating the observed strain offsets associated with surface rupturing creep events are likely occurring in the along-strike seismic gap beneath XHR.

How Many Locations Produce Slip Bursts?
We have described slip bursts in four regions based on the signs of the strain offsets they produce.There is no overlap between those regions (Figure 6), and it is clear just from the signs of the strain offsets (Figure 5) that slip bursts occur in at least four locations.However, our analysis provides no evidence that there are just four localized spots that produce slip bursts.Best-fitting locations vary, and there are significant uncertainties in each event's location; regions of high probability are smeared along several-kilometer regions in Figure 6.
If we select the five most probable scenarios/locations for each slip burst, we are able to see that instead of four distinct non-overlapping regions, we observe a diverse spread of locations.Figure 8a shows this spread in the five most probable locations of all events, as ordered by the best-fitting distance from XHR. Figure 8b shows this spread in the five most probable depths of all events, as ordered by depth.

Relationship Between the Recorded Offsets and Along-Strike Location of Slip
To determine whether the amplitude of the observed strain offsets varies with along-strike location, we compare the median value of along-strike distance of the five best fitting scenarios for each strain offset to its amplitude as shown in Figure 9a.There appears to be little to no correlation between the distance from XHR and the amplitude of the observed ϵ 2EN-na of strain.The ϵ E N-na component may be slightly higher closer to XHR, but there is a large amount of scatter in these values.This would suggest that each of the various slip locations creating the observed strain offsets can produce similarly sized strain offsets.(Waldhauser, 2009;Waldhauser & Schaff, 2008).

Slip Burst Depths
Estimates of slip burst depths are often more uncertain-than estimates of along-strike locations.This uncertainty is evident in Figure 6c where there is a smearing of the relative probability distribution with depth.
Nevertheless, we can assess whether the slip bursts are mostly shallow (0-2 km), intermediate (2-4 km), or deep (>4 km).We chose to use these boundaries for a variety of reasons.First, the shallow region (0-2 km) is chosen as it encompasses the full depth range of the previous estimates of the depth of creep events in the study region (Mortensen et al., 1977). 2 km depth also appears to represent the shallowest depth at which the majority of the seismicity south of XHR occurs (Figure S3 in Supporting Information S1).If the slip bursts predominantly lie above 2 km then we may be able to help explain this apparent lack of earthquakes above 2 km.
The intermediate region (2-4 km) is chosen as 4 km appears to be the shallowest depth at which the majority of the seismicity north of XHR occurs (Figure S3 in Supporting Information S1).If the slip bursts predominantly lie above between 2 and 4 km with a transition to shallower depths moving from north to south, then this might explain the change in the difference in the shallowest earthquake depths on either side of XHR.Between 2 and 3 km depth in the study region also represents the base of the unconsolidated sedimentary sequence deposited in the Hollister region since the Early Pliocene (Allen et al., 1946;Dibblee, 2006;Rogers, 1980).We, therefore, use this intermediate depth range to investigate if this lithological boundary could play a role.The base of this intermediate-depth at 4 km also represents the deepest previously suggested depth for creep events (Bilham et al., 2016).
Finally, we chose the deep region (>4 km) as this would represent events within the depth of seismicity to both the north and south of the XHR seismic gap (Figure S3 in Supporting Information S1).This depth is also below the lithological boundary of the base of the unconsolidated sediments in the region (Allen et al., 1946;Dibblee, 2006;Rogers, 1980).Below 4 km depth, the change in the Green's function of the grid search model at different depths becomes smaller (Figures 4b and 4e), meaning we have less resolution to resolve different depths beyond this point, increasing the uncertainty below 4 km.
We assign a given slip burst to one of these depth ranges based on whether three out of five of the best-fitting locations are within one of the depth ranges (Figure 8b).If the 5 best-fitting values are more spread out, we do not assign the burst a depth classification.The narrower widths of the shallower depth ranges could lead us to underestimate the number of shallow slip bursts.However, many of the slip bursts seem to be clearly assigned a depth range (Figure 8b).With the criteria above, we find that 13%-21% of the assigned slip bursts are shallow, 8%-17% are between 2 and 4 km, and 42%-55% are below 4 km (Figure 10).Uncertainties on the percentage are 70% uncertainties from bootstrapping the events included.We cannot assign a depth to 17%-28% of slip bursts.

Relationship Between the Recorded Offsets and the Depth of the Slip Location
Again, we seek to determine if there is a spatial relationship between the location of slip and the amplitude of the observed strain.Here we compare the median value of the depth of the five best-fitting scenarios for each strain offset to its amplitude as shown in Figure 9b.Just as with the along-strike location, there appears to be little to no correlation between the depth of the slip bursts and the amplitude of the observed ϵ 2EN-na of strain.
If we compare the along-strike location and depth of the slip bursts producing the observed strain offsets, we observe no correlation (Figure 9e).This suggests that there are no spatially distinct regions that produce the observed strain offsets but rather a range of multiple possible locations.

Relationship Between Rainfall, the Strain Offset Amplitude, and the Slip Burst Location
We compare the depths of the slip locations producing the observed strain offsets that occur within 24 hr of rainfall to those that do not.To compare whether rainfall produces a larger strain offset, we group the strain offsets by whether they occur within 24 hr of rainfall.Figure S4 in Supporting Information S1 shows the size of the strain offsets on both the (a) ϵ E N-na and (b) ϵ 2EN-na components of strain.We see no systematic changes in the size of the strain offsets with an association with rainfall.
We also investigate whether the slip bursts producing the strain offsets are systematically at shallower depths associated with rainfall.Figure S5 in Supporting Information S1 shows the slip locations producing strain offsets within 24 hr of rainfall (a) and those that are not (b).The overall slip locations appear similar, and there does not seem to be a significant shallowing of the slip locations during rainfall events.Surface recordings of creep events are known to be associated with rainfall (Roeloffs, 2001;Schulz et al., 1983).As there does not seem to be a relationship between the depth of the slip burst and rainfall, this would suggest that proposed rainfall triggering of creep events could be associated with previously unrecorded slip bursts at depth, which may be driving them as opposed to being rainfall-induced.Further work is needed to determine if this is the case as we do not know the time over which water infiltration to the depth of these slip bursts occurs.

Magnitude of the Slip Creating Creep-Event Associated Strain Offsets
The slip bursts' moments are easier to interpret.In Figure 11, we plot the five best-fitting moments for each slip burst, ordered by their best-fitting moment.Possible moments range from M w 3 to 5, and the central 70% of the moments are between M w 3.2 and 4.1.

Relationship Between the Magnitude of Slip Creating the Creep-Event Associated Strain Offsets and Its Location
In Figure 9f we can see little correlation between the size of the observed strain offset and the magnitude of the slip that created it.However, the reason for this lack of correlation likely arises from the correlation between the magnitude of slip required to produce the observed strain offset and the depth of the location where the slip occurred (Figure 9d).This correlation arises as creating a similar-sized strain offset deeper requires a larger amount of slip than would be required if the slip was shallower.
If we also compare the along-strike location of the slip bursts and their equivalent M w (Figure 9c), much like with the size of the strain offsets (Section 5.1.8)we do not observe any correlation.This suggests that the magnitude of slip required to produce the observed strain offsets is not dependent on the along-strike location but only on the depth of the slip burst.
Although the magnitude of slip for surface rupturing creep events is rarely reported, we can estimate the equivalent M w of creep events using previous estimates of the rupture extent of creep events (e.g., Evans et al., 1981;Gladwin et al., 1994;Goulty & Gilman, 1978;Mortensen et al., 1977).The resulting estimates for the moment released in creep events range from an equivalent M w of 3.0-4.2,with a median value of 3.5 (Evans et al., 1981;Gladwin et al., 1994;Goulty & Gilman, 1978;Mortensen et al., 1977).This range encompasses ≈88% of the five best fitting scenarios for all of the observed strain offsets.This suggests that the slip bursts we observe are a similar magnitude to the surface rupturing creep events which follow them.

Propagation of Slip
In our analyses, we have located bursts of slip that create offsets in the strain record.Often, several slip bursts will occur within a single creep event.By comparing the bursts' location, we can track the propagation of slip.For instance, in the November 2012 creep event, there were three sequential strain offsets (Figure 12a).Each of these three strain offsets is associated with a phase of slip at either XHR or CWN (Figure 12b).The first slip burst is located ≈2 km north of XHR (Figure 12c).It creates a strain offset with a positive ϵ E N-na offset and negative ϵ 2EN-na offset.Slip is also triggered at the surface during this slip burst and is recorded at XHR (Figure 12b).About 8 hr later, another slip burst occurs, this time just beneath XHR (Figure 12d).This burst of slip creates a negative strain offset on both components.Again, slip is triggered at the surface and is recorded at XHR (Figure 12b).Finally, a third slip burst occurs between XHR and CWN about 20 hr after the first slip burst.This burst of slip also creates negative strain offsets on both components.Slip is once again triggered at the surface but is recorded this time at CWN (Figure 12b).
The combination of migrating surface creep observations and migrating slip burst location indicates that this sequence of slip bursts is likely triggered by a southeast-propagating slip front.
We have examined other examples of creep events that encompass slip-burst sequences.We find 20 creep events that have multiple slip bursts, with ≈60% showing evidence of slip propagation.This suggests that the majority of creep events are not made up of the repeated rupture of the same fault patch but are made up of a propagating slip that triggers slip on multiple patches.

Discussion
Our analysis of the strain observations suggests that few-hour bursts of slip frequently occur during creep events, particularly at the beginning of creep events.These bursts occur at a range of locations with magnitudes between M W 3.2 and 4.1.

How Does the Size of the Slip Burst Relate to Creepmeter Observations?
We can observe strain offsets associated with slip bursts before surface recordings of creep events.Figures 13a  and 13b illustrate the relationship between the amplitude of the strain offsets produced by these slip bursts and the slip and duration of the surface rupturing creep event.Figures 13c and 13d shows how the slip and duration of surface rupturing creep events relate to the M w of the slip burst that precedes them.In all four panels of Figure 9 we do not observe any correlation between the slip and duration of surface rupturing creep events and the amplitude of the strain offsets or the M w of the slip burst.Within two-layer models of creep events (e.g., Bilham & Behr, 1992;Slater & Burford, 1979;Wei et al., 2013), slip at depth triggers slip in the upper layers of the fault, and this is recorded as a creep event on creepmeters.The lack of correlation between the equivalent M w of the slip bursts and the slip in the surface rupturing creep events recorded by creepmeters is surprising, as one would expect a bigger slip burst at depth to create a surface creep event with a larger displacement.Further work is required into the time evolution of surface rupturing creep events to better understand the physical controls on their displacement amplitudes.We need to assess whether the displacement amplitude is a function of the rheological properties of a surface carapace in which the creepmeters are installed rather than the size of the slip burst at the depth to which it is related.

How Many Locations Produce Creep Events?
Our observations suggest that, in the Harris Ranch-Cienega Winery area, slip bursts associated with creep events occur at a range of locations (Figure 8).In some creep events, multiple locations rupture sequentially (e.g., Figure 12).This suggests that several, and perhaps many, discrete fault patches can slip to produce creep events.The locations of these slip bursts appear to be within an along-strike seismic gap (Figure 7).
The observation of multiple slip locations contrasts with at least one other creep event observation: that made via strain offsets in San Juan Bautista.All of the slip bursts originate from a single location (Gladwin et al., 1994).If all of the strain offsets have the same sign combination, then they would originate from the same region of the fault shown in Figure 5.This is the case for the strain offsets observed near San Juan Bautista, which also have a similar temporal evolution (Gladwin et al., 1994).If the strain offsets are all self-similar, then it would be reasonable to suggest that they are produced by repeated slip of the same fault patch.We do not find this to be the case.Rather, slip bursts in the Harris Ranch-Cienega Winery region appear to originate from different locations, often showing along-strike propagation in slip burst sequences associated with evidence for surface rupturing creep event propagation (Figure 12).
The difference in fault behavior from San Juan Bautista and the Harris Ranch-Cienega Winery region is interesting as both regions experience multi-phase creep events and have similar displacement amplitudes (Gittins & Hawthorne, 2022).Yet, despite only being 11 km apart, these two regions display different physical models for producing the observed surface creep events, with San Juan Bautista characterized by repeated slip on the same patch (Gladwin et al., 1994) and the Harris Ranch-Cienega Winery region characterized by slip on multiple patches (this study).This may suggest that there are different numbers of patches between the two regions resulting from differing levels of shallow heterogeneity.This shallow heterogeneity may play a role in producing the slip bursts and their associated surface rupturing creep events and has been reported as having a role in producing a variety of slow fault phenomena (e.g., Ando et al., 2010;Luo & Ampuero, 2018;Skarbek et al., 2012;Wei et al., 2013).
The number of patches that slip during creep events could be important for constraining models for creep events using shallow heterogeneity.For instance, the Wei et al. (2013) model requires patches in a certain size range.It may be plausible to have one patch in that size range, but is it plausible to have many patches in the right size range: large enough to accelerate independently but too small to host earthquakes (Ando et al., 2010;Dublanchet et al., 2013;Skarbek et al., 2012).The Wei et al. (2013) model also includes no along-strike variability in heterogeneity, something that is suggested when comparing the number of fault patches in San Juan Bautista (Gladwin et al., 1994) and the Harris Ranch-Cienega Winery region (this study) 11 km away.Further targeted work is required to understand why slip bursts and creep events occur and why they have slip rates of order 10 6 ms 1 .

What Are the Implications of Patches at Depth?
Our location estimates show that many slip bursts are deep; 42%-55% occur below 4 km.This deep slip, which usually occurs at the beginning of creep events, suggests that creep events are initiated spontaneously, driven by weakening on the fault.
The deep slip burst locations appear to contradict one previously proposed model for creep events: the "kick" model.In a kick model, slip accelerates because a velocity-strengthening near-surface fault has been "kicked" by some surface perturbation, such as rainfall (e.g., Helmstetter & Shaw, 2009;Kanu & Johnson, 2011;Marone et al., 1991;Perfettini & Ampuero, 2008;Scholz, 1990;Wei et al., 2013).However, slip is only likely to accelerate significantly if (a) the perturbation is close to the slipping location and (b) the slipping region has low normal stress, so is responsive.But there is no reason to think that faults at >4 km depth are close to surface perturbations or have low normal stress.
Instead, the deep locations of our identified slip bursts, many of which occur at the beginning of creep events, suggest that slip originates in deeper patches.These deep patches may be analogous to the deeper layers of twolayer models (e.g., Bilham & Behr, 1992;Slater & Burford, 1979;Wei et al., 2013).In these models, a conditionally stable layer of the fault would slip at depth (creating the observed strain offset) and trigger slip in an upper surface layer later, which would be recorded by surface creepmeters.
On the other hand, our depth estimates suggest that slip bursts are significantly deeper than those proposed by previous two-layer models (e.g., Bilham & Behr, 1992;Wei et al., 2013) and estimates of creep events' rupture depths (e.g., Mortensen et al., 1977).It may be interesting that the bursts mostly occur near or below the 2-3 km base of the unconsolidated sedimentary sequence deposited in the Hollister region since the Early Pliocene (Allen et al., 1946;Dibblee, 2006;Rogers, 1980).However, the role of this geological boundary in slip dynamics remains unclear.

How Much Moment Do Creep Events Release?
Our analysis suggests that the median moment released in a slip burst is 3.6.If each of the 71 slip bursts has this moment, they will account for 2.2%-6.7% of the moment released by a 10 × 10 km section of the fault over a 10year period, assuming a background slip rate of 10-30 mmyr 1 .This relationship is shown in the cartoon in Figure 14.This percentage is much lower than that suggested from looking at surface recordings of creep events, where >50% of slip is accommodated in creep events (Gittins & Hawthorne, 2022).
However, the percentage could be a dramatic underestimate.We have only considered the 71 clearest slip bursts and ignored less distinct events.If we include the less distinct events, we can visually identify at least 92 bursts of slip at XHR and CWN between 2009 and 2019.Further, we may have excluded many slip bursts because we have only examined the strain records during creep events recorded at the surface.Slip bursts could occur independent of surface slip; propagating deeper slip does not always trigger slip at the surface (Evans et al., 1981;Gittins & Hawthorne, 2022).
Geodetic estimates suggest that the San Andreas Fault currently slips at a rate below the long-term geological slip rate.This low slip rate produces a moment deficit equivalent to a M w 5.2 to 7.2 earthquake every 150 years (Maurer & Johnson, 2014;Michel et al., 2018;Ryder & Bürgmann, 2008).That moment deficit could be  (Waldhauser, 2009;Waldhauser & Schaff, 2008).accommodated by earthquakes.Alternatively, the deficit could be at least partly accommodated by bursts of aseismic creep events (e.g., Khoshmanesh & Shirzaei, 2018).
Accommodating the slip deficit of a M w 5.2 earthquake would require ≈250 extra M w 3.6 slip events.For a M w 6.2 earthquake, we would require a burst of almost 800 more M w 3.6 slip events.For comparison, between 1991and 2001, Gittins and Hawthorne (2022) observed 630 creep events on 17 creepmeters along the San Andreas Fault.If the creep event rate increased by ≈85% for a 10-year period, then the slip bursts associated with them could accommodate a moment deficit equivalent to a M w 6.2 earthquake.Given that aseismic slip is known to occur in larger bursts of fault slip rate (Khoshmanesh & Shirzaei, 2018), this does not appear outside the realms of possibility.

Conclusion
We have identified and located 71 bursts of slip that produce strain offsets associated with 44 creep events at the northern end of the creeping section of the San Andreas Fault.The bursts of slip occur at a range of locations along the fault, rupturing patches separated by several kilometers.They occur at a range of depths; half are below 4 km.The identified bursts of slip have a median moment of 3.6 and accommodate at least 2%-7% of the moment budget of the fault.
The depth of these slip bursts is deeper than many previous estimates of creep event depth.Our observations suggest that creep events are not simply shallow, triggered phenomena.They originate spontaneously from a range of locations, including depths that are comparable to the depths of nearby earthquakes.Creep events could play a significant role in the slip dynamics of the creeping section.To understand the role, we will need to further explore the physical processes that allow creep events to nucleate, propagate, and interact.

Data Availability Statement
B065 is part of the Plate Boundary Observatory.This material is based on services provided by the GAGE Facility (GAGE, 2021), operated by UNAVCO, Inc., with support from the National Science Foundation, the National Aeronautics and Space Administration, and the U.S. Geological Survey under NSF Cooperative Agreement EAR-1724794.The United States Geological Survey provides creep data and the SJT strain data (USGS, 2021).Waveform data, metadata, or data products for this study were accessed through the Northern California Earthquake Data Center (NCEDC), https://doi.org/10.7932/NCEDC(NCEDC, 2014;Waldhauser, 2009;Waldhauser & Schaff, 2008).The plotted fault traces are taken from the Quaternary fault and fold database for the United States, provided by the USGS and the California Geological Survey and is available at (USGS & CGS, 2020).18.12.07.03 of the okada_wrapper used for calculating the Green's function as part of our grid search is available by MIT license and is available at https://github.com/tbenthompson/okadawrapper (Okada, 1985(Okada, , 1992;;Thompson, 2015).The StrainProc software used to clean the strainmeter data is available from (Hawthorne, 2024;Hawthorne, Bostock, et al., 2016).

Figure 1 .
Figure 1.Map of the northern end of the creeping section of the San Andreas Fault (black).Other faults in central California are shown in gray.The USGS creepmeters are colored triangles (XSJ: dark blue, XHR: medium blue, CWN: light blue), and the USGS (green) and PBO (orange) strainmeters are colored circles.Faults plotted are from the USGS/California Geological Survey Quaternary Faults and Folds database (USGS & CGS, 2020).

Figure 2 .
Figure 2. Creep and strain and creep records associated with four creep events: (a) & (e) 02 November 2014 at XHR, (b) & (f) 08 November 2012 at XHR, and (c) & (g) 26 June 2012.(d) & (h) 15 September 2015 at CWN. Panels (a) to (d) show creep records from XSJ (dark blue dashed line), XHR (medium blue solid line), and CWN (light blue dot-dashed line).The gray vertical bar represents the time shown in the corresponding panel below.The orange vertical bar denotes the duration and timing of the identified strain offsets at B065 (described in Section 3), as depicted in the associated panel below.In panels (e) to (h), the orange and brown lines represent the ϵ E N-na and ϵ 2EN-na strains at B065, respectively.

Figure 3 .
Figure 3. Histograms of the observed strain offsets recorded at PBO strainmeter B065 on the (a) ϵ E N-na , and (b) ϵ 2EN-na components, as well as (c) their duration and (d) the time between the strain offset onset recorded at B065 and the creep event onset recorded at either XHR or CWN.

Figure 4 .
Figure 4. (a) Green's functions of the ϵ E N-na component of strain.The vertical section B-B′ shows how the Green's function varies with depth beneath XHR.The B-B′ section is plotted in panel (b) along with the uncertainty distribution (black histogram), which is shifted by an observed strain offset.The horizontal section C-C′ in panel (a) shows how the Green's function varies along strike at a depth of 2.25 km.The C-C′ section is plotted in panel (c) along with the uncertainty distribution (black histogram), which is shifted by an observed strain offset.(d) Green's functions of the ϵ 2EN-na component of strain.The vertical section E-E′ shows how the Green's function varies with depth beneath XHR.The E-E′ section is plotted in panel (e) along with the uncertainty distribution (black histogram), which is shifted by an observed strain offset.The horizontal section F-F′ in panel (f) shows how the Green's function varies along strike at a depth of 2.25 km.The F-F′ section is plotted in panel (f) along with the uncertainty distribution (black histogram), which is shifted by an observed strain offset.The contours in panels (a) and (d) are in five nanostrain intervals.

Figure 5 .
Figure5.Combined Green's functions for the ϵ E N-na and ϵ 2EN-na components of strain.Four highlighted regions show the different sign combinations of the two Green's functions.The orange and brown lines are the zero contours of the ϵ E N-na and ϵ 2EN-na strain components, respectively.The sign of the regions on either side of these contours is highlighted by a + or symbol.

Figure 6 .
Figure 6.Probability plots of each fault location slipping to produce the observed strain offsets, stacked by strain offset origin location (offset sign combinations).The dark blue contours on each plot indicate the locations that are above 20% of the maximum probability, and the gray dots show seismicity from the relocated Northern California Seismic Network catalog (Waldhauser, 2009; Waldhauser & Schaff, 2008).(a) Location 1: ∼3 km north of XHR.(b) Location 2: 5 km north of XHR or 5 km south of CWN.(c) Location 3: just north of XHR.(d) Location 4: between XHR and CWN.

Figure 7 .
Figure 7. (a) Density map of the five best-fitting locations of slip bursts that occurred when data is available from SJT.(b) Density map of the five bestfitting locations for all slip bursts, regardless of SJT's operational status.The gray dots show local seismicity from the Northern California Earthquake Data Center catalog(Waldhauser, 2009;Waldhauser & Schaff, 2008).

Figure 8 .
Figure 8. Illustrations of the relative probability of the five best-fitting scenarios for each slip burst.(a) The relative probability of the five best-fitting scenarios for each slip burst ordered by distance from XHR.The dark blue dashed vertical line at 0 km represents the along-strike location of XHR, whereas the light blue dashed vertical line at 4 km represents the along-strike location of CWN.(b) The relative probability of the five best-fitting scenarios for each slip burst ordered by depth.The orange dashed vertical line represents a depth of 2 km, and the brown vertical dashed line represents a depth of 4 km.

Figure 9 .
Figure 9. (a)The median distance from XHR of the five best fitting scenarios for each strain offset plotted against the observed strain offsets on the ϵ E N-na (orange) and ϵ 2EN-na (brown) components.(b) Median depth of the five best fitting scenarios for each strain offset plotted against the observed strain offsets on the ϵ E N-na (orange) and ϵ 2EN-na (brown) components.(c) Median distance from XHR and equivalent M w of these 5 best-fitting scenarios.(d) Median depth and equivalent M w of these 5 best fitting scenarios.(e) Median distance from XHR and depth of the five best-fitting scenarios for each strain offset.(f) Median equivalent M w of the five best fitting scenarios for each strain offset plotted against the observed strain offsets on the ϵ E N-na (orange) and ϵ 2EN-na (brown) components.

Figure 10 .
Figure 10.Percentage of slip locations in each depth range.Histograms are determined by bootstrapping the events included in the percentage calculation.

Figure 11 .
Figure 11.Magnitude of the top five scenarios for each strain event.Each square is colored by its relative probability.

Figure 12 .
Figure 12.(a) Strainmeter record spanning slip bursts during the 08-09 November 2012 creep event.(b) Creepmeter record over the same interval, showing the creepmeter records of XSJ (dashed dark blue), XHR (solid medium blue), and CWN (dot-dashed light blue).In both (a) and (b), the durations of the three strain offsets are shown as dark, medium, and light green vertical bands.(c), (d), & (e) Illustrations of the probability of slip at different fault locations producing the observed strain offsets.

Figure 13 .
Figure 13.(a) Slip recorded on creepmeters in the surface rupturing creep event following the strain offset plotted against the amplitude of the strain offset.(b) Duration of the surface rupturing creep event following the strain offset plotted against the amplitude of the strain offset.The orange and brown dots in panels (a) & (b) represent the offsets recorded on the ϵ E N-na and ϵ 2EN-na components of strain, respectively.(c) Slip recorded on creepmeters in the surface rupturing creep event following the strain offset plotted again the equivalent M w of the slip bursts that create the observed strain offset.(d) Duration of the surface rupturing creep event following the strain offset plotted against the equivalent M w of the slip bursts that create the observed strain offset.

Figure 14 .
Figure 14.Schematic diagram showing the moment release of a 10 × 10 km region of the fault below XHR (blue region) over a 10-year period compared to the moment released by 71 slip bursts (orange region) which have an equivalent M w of 3.6.The gray dots show local seismicity from the Northern California Earthquake Data Center catalog(Waldhauser, 2009;Waldhauser & Schaff, 2008).