The M _S 6.4 Yangbi earthquake is located in western Yunnan, which is located in the southeast margin of Tibet Plateau (Figure 1). It is the Yangtze paraplatform and northwest Yunnan geosynclinals fold belt junction area, which has strong characteristics of structural tension since the Quaternary period (Huang et al., 2014). In the present study, the vertical component of continuous waveforms (100 Hz) that were recorded in six broadband stations that are located within approximately 100 km of the epicenter of the Yangbi Earthquake between January 2019 and December 2021 was utilized (Figure 1). The area hosting the hypocenter was adequately covered by ten station pairs, and the minimum, maximum, and average distances between these stations pairs are 52, 143, and 97 km, respectively.

The data preprocessing procedure that was utilized in the present study was described in Liu et al. (2014). Raw seismic data were partitioned into 1-day intervals and the vertical component data were then resampled at 10 Hz to improve computational efficiency. Temporal normalization and spectrum whitening in the 1–20 s range were applied to each segment to minimize earthquakes effects. Reference cross-correlation functions (CCFs) for station pairs were obtained via the stacking of CCFs covering the period from 1 January 2019, to 21 April 2021 (a month preceding the Yangbi Earthquake). To improve the signal-to-noise ratio of daily CCFs, these were derived through the stacking of CCFs for 61 d, which included 30 d before and after a target day. To prevent mixing of pre- and post-seismic signals, daily CCFs that were calculated for periods before and after the Yangbi Earthquake were stacked separately (Liu et al., 2014; Liu et al., 2022). For example, the daily CCF on 28 May 2021 is obtained from stacking only post-seismic daily CCFs from 22 May 2021 to 28 June 2021. Thus, the stacking days of daily CCFs is smaller than 61 d within a month before and after the main shock. Daily CCFs for the EYA–YUL station pair in the period band of 1–10 s from January 2019 to December 2021 are shown in Figure 2A, and clear surface and stable coda wave signals are obvious. Owing to the decrease of the coda wave coherence as the timelapse increases, the lapse window for the positive portion of CCFs was determined as 30–130 s after the arrival of Rayleigh waves in this period band (Supplementary Figure S1). The windows in the negative portion were symmetrical to that in the positive portion.

To supplement station pairs near the epicenter, the single-station autocorrelation approach was also considered. Autocorrelation functions mainly reflect changes in the shallow crust near a station, and the associated flow processing is similar to that for the cross-correlation of a station pair. Figure 2B shows the autocorrelation functions for the EYA station using 1–3 s band-pass filter from January 2019 to December 2021. The lapse windows for the autocorrelation were determined as the fixed windows with range of ±(5–55) s.

Large to medium velocity perturbations can be directly obtained by measuring the relative traveltimes of the direct waves. Studies on the detection of small–medium changes as a function of time have focused on seismic coda waves. Coda waves, which are also known as multiply scattered waves, usually arrive later than direct waves. The later the arrival of scattered wave phases at the receiver, the longer the associated propagation paths and the higher the sensitivity to minor perturbations in a medium (Snieder et al., 2002).

Assuming that the change of relative seismic wave velocity (δv/v) is spatially homogeneous, the travel time shift δt between CCFs is proportional to the lapse time t and can be calculated as δv/v = −δt/t (Poupinet et al., 1984; Snieder et al., 2002). Therefore, the measurement of the travel time shift δt is very important for coda wave interferometry. In the present study, the MWCS method was utilized to measure the relative time shift between a reference CCFs that corresponds to the initial state, and a current CCFs that has encountered a velocity change in the medium. With the MWCS method, a series of overlapping time windows are defined in the coda wave, and the time shifts in these windows are estimated by means of the cross-spectrum method. The cross-spectrum X(f) between the reference CCFs and current CCFs is calculated as follows (Clarke et al., 2011): X f = F r e f f ⋅ F c u r * f (1) where F r e f and F c u r are the Fourier-transformed segments of the reference and current CCFs. The asterisk denotes the complex conjugation and f is the frequency. X f can also be expressed by its amplitude X f and phase e i ϕ f . If the time-shift is constant in each window segment, ϕ is linearly proportional to f: ϕ f = 2 π δ t f (2)

The time shift for each window segment is the slope and the associated error e m is caulcuated as e m = ∑ j w j v j ∑ i w i v i 2 2 ∑ j ϕ j − m v j 2 N − 1 (3) where w are weights, m is the slope of ϕ , and v are 2 π f .

