Monitoring Velocity Change Over 20 Years at Parkfield

We monitored the time history of the velocity change (dv/v) from 2002 to 2022 to investigate temporal changes in the physical state near the Parkfield Region of the San Andreas Fault throughout the interseismic period. Following the coseismic decrease in dv/v caused by the 2003 San Simeon (SS) and the 2004 Parkfield earthquakes, the dv/v heals logarithmically and shows a net long‐term increase in which the current dv/v level is equivalent to, or exceeding, the value before the 2003 SS earthquake. We investigated this long‐term trend by fitting the model accounting for the environmental and coseismic effects to the channel‐weighted dv/v time series. We confirmed with the metrics of Akaike information criterion and Bayesian information criterion that the additional term of either a linear trend term, or a residual healing term for the case where the healing had not been completed before the SS earthquake occurred, robustly improved the fit to the data. We eventually evaluated the sensitivity of the dv/v time history to the GNSS‐derived strain field around the fault. The cumulative dilatational strain spatially averaged around the seismic stations shows a slight extension, which is opposite to what would be expected for an increase in dv/v. However, the cumulative rotated axial strain shows compression in a range near the maximum contractional horizontal strain (azimuth of N35°W to N45°E), suggesting that the closing of pre‐existing microcracks aligned perpendicular to the axial contractional strains would be a candidate to cause the long‐term increase observed in the multiple station pairs.