After all time shifts δt are measured in the window segments in an interest range of lapse time, the relative time shift δt/t is estimated by a weighted linear regression passing through zero, and then δv/v can be obtained by −δt/t (Poupinet et al., 1984; Clarke et al., 2011). To evaluate the reliability and accuracy of the method, waveform modeling data that were reported in Yuan et al. (2021) were used, and the +.1% velocity perturbation that was determined in the model was correctly measured (Supplementary Figure S2).

Figure 3 shows an example of the relative time shift estimation using CCFs for the EYA–TUS station pair. Compared with the theoretical waveform in Supplementary Figure S2, CCFs of the EYA–TUS station pair produced a lower coherence and preserved the acausal signal. Regarding a 1–10 s period band and a station pair with distance of 62 km, the timelapse window is between 53 and 153 s and window segments are 18 s with a step of 1 s. Time shifts ( δ t ) between reference and daily CCFs can be measured in each window segments by using the MWCS, and the slope ( δ t / t ) can then be estimated via a weighted linear regression of the time shifts.

Variations in spatial distributions of noise sources and the instrumental clock errors can independently affect the measurement of the travel time shift (Stehly et al., 2007). Clock errors can produce an overall shift in the entire cross-correlation time, thereby increasing traveltimes in the positive portion and decreasing traveltimes in the negative portion, and vice versa. Comparatively to the measurement of velocity changes in a medium, a direct arrival surface wave was used instead of a coda wave to measure instrumental clock errors. Clock errors that was measured from the surface wave using the noise cross-correlation technique can be expressed as follows (Stehly et al., 2007): δ τ t + δ τ − t 2 = D t + ε t + ε − t 2 (4) where δ τ t denotes the variation in the traveltime of the surface wave that is measured in the positive or negative portion. D t is the time shift caused by instrumental clock errors, and ε t is the time shift associated with the spatial variation of noise sources. Therefore, clock errors can be estimated using Eq. 4 by assuming that D t is greater than ε t + ε − t 2 .

Following approaches that were advanced in Section 2.2 and Section 2.3, only the coda wave window was replaced with a surface wave window, that is, it starts 30 s before the surface wave time and ends 40 s after. Supplementary Figure S3 shows an example of clock errors that were estimated using data from the LUS–YUL station pair. The clock errors that were obtained on August 1 and October 20, for example, are .42 and −.04 s, respectively. The clock error on October 20 is less than one sampling rate, and thus, it can be considered as zero. Clock errors for these stations were evaluated for the period from 2019 to 2021, and results for station pairs with possible large clock errors are shown in Figure 7. The EYA station displays a clock drift of approximately −0.2 s throughout 2019, whereas LUS station exhibits a clock drift of −.6 s between July and August 2021. Considering that daily CCFs were obtained by stacking CCFs of 61 days, the estimated clock drifts are likely underestimated, in particular, the clock drifts from the LUS station may reach seconds.

Two groups of teleseismic earthquakes that occurred in the Alaska Peninsula and the Philippines are selected to validate the clock drifts (Supplementary Tables S1, S2; Supplementary Figure S4). Differences in traveltimes of the phases between stations should essentially stable over time for nearby teleseismic earthquakes. But the reference arrivals of earthquakes on July 29 and 14 August 2021, were relatively early at the LUS station compared to that of other three earthquakes (Supplementary Figure S5). If the linear trend in Supplementary Figure S5B is eliminated, more intuitive reference arrivals relative to the distance from the epicenter can be obtained. Supplementary Figure S6A demonstrates that arrivals of earthquakes on July 29 and 14 August 2021, at the LUS station significantly differ from those of the other three earthquakes, and the drift is −1 s. Arrivals of all earthquakes at the other four stations do not show any obvious drift between 22 July 2020, and 11 October 2021. In addition, arrivals of earthquakes at the EYA station in 2019 slightly differ from those of the other three events, and the drift is −.3 s (Supplementary Figure S6B).

Considering the LUS–YUL station pair as an example, the velocity change that was calculated using the least squares fitting MWCS method is less than .02% for a clock error of 1 s (Figure 8). This minimal impact is probably because the slope of dt/t based on the MWCS method is unaffected by such an overall time drift of the cross-correlation time. However, a large clock error reduces the correlation between the reference and daily CCFs, and this affects subsequent calculations.