Introduction
Continuous seismic monitoring of the Earth on decadal timescales allows for the exploration of environmental, volcanic, and tectonic phenomena that control the seismic properties of the subsurface.Observed changes in the seismic velocity of the subsurface often result from a combination of different effects, and decadal surveys are required to untangle the contribution of each mechanism.For example, the recovery from the 1999 Chichi Earthquake is masked by environmental fluctuations such as rainfall (Feng et al., 2021), and hydrological loads overprint volcanic activities in Mount St. Helens (Hotovec-Ellis et al., 2015).
The seismic quantity dv/v provides a relative measure of the volume-averaged perturbation in the seismic velocity of the subsurface and is usually measured using the coda of waves from repeated sources and receivers.Coda waves are often used since they are sensitive to small perturbations in the subsurface (Lobkis & Weaver, 2003;Snieder et al., 2002) and ballistic wave measurements are complicated by varying source locations or ambient field stationarity (Colombi et al., 2014;Takano, Brenguier, et al., 2019).Reliable repeated coda sources and receivers may come from either repeating earthquakes recorded at the same stations (Poupinet et al., 1984) or from repeated ambient noise cross-correlation functions (CCFs) (Sens-Schönfelder & Wegler, 2006).Poupinet et al. (1984) first applied the doublet method using two microearthquakes occurring before and after the 1979 M5.9 Coyote Lake earthquake and found the time delay in the coda of the doublets increased linearly with time, indicating the change in dv/v over the medium.Numerous studies have since shown that the seismic properties of fault zones can be monitored by dv/v using waveforms from any earthquake doublet, or set of repeating earthquakes, with similar enough waveforms to compare the arrival times (e.g., Hotovec-Ellis et al., 2022;Rubinstein et al., 2007;Sawazaki et al., 2015;Schaff & Beroza, 2004;Sheng et al., 2021).
One of the most popular approaches used to extract the continuous measurement of dv/v in time uses the crosscorrelation of ambient seismic noise (Nakata et al., 2019).The advantage of using ambient seismic noise is that one does not have to wait for earthquakes to occur, which enables continuous measurements of dv/v.The dv/v is sensitive to the groundwater level (GWL) (Clements & Denolle, 2018, 2023;Illien et al., 2021Illien et al., , 2022;;Mao et al., 2022;Nishida et al., 2020;Rivet et al., 2015;Sens-Schönfelder & Wegler, 2006) and the thermoelastic deformation induced by the air temperature (Colombero et al., 2018;Gassenmeier et al., 2016;Richter et al., 2014), or even both (Feng et al., 2021;Lecocq et al., 2017;Wang et al., 2017).When interpreting dv/v measurements, it is thus crucial to identify the environmental factors as well as other mechanisms that cause variations in dv/v, such as a drop in velocity associated with coseismic earthquake rupture, logarithmic healing following an earthquake, and the other potential tectonic factors (Clements & Denolle, 2023;Feng et al., 2021;Taira et al., 2015).
Regardless of the type of wavefield used (e.g., earthquakes or ambient field), dv/v measurements vary in time due to the changing tectonic strains associated with the earthquake cycle.Before the earthquake, laboratory studies predict a slight decrease in the seismic velocities (Shreedharan et al., 2021), which has been observed only in Parkfield, CA (Niu et al., 2008) and Italy (Chiarabba et al., 2020).During the earthquake, many studies have investigated the effect of coseismic damage on the seismic velocities (e.g., Wegler & Sens-Schönfelder, 2007), either near the fault due to the extremely high strain rates (Brenguier, Campillo, et al., 2008;Lu & Ben-Zion, 2021;Nimiya et al., 2017;Taira et al., 2015;Wu et al., 2016) or at the surface as a result of the strong ground motions (Bonilla et al., 2019;Gassenmeier et al., 2016;Viens et al., 2018).Studies have found significant correlations between reductions in the seismic velocities with peak dynamic strains or stresses (Hobiger et al., 2016;Richter et al., 2014;Taira et al., 2015;Viens et al., 2018).After the earthquake, the seismic velocity appears to heal logarithmically, comparable to the recovery governed by the slow dynamics (Johnson & Sutin, 2005;TenCate & Shankland, 1996;TenCate et al., 2000).Theoretical arguments suggest that healing should saturate after some time (Snieder et al., 2017).However, the saturation time scale is poorly understood and observed healing times vary tremendously (Clements & Denolle, 2023;Illien et al., 2022;Viens et al., 2018).The postseismic healing is observed to last days, weeks, or months (Clements & Denolle, 2023;Marc et al., 2021;Viens et al., 2018), whereas little attention has been paid to the inter-seismic period between earthquakes.This is partly because inter-seismic loading usually occurs at low strain rates relative to the co-and post-seismic rates (Shreedharan et al., 2021) and is very difficult to identify due to contamination of the signal by various environmental processes and other seismic events.
This study focuses on a particularly well-instrumented segment of the San Andreas Fault (SAF) near Parkfield, California.These instruments were installed by the Berkeley Seismology Lab and have provided continuous seismic data since 2002 (https://doi.org/10.7932/HRSN).The seismic stations are located in the ruptured area of the 2004 Parkfield earthquake.Brenguier, Campillo, et al. (2008) measured the dv/v at Parkfield from 2002 until 2008 and showed a coseismic decrease in velocity associated with two earthquakes and their logarithmic healing, which they interpreted as a postseismic response.Subsequently, Wu et al. (2016) extended the analysis period to 2011 and evaluated the dv/v with different frequency bands to investigate depth-dependent velocity perturbations.They showed the coseismic decrease in velocity is larger with a higher frequency band corresponding to the shallow depth, which could be partially explained by the loss of the dv/v sensitivity for the layered perturbation at depth with the low-frequency range (Obermann, Planès, Larose, Sens-Schönfelder, & Campillo, 2013;Yuan et al., 2021) and potentially depth-dependent damage.examine models that include the effects of precipitation, temperature, earthquake healing (e.g., SS and PF), and the long-term trend term.This latter term is included as either a linear term or a residual healing term in our analysis to quantitatively determine the non-negligible increase in dv/v, which is observed by the greater value of dv/v potentially in 2022 relative to the value before the 2003 SS earthquake.
Non-linear elasticity rheology predicts that dv/v is proportional to dilatational strains (see review in Clements & Denolle, 2023).Many observations of such relation support the theory: tidal oscillations and the volcanic activities (e.g., Donaldson et al., 2017;Hirose et al., 2017;Hotovec-Ellis et al., 2022;Mao et al., 2019;Nishida et al., 2020;Sens-Schonfelder & Eulenfeld, 2019;Takano et al., 2017;Takano, Nishimura, et al., 2019;Ueno et al., 2012;Yamamura et al., 2003).In this study, we evaluated the surface strain field at Parkfield using the Global Navigation Satellite System (GNSS) data set.We compare the GPS-derived strains with the dv/v measurements to constrain the effects of tectonic deformations on dv/v measurements.The strain field around the SAF is accumulated (e.g., Agnew, 2014) and perturbed by creeping on the shallow fault (Bacques et al., 2018) and by the intricate pattern of locked and creeping patches deeper on the fault surface (Jolivet et al., 2015).In order to compare this spatially complex strain field with dv/v measurements, we calculate the spatially averaged cumulative strains in the region of the SAF that includes the borehole seismic stations.
In this study, we describe the methods and processing workflows used to obtain stable dv/v measurements at decadal timescales (Section 2), which required developing new high-performance software.We then present a survey of the dv/v results, with a particular focus on several distinct station-component pair combinations and frequency bands (Section 3).To fit models and quantify their sensitivity to the dv/v time series, we employ Markov chain Monte Carlo (MCMC) for Bayesian inference (Section 4).We then compare the strain field at Parkfield using the GNSS data to the dv/v time series (Section 5).The quantitative contributions of different physical mechanisms of both tectonic and environmental origins to the observed dv/v time series are constrained and discussed (Section 6).The abbreviations used in this article are listed in Table 1.

Methodology
Day-long instrument corrected time series of continuous (20 Hz sampling rate) three-component seismic data from January 2002 to May 2022 from 13 stations of the High-Resolution Seismic Network (HRSN; Figure 1) were used to compute the dv/v time-series.We removed the instrumental response at the beginning of data processing.Each of the three components, rather than just the typically used vertical component, was used in this analysis in order to improve the dv/v quality (Schaff, 2012) and the temporal continuity after channel-weighting (see Section 2.6).
The calculation of a dv/v time series from an ambient seismic noise data set of this size (1.6 TB) is a computational challenge due to the extremely large number of cross-correlations that are required to be computed.To overcome this computationally prohibitive obstacle, which prevents the use of existing seismic processing tools, we developed a Julia-based software SeisMonitoring.jl(Okubo, 2023c) which can be used in highperformance computing environments.Julia is a powerful language to perform computationally intensive processes (Bezanson et al., 2017) such as computing large numbers of cross-correlations in parallel.Jones et al. (2020) developed the Julia-based software tool, SeisIO.jl, to handle the seismic waveforms with metadata in the structure and showed the advantage in using memory and computational speed.Clements and Denolle (2020) developed the Julia modules called SeisNoise.jl,which efficiently computes the crosscorrelations of short-time windows to conduct the dv/v monitoring.The new software tool presented here, SeisMonitoring.jl,uses these packages as dependencies, which greatly reduces the processing time with process parallelization.

Data Quality
We first evaluated the power spectral density (PSD) of the seismic waveforms to investigate the data quality regarding the data continuity and noise levels.We used the Blackman-Tukey method to estimate the PSD, which is based on Wiener-Khinchine's theorem, such that the PSD can be obtained with the Fourier transform of the auto-correlation functions (ACFs).This scheme is efficient for our case compared to Welch's method as we can reuse the monthly stacked ACFs to compute the decadal history of the PSD.We computed the daily ACFs after removing the instrumental response associated with each station and channel.We then applied the Tukey window with 5% on the ACFs and computed the FFT to obtain the discrete Fourier coefficients, which represent the PSD of the day.
Figure S2 shows the PSD associated with the station EADB as a representative station pair.The PSD of the other stations can be found in Data Set S1.Some defects in the waveform data, which did not pass the quality control thresholds, were removed during the noise preprocessing.Sudden changes in the power spectrum amplitude for several stations and times are observed in the time series (e.g., EADB-1 and EADB-2 after 2014; Figure S2 in Supporting Information S1).They are likely due to instrumental changes and artifacts, which are not noted in the metadata.These step-like changes can potentially cause artifacts in the estimation of dv/v, which is discussed in the next section.Notably, there was a prominent peak in the PSD at 1 Hz in the channel EADB-3 around 2012-  (Johanson & Bürgmann, 2010).The stars indicate the hypocenters of major earthquakes; the SS earthquake is obtained from McLaren et al. (2008), and the others from SRCMOD (Mai & Thingbaijam, 2014).The seismic station locations with their name can be found in Figure S1 in Supporting Information S1.SN: South Napa earthquake, RC: Ridgecrest earthquake sequence.
2014, which could also be due to instrumental and/or environmental noise.This type of noise is observed in some other stations, such as FROB, contaminating the ACFs as the coherence of the periodic noise becomes dominant.As explained later in this section, we mitigate this artifact in the noise cross-correlations during the stacking process.
We cross-verified the data quality estimated from the PSDs using the supporting information of Shelly (2017).In this study, we did not correct the swap of channels or the reversed polarity on the raw data reported by Shelly (2017).Instead, we applied multiple thresholds during the processing of the seismic noise associated with the amplitude of hourly and daily CCFs and the Pearson's correlation coefficient (CC) between the reference and current CCFs during the stacking so that the perturbed CFs are mostly excluded from the stacked CCF to minimize the bias caused due to those instrumental issues in the analysis of dv/v.

Removal of Coherent Transient Signals
The  (Allen, 1978) to extract impulsive earthquake and emergent tectonic tremor events.
The kurtosis within the moving window is sensitive to the spike-like transient signals that detect the events.We followed Baillard et al. (2013) and used kurtosis as the characteristic function to pick the events.We calculated the kurtosis for every three-minute-long segment.An event is detected when the kurtosis exceeds a given threshold.The standard STA/LTA algorithm can detect tremor-like events.Unlike for detecting earthquake signals, we chose time window lengths of STA/LTA tuned to detect relatively long-term perturbations such as tectonic tremors: 3 minutes and one day for the STA and LTA lengths, respectively.The threshold of kurtosis and STA/LTA is set to be three.The combination of kurtosis and STA/LTA detection is thus performed to clean the raw data containing various types of transient signals.
We removed the signal in the time windows when either kurtosis or STA/LTA exceeded the threshold by applying the inverted Tukey window.We computed the percentage of taper with respect to the removed time windows such that the taper duration is fixed at 30 seconds. Figure S3 in Supporting Information S1 shows an example of a waveform observed at EADB and the removal of transient signals.Figure S4 shows the data availability before and after removing the transient signals.We computed the fraction of missing data in the daily waveforms associated with each station and channel.A sufficient amount of continuous data is available even after removing the transient signals from the raw data.

Auto-and Cross-Correlation Functions
We computed the CFs for all the combinations of station and component in Julia using SeisMonitoring.jl.We sliced the data into hourly time windows with half an hour of overlap.We applied detrending, demeaning, and tapering on the time windows before computing the FFT and the CFs.The maximum lag time of the correlation is one hundred seconds, which is selected to be long enough to allow for the application of a taper to mitigate edge artifacts from the band-pass filtering.The resolution of the CF can be improved by either spectral normalization (Bensen et al., 2007;Viens et al., 2017) or temporal normalization methods such as the one-bit filter (Bensen et al., 2007;Campillo & Paul, 2003;Durand et al., 2011;Larose et al., 2004;Seydoux et al., 2016;Shapiro & Campillo, 2004;Shapiro et al., 2005), or using non-linear filters (Baig et al., 2009;Hadziioannou et al., 2011;Moreau et al., 2017;Viens & Van Houtte, 2019).In this study, however, we computed the CFs without those normalizations for simplicity in the processing and interoperability of the correlated wavefield.The transient removal would replace the temporal or spectral normalization, although its efficiency remains to be evaluated in this study.
To analyze the dv/v in different frequency bands, we applied the continuous wavelet transform (CWT) to the CFs as accurate narrow-bandpass filters.The theoretical background of CWT is summarized in Torrence and Compo (1998).Mao et al. (2020) showed the wavelet cross-spectrum approach improves the measurement of dv/ v compared to the standard doublet method.Yuan et al. (2021) showed the combination of CWT and stretching method efficiently retrieves the dv/v with various frequency bands.We thus applied the CWT and inverse CWT to reconstruct the CFs in the time domain with the frequency bands of 0.2-0.5, 0.5-0.9,0.9-1.2, and 1.2-2.0Hz.We used the Morlet mother wavelet for the CWT, and the factors used to reconstruct the waveforms are taken from Table 2 of Torrence and Compo (1998).The module is implemented in SeisDvv.jl in Julia, which is translated from the python module pycwt.We applied the tapering on the CFs to avoid the artifacts at the edges of CCFs after filtering with the CWT.Note that the spectral normalization is not applied in this study.
After filtering the CFs into frequency bands, we threshold out bad CFs using the maximum amplitude of hourly CFs before computing the daily stacks.Most of the transient signals have been removed from the waveform during the pre-processing.Still, we continue curating the result by removing the remaining cause of the perturbation of the stacked CFs.We detect and mask the hourly CF stacks with low coherency, which is less likely to improve the stationary phases in the day-stacked CFs.To detect them, we first evaluated the maximum amplitude of CFs associated with each hourly time window and subsequently computed the median of these values over a day to detect the outliers in the set of hourly CFs.We rejected the hourly CFs if the maximum amplitude of CFs is greater than three times or smaller than 10% of the median.This median-mute filter allows for cleaning the daily stacked CFs in terms of the fluctuation in the amplitude caused by the remaining transient signals or the instrumental issues.Note that this filter cannot address the incoherent phases with an amplitude similar to the ambient noise due to the small earthquakes.The artifacts of the incoherent CFs are thresholded out during the stacking phase implemented as the selective stack explained in Section 2.4.
Figure S5 in Supporting Information S1 shows the nine components of CCFs associated with the station pair of LCCB-SCYB for the frequency range of 0.9-1.2Hz.It shows the relatively strong coherence in the ballistic wave arrivals followed by the coda waves.The other station pairs can be found in Data Set S2.The original data of CCFs are available in the cloud storage DASway (Okubo, 2023a).Note that the location codes and channel names have been replaced between November 2010 and September 2011.We thus unified them to obtain the continuous set of CFs, given that they are recorded on the same site.

Stacking
We computed the reference CFs by stacking the daily CFs from January 2010 to May 2022.We assume the reference CFs converge enough to evaluate the dv/v with this period.The bottom waveforms of Figure S5 in Supporting Information S1 show the reference stacks for the cases between 2010 and 2022 and the stack over the entire study period as a comparison.We applied the median-mute filter using the maximum amplitude of daily CFs to the reference stacks.The comparison of reference stacks over the different periods shows nearly identical waveforms, indicating this reference is converged enough to measure the dv/v to the current monthly stack.The distinct improvement in the stacked CFs with the median-mute can be found in other station pairs (e.g., EADB-EADB in Data Set S2).
To evaluate the time history of dv/v, we stacked the CFs over 30 days with a sliding step of 15 days used as the current CFs.Numerous stacking methods have been proposed to improve the coherent signals in the stacked CFs (e.g., Kanasewich et al., 1973;Korenaga, 2013;Nakata et al., 2015;Pavlis & Vernon, 2010;Schimmel & Paulssen, 1997;Ventosa et al., 2017).In this study, we selected the daily time windows to be stacked with the high correlations between the 12-year reference stack and current CFs to exclude further the time windows that degrade the convergence of the stacked CFs.This selection approach to enhance the signal-to-noise ratio (S/N) of CFs has been proposed in previous studies (e.g., G. Liu et al., 2009;Olivier, Brenguier, Campillo, Lynch, & Roux, 2015;Thangraj & Pulliam, 2021).The metric of selection to improve the stacked CFs varies concerning the purpose of stacking.G. Liu et al. (2009) showed the weighted stack of the common midpoint gathers with the local correlation, and the weights were defined as the CC of short-time moving windows.Olivier, Brenguier, Campillo, Lynch, and Roux (2015) used the S/N associated with the S-wave window to enhance the coherence in the stacked CCFs.These metrics aim to improve the S/N of the coherent signals to identify the reflections or the wave arrival times.
In our case, we need to enhance the stability of the stationary phases in the coda part of CFs.Thus, we used Pearson's CC to evaluate the similarity between the reference stack and the daily CFs.We computed the monthly stacked CFs using the daily CFs with which the CC to the reference is greater than the threshold.We set the threshold to zero with a range of [ 1, 1], excluding only the correlations that have shifted as much as it would degrade the stacked CFs.This metric is similar to the global correlation defined in G. Liu et al. (2009).The CC is evaluated over the entire lag time of CFs, though the similarity of high amplitude signals dominates it.It should be mentioned that the threshold based on the CC does not need to be strict because the target phase shifts that are induced by velocity changes decrease the CC.Otherwise, it could cause the underestimation of dv/v as the CC decreases with the phase change.The purpose of this threshold is to exclude the CFs that show a large discrepancy to the reference stack due to the perturbation of the noise sources so that the measurement of dv/v is assumed to be implausible.The comparison of the performance for the selective stacking to the other stacking methods can be found in Yang et al. (2022).In our case, 70% of monthly stacks are not affected by the selective stack, that is, the linear stack is applied.However, only considering the rest of the monthly stacks (the ones that had lower acceptance ratios), an average of around 40% of the daily CFs were not selected, suggesting the selection is needed for the robust measurement of dv/v.Although the selective stack performs similarly to the linear and robust stacks (Pavlis & Vernon, 2010) in their data set, the selective stacking helps stabilize the current CFs for our data set as the transient perturbations due to the small events and the environmental noises might remain in the waveform.

Measurement of dv/v
We estimated the dv/v by the phase difference between the reference and current stacked CFs using two different schemes: the stretching method (Lobkis & Weaver, 2003) and the moving window cross-spectral (MWCS) analysis (Clarke et al., 2011;Poupinet et al., 1984).Using both ways demonstrates the difference in the stability of dv/v measurements (Hadziioannou et al., 2009) and the magnitude of its estimation such that the range of estimated dv/v with the MWCS is more likely to be lower due to the weighted linear regression with the local coherency in the short-time moving window (e.g., Hillers et al., 2019).
For the stretching method, we computed the stretch coefficient ɛ, equivalent to the dv/v, in two steps following Viens et al. (2018).The coda window where we measure the phase change to evaluate dv/v is selected as the following: The start time of the coda window is either set as twice the arrival time of the ballistic wave estimated by the distance of the station pair divided by the wave velocity, or the minimum threshold of 5 s if the distance of the stations is close, and for the cases with the ACFs.We assumed the wave velocity to be 1 km/s.The end time of the coda window is fixed at 40 s.We applied the threshold associated with the length of the coda window such that it is more than five wave periods related to the central frequency of the frequency band.If the station pair did not meet the threshold, we excluded the station pair for the evaluation of dv/v.
We conducted a grid search to find the best stretching coefficient with a spacing of 0.02% to compute the profile of the CC as a function of the ɛ between the reference and the dilated current CF.We used both positive and negative lags to evaluate the CC.We then applied the spline interpolation of the profile to obtain a finer coefficient estimation following Viens et al. (2018).Note that the ACFs should have the identical measure of dv/v associated with a pair of cross-components (e.g., BP1-BP2 and BP2-BP1) due to the symmetry of CFs.However, we obtained slightly different estimations in our analysis caused by the subtle difference in the spline interpolation.We quantified all the component pairs and confirmed that the difference is mainly within a single spacing of the grid search so that it is not critical to the statistics of the dv/v time history.
For the MWCS method, we followed the process flow described in Clarke et al. (2011).The short moving time window is set as 6 s with a step of 3 s (i.e., 50% overlap).Like the stretching method, we used both the positive and negative time lags to estimate the best dt/t, that is, the phase shift over lag time.Note that we used the same criteria for the selection of the coda window with the stretching.The dt/t is the slope of the linear trend of phase shift against lag time, which is obtained using the weighted linear regression along the lag time either with or without imposing the intersection at the origin of the lag time.The latter excludes the offset of dt, which helps mitigate the artifacts of instrumental clock drift in the case of the CCFs.We discuss it more in Section 6.1.2.We used the results with the later scheme to mitigate the artifacts of the clock drift.
We evaluated the quality of the dv/v estimates with the correlation coefficients between the reference and the stretched CFs with the best-fit dv/v value and the estimation of the error derived in Clarke et al. (2011) for the cases with the stretching method and MWCS analysis, respectively.Too large differences between the reference and current CFs may show a much greater magnitude of the dv/v.Still, they may be caused by the source perturbation or instrumental issues rather than the velocity change of the structure.Given the published analyses on Parkfield (Brenguier, Campillo, et al., 2008;Delorey et al., 2021;Wu et al., 2016), we do not anticipate large changes in dv/v.We carefully selected the threshold as it should be soft enough not to cause bias due to removing the outliers.We set the threshold of 0.7% and 0.02% associated with the CC after stretching and the measurement error in MWCS, respectively, for the following analysis of the dv/v.
Figure 2 shows the number of station-component pair combinations after the thresholding of dv/v.The hundreds of pairs, including auto-and cross-component correlations with the ACFs and CCFs, are obtained to evaluate the time history of dv/v.The number of available pairs decreased due to the decommissioning of the RMNB after 2011.The area highlighted in gray indicates the period where the dv/v is scattered due to the clock drift on EADB and GHIB, discussed in Section 6.1.2.Overall, we used around 380 station and component pair combinations to evaluate the time history of dv/v.The datasheets of dv/v and the measurement error associated with the stretching and MWCS methods can be found in Data Sets S3 and S4, respectively.

Channel Weighting
To conduct the model fitting, we computed the channel-weighted time history of dv/v associated with all the possible auto-and cross-correlation station pairs.The time history of dv/v needs to be continuous and long enough to fit the models.We thus computed the weighted average of the dv/v associated with the nine components channel correlations following Hobiger et al. (2014), such as where dvv(t) is the weighted average of dv/v, and c k (t) is the CC after dilating the current CFs with the estimated dv/v.The subscript k indicates the component of the channel pairs.N is the number of channels to be averaged.We include the channels only if the quality of the estimated dv/v exceeds the threshold.The weighted average of CC is written as . (2) We also evaluated the error of the weighted average of the dv/v with the propagation of error as follows: where σ k is the error of kth dv/v.
For the case with the stretching method, the CC is computed with the measurement of dv/v.We used the error estimation obtained by Weaver et al. (2011) as σ k .We used the error derived in Clarke et al. (2011) for the case with the MWCS and simplified the averaging of dv/v by setting the CC to be unity for all the channel pairs.
We computed the fraction of the valid dv/v over the entire period and selected the station pairs with more than 0.7 of the fraction.We removed 33 and 27 pairs from 83 station pairs with this threshold for the cases with stretching and MWCS, respectively.Note that 8 of the 91 possible station pairs were already excluded during the preprocessing steps due to the large interstation distances and the low quality of the data.The 83 station pairs are listed in Table S1 in Supporting Information S1.The measurement of dv/v is sometimes unstable due to, for example, cycle skipping (Mikesell et al., 2015), which would cause bias in the model fitting.We thus applied the threshold on the dv/v measurement to ignore absolute variations with amplitude greater than 0.3%, which are assumed to be uncorrelated to the velocity change of the structures, to remove those outliers.The channelweighted average of dv/v for all the station pairs can be found in Figures S7 and S8 in Supporting Information S1, used for the model fitting, as described in Section 4.

Computational Effort
The total volume of the original waveform data is 1.6 TB.We used a Linux workstation with 48 cores, FASRC Cannon compute cluster (https://www.rc.fas.harvard.edu),and Frontera at the Texas Advanced Computing Center (TACC) (Stanzione et al., 2020) to develop the software tools and to conduct the case studies.We parallelized the processes using the Julia-native function Distributed.pmap to distribute the tasks to the workers.The heaviest computational process is to compute the cross-correlation due to the frequency of file I/O to access the waveforms and the number of station-component pair combinations.We optimized the parallelization by first computing the FFTs associated with a given time chunk for all the stations and channels in a node.We then distributed the FFTs from the master process to the workers in parallel to obtain the CFs.We submitted the jobs for the cases, for example, every 2 years, and iterated the time chunks to complete the jobs so that we obtained the CFs for 20 years with reasonable computational time.The other processes, such as downloading data, removing the transient signals, stacking, and evaluating the dv/v are also parallelized using the pmap.Further discussion can be found in Text S1 in Supporting Information S1.
or error as shown in Figure 3.The color contour shows the number of pairs within a bin of 0.02% dv/v as a proxy of the probability distribution with respect to the time bins.We plotted the median and the first and third quantiles of the dv/v.The temporal resolution is about a month, as the current CFs consist of the 30day stack with a step of 15 days.The shaded area indicates the period, which could contain the clock drift (see Section 6.1.2).The reference period is from January 2010 to May 2022, as annotated in the bottom arrows.
The coseismic velocity decreases, and the subsequent recovery phase, as shown in Brenguier, Campillo, et al. (2008) and Wu et al. (2016), are reproduced in this study up to 2011.The dv/v continues to recover and exceeds the dv/v level prior to the SS earthquake.This long-term increase is of great interest in this study.The fluctuation of dv/v measurement and the absolute change of dv/v are smaller with the MWCS, similar to the comparison of Hillers et al. (2019).
Data consistency, as measured by the progressive decrease in the range of the quantiles, increases over time.Therefore, the interpretation of physical phenomena acting before the earthquakes (SS and PF) is not likely to be robust, given the lack of stability in the dv/v measurements.In contrast, a robust feature is that the current dv/v would eventually be much greater than before the SS earthquake if the dv/v keeps increasing for more decades without the two coseismic drops in dv/v.A relatively steep positive change in the dv/v for the case using the stretching method was found around 2014.This could be most likely the artifacts due to the sudden change in the noise power spectrum (see Section 2.1).
Given that the rapid increase does not appear in the MWCS measurements, we interpret that there must be a specific scattered wave of higher amplitude that dominated the stretching measurements but that was downweighted in the MWCS calculation.While it appears in 2014, it is not strictly aligned with the date of the Napa earthquake.We explored its relation with the reports of non-stationary rates of Low-Frequency Earthquake occurrence (Delbridge et al., 2020), which we further discuss in Text S2 in Supporting Information S1.A few station pairs do not include the transition in their channel-weighted dv/v time histories, as shown in Figures S7 and S8 in Supporting Information S1.We further discuss the artifacts due to the non-tectonic origins in Section 6.1.We note that correcting this step-like change after 2014 does not remove the long-term trend in dv/v, which exceeds the original baseline level prior to the occurrence of either of the SAF earthquakes (e.g., SS).
We analyzed the dv/v time history with different sets of station-component pair combinations.Figure 4 shows the dv/v for the cases with the single-station and cross-stations (i.e., ACFs and CCFs) with the nine-channel correlations and the sets of channel correlations for both ACFs and CCFs associated with the vertical and horizontal components.The long-term trend is nearly identical between the different groups of pairs, except the vertical component of the stretching method shows the relatively quick healing of the dv/v.
Figure 5 shows the dv/v with the different frequency ranges.We used the same length of the coda window and the thresholds in the CC or error with them.We do not show the results associated with the lowest frequency band of 0.2-0.5 Hz from the case with the stretching, as the estimation of dv/v was unstable.The coseismic drop in velocity increases with frequency, as seen by Wu et al. (2016).We also observed the seasonal perturbations, which need to be excluded to evaluate the dv/v healing.The MWCS measurements in the frequency band 0.2-0.5 Hz show neither the coseismic velocity decrease nor the continuous healing of dv/v as they were likely to be masked by the large variation in the dv/v measurement.The quality of dv/v measurement at low-frequency ranges might be improved by optimizing the processing parameters, for example, the stacking period, coda window length, and the short-time window length associated with the MWCS.Wu et al. (2016) first reported the increase in dv/v drop with seismic frequency in Parkfield and interpreted it as a depth-dependent damaged structure.However, it also comes as a natural decay of the sensitivities for the layered perturbation of the medium with frequencies (Obermann et al., 2015;Obermann, Planès, Larose, Sens-Schönfelder, & Campillo, 2013;Yuan et al., 2021), which can be interpreted as a greater spatial sensitivity (e.g., of the unperturbed medium) with lower frequencies.Interestingly, the rate of long-term increase is also observed to increase with frequencies (Figure S9 in Supporting Information S1).Yuan et al. (2021) showed that only uniform (i.e., depth independent) change in velocity affects frequencies equally.Therefore, if the variation of dv/v associated with the seasonal effects is uniform with frequencies, while the rate of long-term increases with them, it has a potential to separate those effects such that the tectonic transient signals and long-term in dv/v are shallow perturbations, and seasonal effects affect a greater depth range.
We use the frequency band of 0.9-1.2Hz to conduct the model fitting and the comparison to the strain in the following sections.The ∼1 km depth of structural change corresponding to the frequency band is inferred from the Rayleigh wave depth sensitivity kernel (Figure S10 in Supporting Information S1).Note that both the Love and Rayleigh waves could contribute to the channel-weighted dv/v, whereas only the Rayleigh wave can affect the vertical component shown in Figure 4.

Multi-Factor Model
In order to investigate the source of the observed long-term increase in dv/v, we fit the time series with several different models by decomposing the time series of dv/v at each channel-weighted station pair into different model components using the observed environmental factors near the Parkfield section of the SAF.In California, there have been numerous reports of seasonal effects on shallow measurements of dv/v (Hillers et al., 2015; We used the monthly precipitation time series recorded at Parkfield with the Remote Automated Weather Stations (RAWS) provided by Western Regional Climate Center.We resampled every 15 days by linear interpolation to synchronize with the time history of dv/v.
The time series of temperature at Parkfield is downloaded from the National Oceanic and Atmospheric Administration (NOAA).We removed the mean offset from the daily averaged temperature time series to obtain the temperature anomaly.We then applied the low-pass filter with a cut-off period of 20 days and downsampled it along the time history of dv/v.

Base Model
To model the dv/v time series, the non-tectonic factors, such as the precipitation and the thermoelastic strain caused by changes in atmospheric temperature, are included.We develop a dv/v time series model associated with the non-tectonic factors along with the coseismic decrease in dv/v and the logarithmic healing model proposed by Snieder et al. (2017) which we refer to as "the base model."The base model is formulated as where the parameters with numerical indices (i.e., 0, 1, and 2) indicate parameters which are fit to the dv/v time series and are defined in Table 2.
The first term a 0 is the constant average level of the dv/v time series.The second term describes the hydraulic effects on the dv/v time series and is comprised of a proportionality constant p 1 times the change of GWL (ΔGWL).The change in GWL is derived from a model of the pore-pressure diffusion of the observed rainfall (Akasaka & Nakanishi, 2000;Sens-Schönfelder & Wegler, 2006), which excludes storage in large spatially confined reservoirs such as aquifers and lakes, such that where p(t i ) is the precipitation at time t i , ϕ is the porosity, and α 0 is an exponential factor, which controls the hydraulic decay rate.The summation over the index n indicates that the value of the change in GWL at time t i is dependent on the previous time steps of the precipitation time history.The ΔGWL time series is then calculated using the RAWS time series of precipitation to estimate p(t i ) and the values of α 0 and ϕ as schematized in Figure 6a.Due to the less constraint in α 0 from the dv/v time series, we fixed it during the MCMC analysis described in Section 4.4.We also fixed the porosity ϕ to be 5%; however, any error in the assumed porosity is automatically absorbed by the parameter p 1 .The best fitting model factor p 1 is then estimated from the dv/v time series using this change in GWL time series ΔGWL(t i ).The factor p 1 is constrained to be negative to reflect the physical constraint that the dv/v time series values will decrease for an increase in ΔGWL.
The third term is associated with the atmospheric temperature (Figure 6b).Berger (1975) and Ben-Zion and Leary (1986) described the thermoelastic response at depth due to the slow diffusion of a temperature change at the surface (e.g., with time delay).We simplified this effect by shifting the time history of temperature by t shift 0 with the constant factor p 2 .Note that we constrain the t shift 0 to be positive to meet the causality of the thermoelastic deformation.
The fourth and fifth terms show the coseismic velocity decrease and the healing associated with the SS and PF earthquakes, respectively.The constants of proportionality of s 1 and s 2 , which relate the modeled healing to the observed changes in dv/v.We used the logarithmic healing model of Snieder et al. (2017) given by L t,τ min ,τ max ,t EQ ) = where the τ min is the minimum relaxation time corresponding to the initial healing rate, τ max indicates the period when the healing is completed, and t EQ corresponds to the occurrence time of the earthquake (Figure 6c).The magnitude of coseismic decrease s i L t = 0,τ min i ,τ max i ,t EQ ) is equivalent to s i ln τ max i / τ min i ) , which can be related directly to the observed coseismic decrease in the dv/v timeseries.We improved the computational efficiency of the numerical integration using the pre-compiled library as described in Text S3 in Supporting Information S1.

Model With Linear Trend
We synthesized the model with a linear trend term such that where b 0 is the slope of the linear trend (Figure 6d).Note that the intersection of the linear trend is included in the term a 0 of the base model y base .The linear trend term is applied over the entire study period such that we assume it is caused by the background tectonics or some other factors affecting the dv/v on the entire period.Ikeda and Tsuji (2018) included this linear trend term in their model as they found the long-term increase in the dv/v from the observation near the Nankai trough off the Kii Peninsula, Japan, which is discussed in Section 6.2.Ermert et al. (2023) also included a necessary trend correlated with Mexico City basin subsidence.

Model With Residual Healing
We considered an alternative model that includes a residual healing term, which assumes that the healing from other earthquakes that occurred prior to our observational period affects our perceived baseline dv/v values before the SS earthquake.Potential candidates from the USGS catalog are the M6.7 Coalinga earthquake on May 1983 (e.g., Stein & King, 1984), the M4.8, 8 km NW of PF on November 1993 with a source depth of 10.9 km, and the M4.9, 3 km NW of PF on December 1994 with the depth of 8.3 km.The magnitude of the latter two earthquakes is relatively small, but their location is close to the seismic stations.Thus, if the coseismic decrease in dv/v occurred due to those events, and the healing still took place until the SS earthquake, a residual of the healing could be included in the dv/v.This situation is analogous to the discussion in Vidale and Li (2003) for the interruption of velocity healing due to the 1999 Hector Mine earthquake after the 1992 Landers earthquake.
To formulate this model, we added the term associated with the residual healing to the base model such that where c 0 and H are the factors of residual healing and the Heaviside function at the date of the SS earthquake, respectively (Figure 6e).