Positive and negative velocities linked to earthquakes may correspond to regions of increased and decreased stress, respectively (Rubinstein and Beroza, 2004). To explain observations from the noise cross-correlation, static stress changes caused by the earthquake were calculated using Coulomb 3.0 (Lin and Stein, 2004). The static slip distribution of the Yangbi Earthquake provided by Xu Zhang was utilized (https://www.cea.igp.ac.cn/kydt/278248.html, see Supplementary Figure S7), whereas the Poisson’s ratio and shear modulus were set to .25 and 32 GPa, respectively. The calculated volumetric strains caused by fault slips at a depth of 5 km exhibit an obvious four-quadrant distribution (Figure 9). The EYA and HEQ stations are in the dilatation region, where a decrease in velocity is anticipated, in fact, the velocity changes for the EYA–HEQ station pair decreased by .02%. The other four stations fall within the contraction region, where a velocity increase is expected, in fact, velocity changes for the BAS–TUS, TUS–YUL, and LUS–TUS station pairs corresponding increased by .05%, .04%, and .08%. However, how to decide if two stations are in different strain regions. For example, the velocity changes for the EYA–TUS station pair increased by .06%, whereas that for the EYA-YUL pair decreased by .03%. Therefore, an analysis of the spatial sensitivity distribution of coda waves is necessary.

Based on numerical simulations of seismic waves in 2D and 3D heterogeneous elastic media, Obermann et al. (2013) suggested that the sensitivity of coda waves can be modeled as a linear combination of the sensitivities of body and surface waves. They indicated that early coda waves are dominated by fundamental surface waves, which mainly reflect shallow perturbations, whereas later coda waves are dominated by body waves. In the present study, the timelapse windows for CCFs were determined as 30–130 s after the arrivals of Rayleigh waves, and this contained sufficiently long signals. The sensitivity kernel that is expressed as follows can be used (Pacheco and Snieder, 2005): K S , R , r 0 , t = 1 p S , R , t ∫ 0 t p S , r 0 , t ′ p r 0 , R , t − t ′ d t ′ (5) where S and R are the positions of the source and receiver, respectively; r₀ is the position of the local velocity variation; t is the center timelapse for doublet analysis; and p (s, r, t) is the probability that the wave has traveled from s to r during t. This probability can be approximated using the full-space solution of the diffusion equation, which is expressed as follows (Mao et al., 2019): p s , r , t = 1 ( 4 π D t ) 3 / 2 e − ‖ s − r ‖ 2 4 D t (6) where D is the diffusion constant. The multiply scattered waves that propagate in 3 dimensions can be described by D: D = C E ∗ l ∗ 3 (7) where C_E is the energy velocity and l* is the scattering mean free path. Considering that S waves account for most of the energy in coda waves, a ratio of 9:1 was used to calculate C_E: 1 C E = 0.89 V s + 0.11 V p (Obermann et al., 2016). Therefore, the key to determine the sensitivity kernel is to estimate the scattering mean free path l*. In general, at a larger scale (e.g., crustal) the mean free path is fairly constant relative to the frequency, but this may not be applicable at local scales. Chaput et al. (2015) estimated that the scattering mean free path for the Erebus volcano in Antarctica at 1.5 Hz is ∼2 km, and values slowly decreased as the frequency increased. Data for the scattering mean free path for the Yangbi area are scant, but theoretically values that involve 5%–10% heterogeneity are in the range of 2–10 km (Obermann et al., 2013). Here, l* = 5 km was considered in the period band of 1–10 s for analysis. If the distance between two stations is 60 km and the center time of coda wave window is 100 s, the normalized depth sensitivity of the scattering waves can be obtained based on Eq. 5 (Figure 10). The sensitivity of scattered waves is high near the source and receiving points and relatively low in the middle portion. It decays with increasing depth, and at ∼12 km it reduces to 10% of the value at surface. The depth sensitivity for station pair LUS–TUS with path crossing through the fault zone shows that the epicenter of the Yangbi earthquake (marked by the red star in Supplementary Figure S8) is located in the weak sensitivity zone (smaller than 10%). The scatter waves for the existing station pairs may not be able to sample the major rupture area of the Yangbi earthquake. Velocity changes obtained using the noise cross-correlation technique can be considered weighting effects along propagation paths, and the largest weight coefficient is obtained near the related two stations.

Considering that the static stress rapidly decays as the distance from the epicenter increases, if two stations of a station pair are in different stress-loading regions, the static stress of the station closer to the epicenter exerts a greater impact on the velocity change. For example, the velocity changes for the EYA–TUS, EYA–YUL, BAS–EYA, and HEQ–YUL station pairs were attributed mainly to static stress of the stations closer to the epicenter. According to the meteorological observation data of the Dali Center of the China Earthquake Science Experiment Field, there is no obvious change in the rainfall before and after the main shock (Su et al., 2022). These results suggest that velocity changes are likely linked to a combination of near-surface physical damage and static stress changes. This also explains the inconspicuous drop in the coseismic velocity that was obtained using the autocorrelation for the TUS station, which is closest to the epicenter, compared to those of other nearby stations.