Methodology of MCMC Analysis
To conduct the model fitting with the time history of channel-weighted dv/v, we used the Python-based software tool emcee (Foreman-Mackey et al., 2013).It provides the modules of the MCMC method with various advanced sampling algorithms.We selected the one from them called the stretch move proposed by Goodman and Weare (2010), which updates the model parameters using a set of walkers.We set the number of walkers and the steps of iterations as 32 and 20,000, respectively.The log-likelihood function with a set of model parameters θ is defined as: where dvv is the estimated dv/v time series, and σ2 n = σ 2 n + f 2 0 with the σ n error of the dv/v estimation, and f 0 is a model parameter associated with the additional variance of the dv/v measurement, which is searched during the MCMC analysis similar to the other model parameters.Note that we use a different form of σ2 n from what is documented in emcee as σ2 n = σ 2 n + f 2 0 y 2 model such that the additional variance is added uniformly over the study period (L.Ermert, pers. comm., 2022).We performed the model inversion using MCMC to evaluate the contribution of each model parameter as well as the trade-off between the parameters.The lower and upper bounds of the model parameters are listed in Table 2.
During the preliminary MCMC parameter search, the largest trade-off is observed between the parameters s i and τ min i (i = 1, 2), degrading the model's convergence as summarized in Text S4 and Figure S11 in Supporting Information S1.We thus computed the median values of the maximum likelihood parameter associated with τ min i for the station pairs in the preliminary case study and used them as representative values to fix the model parameters for the present analysis.We also fixed α 0 as some of the station pairs are less sensitive to the ΔGWL.The values of the fixed parameters are listed in Table 3.Note that the fixed value of α 0 is larger than the other applications (e.g., Sens-Schönfelder & Wegler, 2006), which might indicate higher drainage conditions.
The parameter search range associated with the p 1 , s 1 , s 2 , and c 0 is constrained in the positive or negative side for the consistency with the sense of change in dv/v.The bounds of t shift 0 are set to be up to 90 days due to the trade-off between the p 2 and the days of shift such that the time shift of half a year with the flip of the sign shows almost identical time series considering the seasonal variation.We also set a threshold such that the magnitude of s 1 associated with the SS earthquake is to be smaller than half of s 2 with the PF earthquake.The more adequate criteria could be the comparison with the coseismic velocity decrease in the logarithmic healing model, that is, s i ln τ max i / τ min i ) , whereas we used the former criteria to simplify the constraint condition.

Comparison of Best Likelihood Model to the Data
We compared the model fitting associated with the base and linear term models for all the available station pairs.We performed the statistical analysis to show the contribution of the linear trend term in the dv/v as well as the distribution of the other model parameters.We also conducted the MCMC analysis using the model associated with the residual healing term for a representative station pair to investigate its role in the dv/v.
Figure 7 shows the marginalized 1D and 2D posterior probability distributions of the model parameters fit using the linear model to the dv/v time series obtained using the stretching method for the channel-weighted pair of LCCB-SCYB.The marginalized 1D posterior probability distributions of all the model parameters except τ max 1 are strongly peaked and unimodal, indicating that each of the parameters, and their error estimates, are well constrained (gray histograms, Figure 7).The 2D marginalized posterior probability distributions associated with τ max 1 (red 2D histograms, Figure 7) reveal a tradeoff between its value and the value of each of the other parameters (i.e., much larger covariances).The parameter τ max 1 determines the period of healing following the SS earthquake and is not well constrained due to the short observational period T SS obs between the occurrence of the SS and PF earthquakes (i.e., T SS obs ≪ τ max 1 ).We note that the 1D marginalized posterior probability distributions indicate that the parameters s 2 , τ max 2 and b 0 are well constrained, despite the trade-off between parameters pairs (e.g., s 2 -to-τ max 2 , s 2 -to-b 0 , and τ max 2 -to-b 0 ) indicated by their 2D marginalized posterior probability distributions.We omitted to show the uncertainty f 0 as it shows the unimodal distribution, not interfering with the other model parameters.
The results obtained with the MWCS-derived dv/v time series (Figure 8) are similar to those obtained using the stretching method.The main difference is in the parameter t shift 0 , which is shifted from ∼30 days to near zero.This discrepancy between the stretching and MWCS results is likely due to the weighting of linear regression in the MWCS, which could decrease the sensitivity to large fluctuations of phases from temperatures.
Figure 9 compares the observed dv/v time series and the model with the set of best likelihood parameters.Each row shows the contributions of model factors to the dv/v and the preprocessed time series of precipitation and temperature.The dv/v models are synthesized using the maximum likelihood parameters, which are the set of parameters that achieve the best likelihood during the MCMC sampling.For the station pair LCCB-SCYB, the best fitting models to both the stretching and MWCS-derived channel-weighted dv/v time series show that the long-term increase after the PF earthquake is better reproduced with the linear trend term (Figures 9 and 10).

Statistical Analysis of Model Parameters
We estimated the best likelihood parameters for the available station pairs, which can be found in Data Set S5.We show the model fitting of all the station pairs in Figures S7 and S8 in Supporting Information S1 for the cases with the stretching and MWCS methods, respectively.We quantified the convergence of the model using the variance of the residuals such that we selected the station pairs with the variance of residuals smaller than the threshold of 0.002.
We also removed the station pairs and ACFs associated with the FROB, as it shows the strong, suspicious harmonic noise in the raw data, which can cause bias in the CFs.We finally used the 32 and 29-channel-weighted station pairs for the cases with the stretching and MWCS, respectively, for the statistical analysis of model parameters.
We computed the two metrics to evaluate the quality of models: the Akaike information criterion (AIC, Akaike, 1974) and the Bayesian information criterion (BIC, Schwarz, 1978) where N is the number of data, and k is the number of model parameters.Figure 11a compares the AIC and BIC values for the case with the two methods of dv/v measurement.ΔAIC = AIC linear trend AIC base , and same for ΔBIC.Most of the station pair shows the negative ΔAIC and ΔBIC, which supports the additional model term and model complexity.Thus, the linear trend term shows a nonnegligible contribution to the dv/v.Figure 11b shows the estimated value of the slope b 0 %/year.The median for the cases with the stretching and MWCS are 0.0048% and 0.0027%/year, respectively.The coefficient with the temperature p 2 varies in both negative and positive values.In contrast, the median is positive, indicating the increase in temperature corresponds to the rise in dv/v, which is consistent with the model of Richter et al. (2014).The negative sensitivity of dv/v to the temperature could be caused by the artifacts of the model fitting to the station pairs, which are less sensitive to the temperature, while it remains to be identified.G.
Li and Ben-Zion (2023) showed the seasonal variation of the velocity change ranging from 0.4% to 1.2%, which is more significant than our case.This would be caused by the different frequency ranges-1-6 Hz in their case-, while we use 0.9-1.2Hz.
Maximum healing times τ max i are better constrained with the base model than the model with linear trend, which could be caused by the trade-off between τ max i and b 0 (Figures 7 and 8).However, the residuals between the model and data are smaller, with the case using the linear trend term for most of the station pairs.Instead, we need to subtract the long-term increase from the dv/v, if it is plausible, to adequately evaluate the maximum time of the logarithmic healing after the earthquakes.We generally find τ max 2 to be about 10 years, a similar scale to that obtained by Clements and Denolle (2023) in southern California earthquakes.
The other parameters, such as p 1 and t shift 0 , are more scattered when incorporating the linear trend term.These parameters may be less constrained due to their small sensitivity to the dv/v: the effects of precipitation and temperature are not evident in the channel-weighted time history of dv/v in Parkfield (e.g., the seasonality is weak) compared to the other areas (Clements & Denolle, 2023;Sens-Schönfelder & Wegler, 2006).
In general, the non-tectonic parameters are not well constrained since some station pairs do not exhibit strong seasonal signals (e.g., low sensitivity of dv/v to these factors).The posterior distributions of model parameters are narrower for the base model (no linear trend), whereas the addition of the trend term widens the posterior distribution, even if it is well justified by the negative ΔAIC and ΔBIC.

Model Fitting With the Residual Healing
We performed the model fitting with the residual healing term for the LCCB-SCYB as a representative pair of dv/ v time history showing the logarithmic healing and the long-term increase.Figure S12 in Supporting Information S1 shows the marginalized posterior probability distributions for the case with the residual healing model using the MWCS method.The residual healing coefficient c 0 shows the trade-off with a 0 , the static term, such that the a 0 increases with the magnitude of c 0 to fit with the period before SS. Figure 12 shows the contributions of the model factors and the comparison to the data.We found c = 0.05% of the residual healing before the SS earthquake and τ max 2 = 25 years in this station pair.The fitting of the model could reproduce the channel-weighted dv/v by complementing the long-term increase in dv/v with the combination of c 0 , a 0 , and the logarithmic healing.
The AIC and BIC are smaller in the case with the base model, whereas they are larger in the case with the linear trend term shown in Figure 10, suggesting that the model with linear trend term is preferred over the residual healing for this station pair.However, the comparison of the cases with the other station pairs might show different results, and we eventually cannot evaluate how much the comparison of the AIC and BIC between the models with linear trend and residual healing terms is statistically significant as both could reproduce the longterm increase in the dv/v.Thus, the choice of models with the linear trend or the residual healing terms remains to be determined from the model fitting analysis.Instead, if the residual healing model is suitable, we could deduce that the steady state of the velocity corresponding to the a 0 is larger than the current state of dv/v, and the logarithmic healing takes place at least more than 18 years.

Variation of the Linear Term With Fault Normal Distance
We analyzed the slope of the linear trend term b 0 as a function of the fault normal distance to investigate the spatial characteristics of the long-term increase in dv/v.We computed the fault normal distance of the stations by defining the approximated planar fault along the SAF, as shown in Figure S1 in Supporting Information S1.We averaged the distances for the two stations in the cases of CCFs.The b 0 is obtained with the best likelihood parameters.method than the MWCS, which could be caused due to the characteristics of the methodology in the measurement of the dv/v as the MWCS is more likely to underestimate the magnitude of dv/v (e.g., Hillers et al., 2019).Except for those outliers described below, the distribution of b 0 is relatively uniform, or weakly anti-correlated, with the fault normal distance.The effect of the shear localization or the damage zone in the fault core, as inferred by Delorey et al. (2021), is not clearly apparent in b 0 , even though some station pairs near the fault show relatively large values.To investigate if the observed dv/v variations revealed by the variation in b 0 with fault distance reflects the velocity perturbation in and around the fault core, the sensitivity kernels would need to be evaluated carefully (Obermann, Planès, Larose, Sens-Schönfelder, & Campillo, 2013).The high-frequency body waves with the known noise sources would be able to improve the spatial resolution of the dv/v (Sheng et al., 2022) as the direct P wave can be extracted from the noise correlations (e.g., Roux et al., 2005).Some station pairs near the fault show larger values in b 0 relative to station pairs far from the fault (Figure S13 in Supporting Information S1); however, the b 0 with GHIB would be infeasible as the associated time history of dv/v is unstable with larger variances than the other stable station pairs (Figures S7 and S8 in Supporting Information S1).The negative b 0 obtained for station pair LCCB-VCAB indicates a long-term decrease, which is not visible in the dv/v time history (Figure S7 in Supporting Information S1).This negative b 0 value is likely caused by the trade-off in the model fitting between b 0 and the healing term of the PF earthquake rather than the negative trend in dv/v.Note that the slope for this station pair is only obtained for the stretching method, and a positive value is estimated using the MWCS-derived dv/t time series.

Cumulative Strain Field at Parkfield
In order to investigate whether the origin of the long-term increase in dv/v is associated with the accumulation of regional strain on the SAF near Parkfield, we estimate the regional strain at Parkfield from 2009 to 2022 using GNSS data.
The dilatational strain in rock is thought to affect the observed dv/v time series by altering the average wave speeds of the rock with contraction and extension increasing and decreasing the velocities, respectively.This relation is formulated in nonlinear elasticity (Ostrovsky & Johnson, 2001) and is supported by observations (e.g., Brenguier, Shapiro, et al., 2008;Donaldson et al., 2017;Richter et al., 2014;Yamamura et al., 2003).Typically, the sum of axial strain components, that is, dilation, correlates with dv/v observed with the ambient seismic noise (Donaldson et al., 2019).However, Hotovec-Ellis et al. ( 2022) showed the velocity can be sensitive to a single component of the strain rather than the dilation when the pre-existing cracks are oriented in the suitable direction.In the case of Hotovec-Ellis et al. ( 2022), the cracks are generated perpendicular to the radial direction of a caldera, forming the vertical ring fractures such that the opening and closing of them are governed by the radial strain component.Thus, if the cracks are aligned in a preferable orientation, the velocity could be effectively changed with the extension or contraction of the axial strain (Hotovec-Ellis et al., 2022).
The dv/v should be evaluated with the spatial average of the interaction between the dominant crack orientations and the cumulative axial strain, as well as the efficiency of the crack opening and closing with the strain (Sayers & Kachanov, 1995).In this study, however, we computed the spatial average of azimuthal axial strains around the seismic station to investigate the contractional directions of the strain rather than encompassing all the spatial variation of scattered crack orientations with depth-varied maximum compressive stress S Hmax (Hickman & Zoback, 2004), which is less feasible in the resolution of the observations and the uncertainty of the cause of anisotropy.

Processing of the GNSS Data to Estimate the Temporal Evolution of Strain Field
We downloaded the GNSS data processed by the NASA MEaSUREs ESESES Project (Bock et al., 2021).We used the daily displacement time series in the region of western North America, where the outliers and nontectonic jumps have been removed.
We computed the in-plane strain with the triangular elements with the nodes of GNSS stations generated by the Delaunay triangulation.The deformation of an element caused due to the plane strain without the rigid motion is written as follows: where u k x and u k y (k = 1,2) are the displacements of the station.
x k and y k are the relative locations of nodes from the origin that is chosen from one of the nodes in the element.ɛ ij and ω are the strain and the rotation, respectively.This representation of deformation can be found in Shen et al. (1996) and Crowell et al. (2013).Given the displacement of stations, we obtain the strain tensor, which is assumed to be constant in the element.The convention of the sign is positive in extension for the axial strain components and the dilation.We used the software tool PyTAGS (Crowell, 2019) to conduct the mesh discretization with triangular elements and compute the strain components associated with the elements.
We first obtained the strain components oriented to the East and North, corresponding to the x and y directions in Equation 11, respectively.We then computed the dilation ɛ 1 + ɛ 2 and the max shear (ε 1 ε 2 )/2 of the element, where the ɛ 1 and ɛ 2 are the maximum and minimum principal strains, respectively.
Figure 13 shows the temporal evolution of the dilation and max shear.The max shear is positive and localized on the fault, while the dilation is not uniformly distributed along the fault.It should be noted that the discontinuity (i.e., slip) crossing the element causes bias in estimating strain as the element is assumed to be a continuum, and the strain is constant within the element.Therefore, the estimated strain near the fault can be apparent, largely distorted by the slip on the fault.One of the potential causes for the intricate pattern in the dilation could be the localized slip on the shallow surface (Bacques et al., 2018).
The magnitude of cumulative strain far from the fault is in the order of 1μɛ over 10 years.Considering the macroscopic strain rate reflects the region far from the fault, the strain rate should be in the order of 0.1μɛ/year, which is consistent with other studies (Klein et al., 2019;Ross et al., 2022), although the detailed strain distribution cannot be compared due to the different spatial resolutions of the aforementioned studies.

Time History of Cumulative Strain With Gaussian-Weighted Average Around the Seismic Station
To evaluate the sensitivity of velocity change to the strain by comparing the long-term increase in the dv/v with the cumulative strain associated with the dilation and the max shear, we computed the Gaussian-weighted average of the strain field associated with the seismic stations.We only used the channel-weighted dv/v associated with the ACFs.We estimated the time history of the cumulative strain using the weighted average of the strain snapshots.The weight is computed with the numerical integration of two-dimensional Gaussian distribution over the surface of the triangular element as follows: where w i is the weight on the ith triangular element, Ω i is the surface of the triangular element, and ξ is the relative coordinate from given seismic station.σ defines the spatial extent of the weight around the seismic station.We set σ = 5 km.We normalized the weights before computing the weighted average such that the sum of the weights is unity for the cases where the strain data is missing on some triangular elements due to the lack of GNSS data.
We applied the threshold of the missing data of GNSS near the seismic station to exclude the case if the cumulative strain is estimated with the triangular elements only far from the station.We skipped the time step of cumulative strain if the triangular elements within the distance of 5 km were missing due to the gap in GNSS data.

Comparison of Long-Term Increase in dv/V to the Cumulative Strains
We pre-processed the channel-weighted dv/v before comparing them with the time history of the cumulative strains.We subtracted the model components of ΔGWL and temperature with the best likelihood model parameters from the dv/v and removed the offset with the mean between 2008 and 2010, considering the reference time of the cumulative strain.Indeed, we could also remove the logarithmic healing terms from the dv/v to focus on the linear trend.However, the τ max 2 shows the trade-off with the b 0 as shown in Figures 7 and 8, which could cause the overestimation of the long-term increase in dv/v.Therefore, the logarithmic healing terms are kept in the dv/v to avoid bias due to subtracting them.To compare the strain and dv/v, we excluded the stations of GHIB, JCSB, and VARB as the measurement of the dv/v is unstable compared to that of other stations (see Figures S7 and S8 in Supporting Information S1).Therefore, we used the seven stations (CCRB, EADB, LCCB, MMNB, SCYB, SMNB, and VCAB) for the following analysis.We also set the period of comparison until 2020, as many GNSS data are missing after that.
Figure 14 compares the dilation and max shear to the dv/v of ACFs for the cases using stretching and MWCS, respectively.We compared the cumulative strain and the dv/v associated with the ACFs of the seven stations, where the estimations of dv/v were relatively stable from 2009 to 2020.The time steps are synchronized between the cumulative strain and the dv/v, and we compared them every 90 days from the reference time snap.
The Gaussian-weighted time history of the cumulative dilation around the seismic stations shows the subtle extension with time on average of the seven stations.The comparison with the dilation shows weak positive correlations to the dv/v, which is the opposite relation of what one would expect from nonlinear elasticity or from the opening of microfracture in micromechanics.Note that the dilation associated with EADB shows slight contractional accumulation, whereas the magnitude of strain is smaller than the other stations, which would be The lines in the elements show the orientations of the maximum contractional strain.The color contour shows the amplitude of cumulative strains, which is obtained by subtracting the initial state of strains on 3 January 2009, from the current time snapshot.Max shear can be negative as it is the relative value from the initial state.As we evaluate the strain field with the assumption of the continuum in the elements, the localized slip can cause bias in the estimation of strain.
insufficient to evaluate the sensitivity with the dv/v.The similar discrepancy in the sense of dv/v and the dilation is also addressed in Hotovec-Ellis et al. (2022).
We thus focused on the azimuthal axial strain and computed their Gaussian-weighted average around the stations.Instead of summing up all the local strains along S Hmax , whose orientation varies quite significantly with a range of ±45°as shown in Figure 13, we uniformly rotated the strain around the stations to obtain the first-order evaluation of the axial strain to investigate if it shows contraction to cause the closing of given aligned cracks, which increases the dv/v.
Figure 15a shows the mean and standard deviation of the station-averaged axial cumulative strain as a function of azimuth.The cumulative strain is evaluated between the reference in January 2009 and November 2019.The estimation of maximum contractional orientation inferred from the GNSS data is consistent with the range of S Hmax estimated by Hickman and Zoback (2004).It is notable that the azimuthal strain can be in contraction in the range from N35°W to N45°E although the mean value of axial strain, which is equivalent to the half of the dilation ɛ 1 + ɛ 2 , shows the extension.b) max shear to the dv/v with the frequency band of 0.9-1.2Hz.We used the dv/v of the seven stations listed in the legend, while we excluded the GHIB, JCSB, and VARB due to the low quality of the dv/v measurements.We compared the strain and dv/ v at the synchronized time periods every 90 days.The reference slopes are indicated for the sensitivity of dv/v to the strain.The slope and standard error of the linear regressions using all seven stations are annotated in the panels.
We performed a linear regression between the cumulative strain and the dv/v averaged over seven stations to obtain the sensitivity at each azimuth, which is shown in Figure 15b.Figure 15c shows the estimations of negative sensitivity evaluated at N5°E, which are 0.011 ± 0.001%/μɛ and 0.007 ± 0.001%/μɛ for the cases with stretching and MWCS, respectively.
The sensitivity of dv/v to the strain component is smaller than reported previously, as generally referred to in the order of (0.1-1.0)%/μɛ.Although most studies derive the strain sensitivity to the dilatational strain rather than a single extension strain component, and the condition of the tectonic setting, the frequency range of the CFs, and the installation of seismometers vary, we can still compare with other values found in the literature, which is summarized in Text S5 in Supporting Information S1.
Further investigation is necessary to fully explain the magnitude and direction of sensitivity of dv/v to strain.Nonetheless, the preferable interactions between the aligned microcracks, which could be non-uniformly distributed at Parkfield as inferred from the shear wave splitting analysis (see Text S6 in Supporting Information S1), and the contractional strains could be a candidate to explain the long-term increase in dv/v on the multiple seismic station pairs.

Potential Artifacts in dv/v
In this section, we list the potential artifacts to the dv/v caused by instrumental issues rather than tectonic or environmental origins before discussing the factors associated with the long-term increase in dv/v.Zhan et al. (2013) shows the potential of the apparent change in the velocity due to the evolution of the noise source spectrum.This artifact is still controversial as Mao et al. (2020) conducted the numerical experiments of the coda interferometry with multiply scattered waves using the noise sources with different source spectra, indicating the measurements of dv/v are relatively insensitive to the perturbation of the noise source spectrum.

Change of Noise Source Spectrum
Over time, sensors get updated, and instrumental responses may change.We find an increase in the power spectrum at higher frequencies around 2014, possibly due to a change in the sensitivity of stations as discussed in Sections 2.1 and 3.This change appears suddenly, and it might be correlated to the bump of the dv/v shown in the case of stretching in Figure 3, but it cannot explain a long-term trend.
In our case, if the noise source perturbation is the origin of the apparent long-term increase in the dv/v, the noise source spectrum should gradually change over 20 years, uniformly causing the long-term increase in dv/v on the difference station and component pairs.We remark that the seasonality of microseismic noise is not strong at the frequencies of interest in our study around 1 Hz, as shown in the supplementary materials associated with the PSD.We need to investigate if the variation of source spectra can cause a uniform change in the measurement of dv/v over multiple station pairs.Stehly et al. (2007) pointed out that clock drift causes the artifacts in the travel time estimation obtained with the stacked CFs.Various methods have been proposed to correct the time shift (e.g., Gouédard et al., 2014;Hable et al., 2018;Hirose & Ueda, 2023;Sens-Schönfelder, 2008).In our data set, we found the time shift in daily CFs associated with the EADB around 2017 in the plot of CCFs (see the CCFs for EADB-SCYB in Data Set S2), although it seems to be corrected within the same year.We also found the transient time shift in CCFs associated with the GHIB around 2016.

Clock Drift
We can remove the bias of the clock shifts to the dv/v measurement by correcting them as implemented in for example, Brenguier, Campillo, et al. (2008).Instead of applying the proposed correction method to our CFs, we used the MWCS such that we conducted the linear regression without imposing the intersect passing the origin associated with the lag time and the time delay to mitigate the artifacts due to the clock drift.This metric is also implemented in the MSNoise (Lecocq et al., 2014).The temporally shifted CFs cause the artifacts in the cases using the stretching method, as the origin of the stretching should be corrected.The dv/v using the MWCS with imposing the intersect crossing the origin is also biased due to forcing the linear regression to the time delay with non-zero intersect.On the contrary, we could minimize the effect of clock drift without imposing the intersect crossing at the origin with MWCS because the clock drift causes the shift of the intersect but does not change the slope of the time delays, reflecting the structural change of the velocity (Gouédard et al., 2014), which we used in this study.Moreover, it is worth noting that the clock would not influence the results with the ACFs drift as the rate of time shift should be negligible compared to the coda length even if the clock drift occurs.Thus, the clock drift would not be the dominant factor to characterize the linear trend shown in the model fitting.
It should also be noted that the distinctive decrease in the dv/v is indeed observed in the case with both stretching and MWCS around 2017, whereas it also coincides with the large precipitation as shown in Figures 9 and 10.Therefore, it can reflect the velocity change of the structure rather than purely the artifacts of the crock drift.

Phase Shift in the Data Sampling
Before replacing the location codes and channel names around 2011, we noticed the sampling time of the raw waveform was not aligned on a uniform time vector, that is, 0, 0.05, 0.1,…[sec].To compute the correlation functions, we shifted, rather than resampled, the waveform to the closest sampling point on a uniform time vector.According to our tests using both obspy and SeisIO (SeisIO.jl is currently maintained as SeisBase.jlat https:// github.com/JuliaSeismo/SeisBase.jl) download functions, the value of the phase shift seems arbitrary with each file, though it is always less than half of the sampling rate (25 ms).Thus, the dv/v measurements using the stretching method could be biased by the phase shift.However, similar to the discussion about clock drift, the artifacts due to the phase shift are mitigated in the cases with the MWCS.The phase shift of the data set has been corrected after replacing the location codes and channels so that the artifact is not the case for the later part of the dv/v, from November 2010 to September 2011.
We note that some of the daily stacked CCFs lack the first half slice of the short-time cross-correlation windows with a length of half an hour due to the fluctuation of the initial time of the data samplings in the present analysis.While this needs to be corrected for the study's completeness, it would not primarily affect the daily and monthly stacking of the CCFs.

Change in the Instrumental Response
We must be careful about the long-term changes of instrumental responses as we discuss the dv/v over decades.Ueno et al. (2015) showed the pseudo perturbation of the CFs caused by changes in the instrumental response, which has led to biases in the estimation of dv/v.The instrumental response available at the Northern California Earthquake Data Center (NCEDC) is the best sensor and digitizer metadata knowledge.Reassessing this response may give sudden sensitivities changes, as seen in Figure S2 in Supporting Information S1.However, we cannot correct the minor change in the response due to site conditions, such as sensor-to-ground coupling changes or the instruments' aging.
The large fluctuation in the instrumental response is more likely to be thresholded out by the median filter applied during the computation of CFs or the quality of the dv/v estimation during the post-processing.In contrast, the minor change cannot be corrected and should be included in the analysis.However, identical to the discussion of the noise source perturbation, if the difference in the instrumental response causes the long-term increase, the CFs should be distorted such that the changes in the time delay consistently cause the increase in the dv/v over the different station and channel pairs.We would need to investigate instrumental aging when evaluating the longterm dv/v over decades.
Figure 15.The rotated axial strain and the variation of sensitivity with the azimuth.(a) Cumulative axial strain parallel to the given rotation angle evaluated between January 2009 and November 2019.We computed the Gaussian-weighted axial cumulative strain parallel to every five degrees as the circles.The filled band shows the standard deviation of the strain associated with the seven seismic stations used in Figure 14.The vertical dashed lines show the parallel and normal orientations to the San Andreas Fault.The highlighted area shows the range of S Hmax at depth from 0.8 to 2.2 km estimated by Hickman and Zoback (2004).The horizontal dashed line indicates the mean of axial strain over the rotated angle.(b) The sensitivity of dv/v with the frequency band of 0.9-1.2Hz to the axial strain.The solid lines with circle and square markers indicate the sensitivities associated with the stretching and moving window cross-spectral, respectively.The error band indicates the 95% confidence interval of the sensitivity associated with the linear regression of the slope.We excluded the azimuth between ±30°-60°as the sensitivity diverged due to the small axial strain.The dotted line shows the rotated axial strain similar to (a).(c) Sensitivity of the velocity to the axial strain nearly parallel to the S Hmax .The negative sensitivity indicates the dv/v increases with the contractional strain.

Potential Factors to Explain the Long-Term Increase in dv/v
While the linear trend would appear due to non-tectonic factors or instrumental issues, an accumulation of the axial strain can close the aligned and open micro-cracks in the medium, causing a positive dv/v.However, the contribution of the latter tectonic factor associated with the cumulative strain is still debatable.Our analysis did not pursue a 2D, improved spatial imaging of dv/v (Obermann, Planès, Larose, & Campillo, 2013) and our crude knowledge of micro-cracks orientation in the fault zone.Given that the dv/v estimated from the single-station correlations reflects the structural change around the station, the comparison of it to the spatially averaged strain accumulation shows the sensitivity of the dv/v to the strain averaged in space.Besides, the estimated strain field is limited to the two-dimensional horizontal components; the vertical variation of the strain is not considered in the analysis.Therefore, we cannot quantify the sensitivity of dv/v to strain.On the other hand, if the given pattern of regional strain accumulation and the aligned cracks cause the long-term increase in dv/v, we could explain the linear trend shown in dv/v associated with the multiple station-component pair combinations for both ACFs and CCFs.Further, the fact that our measurements of shallow and near fault-zone dv/v exhibit a long-term compression may come from the fact that the 2004 Parkfield earthquake did not rupture at the surface; therefore, the shallow portion may experience partial loading, even if small slow creep may reduce the loading rate (Bacques et al., 2018).
Several studies have found long-term increases in dv/v in the accretionary prism of the Nankai subduction zone.Ikeda and Tsuji (2018) found a linear trend in dv/v near the Nankai trench with the range of 0.01%-0.03%/yearduring a few years, which was extended by Tonegawa et al. (2022) over a decade.These linear trends contrasted with the western part of the accretionary wedge (onshore, further from the plate boundary), where no increase was measured.Authors interpreted such variations in terms of variable compression and fluid drainage from the subseafloor.Offshore observations are different than Parkfield's case, especially in geohydrology, lithology, and microcrack structure.Nevertheless, we draw an analogy for velocity perturbations in the interseismic period near a plate boundary.
The model of residual healing is one of the other candidates to reproduce the long-term trend in the dv/v.Indeed, the quality of the measurement of dv/v is lower before the SS earthquake, as shown by the significant variance in the estimation of the dv/v in Figure 3.The change of preamplifier gain documented in the supplementary material of Brenguier, Campillo, et al. (2008) can be one of the reasons for the change of the variation.Therefore, the duration of the analysis period and the quality of the measurement in dv/v are insufficient to discuss if the healing had occurred before the SS earthquake.The difficulty in quantifying the contribution of residual healing is that the logarithmic healing due to the slow dynamics can be accompanied by the other factors for the long-term trend described above, causing the trade-off.The laboratory experiments with the servo-control of the shear loading rate (e.g., Shreedharan et al., 2021) would have the potential to isolate the logarithmic healing due to the slow dynamics to evaluate the other factors contributing to the dv/v.Another potential scenario of the continuous increase in dv/v is that the baseline level of dv/v increased due to the strain change caused by the PF earthquake.Olivier, Brenguier, Campillo, Roux, et al. (2015) monitored the dv/v associated with the blast in the underground mine and found the change in the baseline level of dv/v due to the static stress change around the station.Our case can be analogous to the situation in their study.We computed the 3D dilatational strain distribution of the Parkfield caused by the coseismic slip associated with the PF earthquake (see Text S7 and Figure S15 in Supporting Information S1), suggesting it is not obvious if the spatially averaged coseismic strain change shows the contraction or extension around the seismic stations.We further need to constrain the shallow coseismic slips, which would have a large effect on the dv/v considering its depth sensitivity, and to untangle all the contributions from the postseismic slips (e.g., Barbot et al., 2009;Johanson et al., 2006), viscoelastic relaxation (Bruhat et al., 2011) and slow-slip events on SAF (Delbridge et al., 2020;Rousset et al., 2019) to investigate if the spatial average of the resultant strain around the seismic station shows contraction such that the baseline is elevated.Lecocq et al. (2017) monitored dv/v, focusing on the thermoelastic and the hydrological effects over 30 years in the region without major earthquakes or volcanic and geothermal activities.They found a long-term increase in the dv/v, which coincides with the change of thermoelastic strain induced by the rise in air temperature.Thus, the long-term change of thermoelastic deformation can be a candidate for the potential factor of the increase in dv/v.The rate of growth in dv/v, which is equivalent to the slope of b 0 in our model, was estimated as ∼0.001%/year in (Lecocq et al., 2017), while we obtained the higher rates of 0.0048%/year and 0.0027%/year for the cases with stretching and MWCS, respectively.It should be noted that the air temperature cannot be the dominant factor of the long-term increase in our case; the long-term increase in the air temperature is ∼0.5°C for 20 years, while its seasonal variation is ±10°C.If the air temperature causes the increase of dv/v with 0.05% for 20 years, the range of seasonal variation should be ∼1%, which is not observed in our measurement.Thus, if the long-term increase in dv/v is caused by thermoelastic deformation, another candidate could be the continuous change of in-situ background temperature at depth.The relative locations between sensors and seismic scatterers can be altered by accumulating the strains or the slip on the fault.Snieder et al. (2002) showed the averaged travel time over multiply scattered waves in the coda window is not influenced by the independently perturbed scatters or the shift of source locations.If this is the case in our analysis, the measurements of dv/v only reflect the velocity change of the medium.However, the conditions associated with the number of scatterers and the frequency ranges with our analysis could differ from the case where the averaged time delay is canceled out.Thus, we further need to investigate if the long-term shift in the seismic stations and scatters can cause the apparent dv/v.
The dynamic earthquake ruptures can generate coseismic damage in the medium around the fault (Andrews, 2005), which forms the flower-like structure of the damage zone with depth; that is, the damage zone area is wider at shallow depth with lower confining pressure, and it becomes narrower at depth (Ben-Zion et al., 2003;Y.-G. Li & Malin, 2008;Ma, 2008).The orientation of off-fault fractures is governed by the stress field around the crack tip (Poliakov et al., 2002;Rice et al., 2005).The angle of the off-fault tensile crack varies with the rupture velocity (Griffith et al., 2009), whereas they are activated near parallel to the maximum compressive horizontal stress to the fault (Okubo et al., 2019;Thomas & Bhat, 2018;Yamashita, 2000).Therefore, the cumulative contractional strain parallel to the S Hmax would not effectively close those cracks.However, the area of the coseismic damage zone would be smaller than the spatial extent of the scattered coda wave.The rupture of the PF earthquake caused a small amount of slip near the surface (Johanson et al., 2006), whereas given the S Hmax in the order of 100 MPa at 1-2 km depth (Hickman & Zoback, 2004), the damage zone width is estimated in the order of hundreds of meters around the fault plane (Okubo et al., 2019).Besides, Y.-G.Li and Malin (2008) also showed the damage zone width is ∼200m at SAF inferred from the fault-guided wave.Hence, the orientation of coseismic tensile cracks would not be dominant in the kilometric extent of the sensitivity associated with the scattered wave.
Notably, the coseismic change in the anisotropy during the PF earthquake has been investigated by Durand et al. (2011) using the rotation of correlation tensor.They showed the substantial rotation after the PF earthquake at station VARB and pointed out that opening and closing of the cracks in the upper layer of ∼3 km would play a role in the observed rotation.Löer et al. (2018) also explored the temporal stability of the anisotropy with beamforming analysis using the ambient seismic noise at Parkfield.While we focused on the long-term trend of the dv/v in the present study, the time history of the anisotropy over decades should also attract research interest.
The weighted average of strain could be largely biased due to the slip on the fault causing the apparent strain change within the triangular element.Thus, the evaluation of dilation shows a primitive estimation around the stations.Nevertheless, it was insufficient to explain the long-term increase of dv/v in the context of compression with the isotropic medium.We thus investigated the accumulation of axial strain as an alternative candidate associated with the strain and the change of dv/v.The distribution of crack orientations cannot be simply retrieved from the fast polarization of shear wave splitting as it can be caused by various effects, such as the aligned minerals and grains rather than the cracks.However, the scattered polarization directions revealed by Y. Liu et al. (2008) would indicate the distribution of cracks is likely to be non-uniform around the fault.The local velocity change would be determined by the interactions of the pre-existing cracks and the cumulative strain, and the dv/v reflects their spatial integration with the sensitivity kernels (Obermann, Planès, Larose, Sens-Schönfelder, & Campillo, 2013).Those mechanisms can potentially cause a long-term increase observed in the dv/v time history.

Conclusion
We monitored dv/v from 2002 to 2022 to investigate the temporal changes near the Parkfield Region of the SAF and constrain the evolution of the SAF throughout the interseismic period.Our analysis revealed that the time history of dv/v was modulated by multiple factors which affected the physical state of the fault and the Journal of Geophysical Research: Solid Earth 10.1029/2023JB028084 surrounding subsurface, which can be separated into two distinct groups: environmental impacts (e.g., temperature variations and changes in the GWL from precipitation) and tectonic phenomena (e.g., earthquakes and interseismic loading).The effect of both the SS and Parkfield earthquakes on the dv/v time series is clearly observable with a similar level of previously reported values (Brenguier, Campillo, et al., 2008;Wu et al., 2016).However, the larger observational period used in this study, as well as the improved quality of the dv/v measurement from using all 3 components, allows us to robustly determine the influence of additional and more subtle environmental and tectonic factors.
In order to achieve this enhanced quality and longer duration of the dv/v time history, which enable statistically robust measurements of the role of inter-seismic tectonic loading, we have developed a Julia-based software tool to efficiently compute all the auto-and cross-correlations of the 3 components, of the 13 stations in the HRSN operated by the UC Berkeley Seismological Laboratory.It provided around 380 station-component pair combinations with good quality to evaluate dv/v over the study period.This software package SeisMonitoring.jl, which enables efficient process parallelization and greatly reduces the processing time, is openly available and will enable the calculation of dv/v for previously prohibitively expensive computationally intensive processes.
Following the 2004 Parkfield Earthquake the dv/v heals logarithmically for almost all of the station pairs and shows a net long-term increase in which the current dv/v level is equivalent to, or exceeding, the state before the 2003 SS earthquake.This long-term increase is investigated using models that are fit to the derived channelweighted dv/v time series.The full posterior probability distributions for these are determined using MCMC sampling.In order to isolate the effect of the inter-seismic tectonic loading, we account for environmental effects such as temperature variations, precipitation load and diffusion, as well as the logarithmic healing model as the aftermath of two earthquakes.Models that include a long-term linear trend or residual healing provide a statistically robust improved fit to the data.Both the AIC and BIC metrics confirm that the long-term trend observed in the dv/v time history is a non-negligible factor, and our analysis prefers the inclusion of either the long-term trend, which we interpret to be associated with the inter-seismic tectonic loading, or the residual healing such that the logarithmic healing has not been completed yet.Whatever factor controls this long-term trend should be spatially uniform near the fault, as this increase is observed in multiple station-component pair combinations.
To investigate the role of inter-seismic tectonic loading on the long-term increase of dv/v, we analyzed the GNSSderived strain from 2009 to 2020 in the Parkfield Region of the SAF.The spatially weighted average of dilatational strain around the seismic stations shows a slight extension, which is the opposite relation of what we expect for the increase in dv/v.However, the rotated axial strain reveals compression in a range near the maximum contractional strain (azimuth of N35°W to N45°E).The observed increase in dv/v can be explained by the pre-existing microcracks aligned perpendicular to the contractional strains, which may efficiently be closed to cause an increase in rigidity.
This study shows that by minimizing the uncertainties associated with the environmental factors through modeling their effect on the observed dv/v time series and by mitigating the non-tectonic artifacts such as the perturbation of the noise sources or the clock drift potential through improved processing methods, the dv/v can be used as an effective tool to monitor and constrain the evolution of the physical state in an active fault zone.

Data Availability Statement
We maintain the Julia-based software SeisMonitoring.jl(Okubo, 2023c) by the continuous integration testing on GitHub to be executed in different machine environments.The input files to conduct the ambient noise processing and the Jupyter notebooks for the post-processing are documented in the GitHub repository Seis-Monitoring_Paper (Okubo, 2024) to reproduce the results from scratch.The intermediate outputs, including the CCFs, are accessible in the cloud storage DASway (Okubo, 2023a).The minimal working example of the SeisMonitoring.jlto conduct from downloading data to measuring the dv/v using a docker container is available in SeisMonitoring_Example (Okubo, 2023b).

Figure 1 .
Figure 1.Location of the borehole network at Parkfield.The yellow triangles show the 13 borehole stations deployed and maintained by the UC Berkeley Seismological Laboratory.The black solid line shows the surface traces of the San Andreas Fault, which is obtained from the Quaternary Fault and Fold database, USGS.The background topography is downloaded from USGS 3D Elevation Program (3DEP) Data sets from the National Map.The blue circles show the LFE family (Shelly, 2017).The red line shows the ruptured fault associated with the 2004 Parkfield earthquake.The dashed line indicates the main fault plane of the rupture associated with the 2003 San Simeon (SS) earthquake (Johanson & Bürgmann, 2010).The stars indicate the hypocenters of major earthquakes; the SS earthquake is obtained from McLaren et al. (2008), and the others from SRCMOD(Mai & Thingbaijam, 2014).The seismic station locations with their name can be found in FigureS1in Supporting Information S1.SN: South Napa earthquake, RC: Ridgecrest earthquake sequence.
continuous seismic waveform contains transient signals such as earthquakes, tectonic tremors, environmental signals, and instrumental noises.The transient signals degrade the stability of the stacked correlation function (CF) as they induce strong, coherent signals on the CF that are typically absent.Zhou et al. (2020) implemented the process of asynchronous temporal flattening, which mutes the time windows containing the large amplitudes caused by the transient signals to avoid the perturbation of CFs.In this study, we also muted the transient signals from the continuous waveform using the kurtosis and classic STA/LTA algorithms Figure S3d in Supporting Information S1 shows the waveform after removing transient signals.After transient removals, the noise amplitude is balanced enough to mitigate the artifacts of transient signals on the stacked CFs.

Figure 2 .
Figure 2. Number of station-component pair combinations after selecting good quality dv/v measurements.The top and bottom figures are associated with the stretching and moving window cross-spectral methods, respectively.The vertical dashed line indicates the dates associated with the San Simeon and PF earthquakes.The shaded area in gray indicates the unstable period of dv/v, which could be caused due to the clock drift.The bottom arrows annotate the reference period from January 2010 to May 2022.

Figure 3 .
Figure 3.Time history of the dv/v for all the station-component pair combinations with the frequency band of 0.9-1.2Hz.The color contour shows the number of dv/v measurements meeting the threshold within 0.02% of the dv/v bin with respect to the time bins.The solid tick and thin lines indicate the median and the first and third quartiles, respectively.The shaded area in gray indicates the period where the clock drift was observed.The red dashed lines are the dates associated with the San Simeon (SS) and PF earthquakes.The offset of dv/v values is shifted by the mean of the median dv/v time history for the period before the SS earthquake occurred.

Figure 4 .
Figure 4. Time history of dv/v with the different set of components.The solid lines and the highlighted bands indicate the median and the first and third quartiles, respectively.The shaded area in gray indicates the same with Figure 3.The single-station and cross-stations contain the dv/v measurements associated with the nine components of auto-correlation functions (ACFs) or cross-correlation functions (CCFs), while the vertical and horizontal components are for the case with both the ACFs and CCFs, respectively.

Figure 5 .
Figure 5.Time history of dv/v with different frequency bands using all components.The lines and highlighted bands indicate the same with Figure 4.The lowest frequency band of 0.2-0.5 Hz associated with the stretching method is removed due to the large fluctuations in the estimation of dv/v.

Figure 6 .
Figure 6.The schematic of fitting model components.The base model consists of the factors associated with (a) ΔGWL, (b) temperature, and (c) logarithmic healing.To better fit the long-term increase shown in the time history of dv/v, we added either (d) linear trend or (e) residual healing terms to the base model.

Figure 7 .
Figure 7.The 1D and 2D marginalized posterior probability distributions of the model parameters obtained from the Markov chain Monte Carlo analysis for the station pair of LCCB-SCYB with 0.9-1.2Hz for the case with the stretching method.The diagonal panels show the probability mass functions where the sum of bars is normalized to be unity.The vertical solid lines indicate the best likelihood value.The color contour shows the 2D histogram of the pair of model parameters.The darker colors indicate larger probabilities.The circles with blue show the best likelihood model parameters.Note that the best likelihood parameters do not always correspond to the peak of the probability distributions.

Figure
Figure 11c summarizes the other model parameters.The a 0 is estimated to be larger with the base model than the model with linear trend term.This is due to the complementary of the long-term healing by the positively large offset as shown in Figures 9 and 10.
Figure 11c summarizes the other model parameters.The a 0 is estimated to be larger with the base model than the model with linear trend term.This is due to the complementary of the long-term healing by the positively large offset as shown in Figures 9 and 10.

Figure 8 .
Figure 8. Same to Figure 7 for the case using the moving window cross-spectral derived dv/v time series.

Figure 9 .
Figure 9.The contributions of the model factors and the comparison to the data for the cases with the base model (left) and the model with linear trend term (right) using the stretching method.Panels (a-c, g-i) show the factors associated with ΔGWL, temperature, and logarithmic healing.Panels (d, j)  show the contribution of the linear trend term, which is fixed at zero for the base model.Panels (e, k) show the comparison between the time history of channel-weighted dv/v (black) and the model with the best likelihood parameters (red).The green dashed line indicates the offset of dv/v associated with a 0 .The mean removes the offset until the San Simeon earthquake of the observation.Panels (f, l) show the residuals between the best likelihood model and the observed dv/v time history.Gray and black lines indicate the raw and the smoothed time history of the residuals, respectively.

Figure 10 .
Figure 10.Same to Figure 9 for the case with the moving window cross-spectral.

Figure
Figure S13 in Supporting Information S1 shows the distribution of b 0 with the fault normal distances.We categorized the ACFs and CCFs in the Pacific and North American sides and the station pairs crossing the fault, following Malagnini et al. (2019) and Delorey et al. (2021).The value of b 0 is generally larger with the stretching

Figure 11 .
Figure 11.Statistical analysis of the model parameters.(a) ΔAIC: AIC linear trend AIC base and ΔBIC: BIC linear trend BIC base for the cases with the stretching and moving window cross-spectral (MWCS).The box plot shows the first and third quartiles with the median value indicated by the horizontal bars in the box associated with the available station pairs meeting the variance threshold with the residuals.The black dots superimposed on the box plot show the individual values.(b) The slope of the linear trend b 0 .(c) The best likelihood model parameters.In each panel, we showed the statistics of parameters associated with the base model and the model with the linear trend for the cases with the stretching and MWCS, respectively.

Figure 12 .
Figure 12.Same to Figure 10 for the case with the residual healing term.Note that the a 0 is relatively large to complement the long-term increase of dv/v.

Figure 13 .
Figure 13.Cumulative strains associated with (a) the dilation and (b) the max shear.Note that the dilation is positive in extension.The triangular elements are formed with the nodes of Global Navigation Satellite System stations indicated with the blue squares.The markers with white triangles indicate the seismic stations.The green square shows the location of SAFOD.The markers with white triangles indicate the seismic stations.The arrows in the left panels show the sense of slip as the rightlateral.The lines in the elements show the orientations of the maximum contractional strain.The color contour shows the amplitude of cumulative strains, which is obtained by subtracting the initial state of strains on 3 January 2009, from the current time snapshot.Max shear can be negative as it is the relative value from the initial state.As we evaluate the strain field with the assumption of the continuum in the elements, the localized slip can cause bias in the estimation of strain.

Figure 14 .
Figure14.Comparison of the cumulative strains associated with the (a) dilation and (b) max shear to the dv/v with the frequency band of 0.9-1.2Hz.We used the dv/v of the seven stations listed in the legend, while we excluded the GHIB, JCSB, and VARB due to the low quality of the dv/v measurements.We compared the strain and dv/ v at the synchronized time periods every 90 days.The reference slopes are indicated for the sensitivity of dv/v to the strain.The slope and standard error of the linear regressions using all seven stations are annotated in the panels.

Table 2
Model Parameters and the Ranges Used for the Markov Chain Monte Carlo Sampling a The s1 is constrained up to half of s2 for each iteration.