Diurnal and Semidiurnal P‐ and S‐Wave Velocity Changes Measured Using an Airgun Source

Monitoring the subtle seismic velocity changes provides a promising tool to study the Earth's dynamic processes. However, the precise monitoring of seismic velocity with high temporal resolution is still challenging especially at large distances. Here, we continuously monitor P‐wave and S‐wave velocities with a precision of up to 10−5 using an airgun source. During a 1‐week experiment conducted in the Yunnan Province, Southwestern China, unambiguous diurnal and semidiurnal velocity changes with amplitudes of ~10−4–10−3 were observed for both P‐ and S‐waves at distances up to 21 km. The amplitudes of P‐wave velocity (Vp) change are 1.2 to 3.8 times the amplitudes of S‐wave velocity (Vs) changes. The velocity variation decreases with the increase in distance and is larger in the basin region than that in the mountainous region. Seismic velocity change shows less correlation with volumetric tide strain but strong correlation with changes of air temperature and air pressure. We propose that the thermal strains from air temperature changes are the primary cause of the observed diurnal and semidiurnal seismic velocity changes. The stronger Vp than Vs change is explained by the higher Vp sensitivity to the saturation rate. The velocity changes in mountainous area are ascribed mainly to the crack density changes, while both crack density change and saturation change are required to account for the velocity changes in the basin.


Introduction
The seismic velocities of rocks change with the applied stress owing to the opening/closing of microcracks or pores, which were well documented in laboratory (e.g., Nur & Simmons, 1969) and field observations (e.g., Larose et al., 2015). Monitoring seismic velocities provide a promising approach for estimating the temporal variations and promoting our understanding on the dynamic evolution of the subsurface (Froment et al., 2013;Mordret et al., 2016;Niu et al., 2008). The seismic velocity is observed to change on different temporal scales. Short-term seismic velocity changes are associated with hazardous events (Brenguier et al., 2008;Liu et al., 2014;Niu et al., 2008;Olivier et al., 2015;Taira et al., 2015). Long-term (seasonal) seismic velocity changes are generally attributed to changes in the air temperature (Meier et al., 2010) and hydrologic loading (Mordret et al., 2016;Tsai, 2011).
In addition to the responses to strong external loading, such as earthquakes and seasonal meteorological loadings, the seismic velocity also changes due to daily loadings. Weak diurnal and semidiurnal seismic velocity changes have been well observed via active source measurements (e.g., Reasenberg & Aki, 1974;Silver et al., 2007;Wang et al., 2008). Recently, diurnal, semidiurnal, and terdiurnal components in the spectrum of seismic velocity changes have been reported via ambient noise interferometry (Mao et al., 2019). Diurnal and subdiurnal velocity changes may be caused by Earth tides (e.g., Reasenberg & Aki, 1974;Yamamura et al., 2003), barometric loading (e.g., Silver et al., 2007;Wang et al., 2008), or a combination of tides and solar radiation (Mao et al., 2019). Monitoring the diurnal and subdiurnal seismic velocity changes also benefits the extraction of possible short-term precursory signals prior to earthquakes (e.g., Niu et al., 2008). However, the relation between the seismic velocity changes and external loadings is not completely understood due to limited high-precision and high-temporal resolution seismic velocity measurements.
The precision of seismic velocity measurements strongly relies on the repeatability of seismic sources (e.g., Cheng et al., 2010;Zhan et al., 2013). Various repeatable active sources have been used to monitor the seismic velocity changes (Maeda et al., 2015;Niu et al., 2008;Reasenberg & Aki, 1974;Yamamura et al., 2003), and most previous sources have focused on single phases at short distances (<km). Recently, large volume The airgun source (black star) was installed in a water reservoir (gray shadow). A source station (open red triangle) is deployed 50 m from the airgun to record the source time function. Clear diurnal and semidiurnal velocity changes are observed at some stations in basin (or fault zone) and mountain marked with solid red and blue triangles, respectively, and the size of solid triangles are scaled to the amplitude of the S-wave velocity changes. The open triangles indicate seismic stations without observed diurnal or semidiurnal velocity changes. Local faults are denoted by black lines. The dashed lines present the approximate basin boundary, which is manually extracted from the elevation model. airgun arrays are used to explore the crustal velocity structure (Chen et al., 2007(Chen et al., , 2017. The airgun source is highly repeatable with strong P and S phases propagating to long (up to hundreds of kilometers) distances , which provides an opportunity to monitor both Pand S-wave velocities. The small volume airgun has been used to monitor the P-wave velocity changes along short baselines (e.g., Reasenberg & Aki, 1974). In this paper, we report the diurnal and semidiurnal Vp and Vs changes at large distances (>20 km) observed from a continuous active source experiment using an airgun source.

Field Experiment
In 2011, a permanent seismic source was built in a water reservoir in Binchuan County, Yunnan Province, Southwestern China, to monitor the subtle seismic velocity changes in the Red River fault zone and the proximity with the anticipations of better understanding the dynamical evolutions (Chen et al., 2017;Wang et al., 2012). The source is comprised an airgun array with four 2,000 in 3 LongLife airguns manufactured by Bolt Co. More details of the source were provided by Wang et al. (2012Wang et al. ( , 2018. The airgun source was put into routine operation since September 2012 (Xiang et al., 2019).
Around the airgun source, 40 portable three-component short period (0.5-50 Hz, 100 samples per second) seismic stations were deployed ( Figure 1). One of these stations (called the source station, hereafter) was installed~50 m from the airgun source ( Figure 1) to record the source signature. All stations were synchronized to the Global Positioning System (GPS) once in every 2 hr to calibrate their internal clocks.
To evaluate the precisions of the velocity changes measured with the airgun source and portable seismic stations, we performed a 163-hr continuous experiment from 6 p.m. on 2 November to 1 p.m. on 9 November 2012 (local time). During this experiment, we fired a total of 138 airgun shots at an interval of 1-2 hr. No precipitation was observed during the experiment.
The source repeatability is of key importance in measuring the subtle subsurface velocity changes (e.g., Zhan et al., 2013). The airgun signals are affected by several dominant factors, such as water level, towing depth, and working pressure (Chen et al., 2014). During the experiment, the water level was continuously monitored with a precision of 0.1 cm and exhibited no discernable change. In addition, we maintained the other influencing factors with the towing depth and working pressure fixed to 12 m and 15 MPa, respectively.

Source Characteristics
We first cropped the 20-s-long airgun signals starting from the origin times provided by the GPS synchronized airgun firing system. The seismic signals recorded by the source station ( Figure 2) present typical signatures of an airgun array (e.g., Caldwell & Dragoset, 2000;Wang et al., 2018). The source signal comprises two main parts, an impulsive higher-frequency (7-30 Hz) signal originated from the high-pressure air release and long-lasting (~1.5 s) low-frequency (2.5-5 Hz) signal from the bubble oscillation (e.g., Caldwell & Dragoset, 2000;Wang et al., 2018). The slowly decaying low-frequency bubble signals are used in this paper to achieve high signal-to-noise ratios (SNRs) at large distances (e.g., Chen et al., 2007;Wang et al., 2018).
To estimate the stability of the airgun source, which is of key importance in velocity monitoring, we first filtered the vertical component source signals to the frequency of 2-8 Hz a little broader than the band of the bubble oscillation. Then, we calculated the cross-correlation between each filtered trace and the linearly stacked reference trace recorded at the source station. Generated under the same conditions, the seismic signals from the different shots were nearly identical (Figures 3a and 3b) with cross-correlation coefficients >0.99. High cross-correlation coefficients reaffirm the high repeatability of the airgun source (Chen et al., 2007).  Figure 1). The source is characterized by two parts: a high-frequency (7-30 Hz) pressure pulse signal and a low-frequency (~4 Hz) bubble oscillation signal lasting for~1.5 s . For clarity, we shifted the waveform and spectrum 0.5 s to the right.

Retrieving the Green's Function
The low-frequency airgun bubble signal lasts for~1.5 s ( Figure 2); the signal duration may cause phase interferences and affect the velocity measurements. To suppress possible phase contaminations, we retrieved the impulsive response between the source and each far field station, which is known as the Green's function (GF), using a deconvolution technique. This technique has been used in teleseismic receiver function (Langston, 1979) and active source monitoring (e.g., Maeda et al., 2015;Yang et al., 2018). We adopted a water level deconvolution method in the frequency domain (Langston, 1979) and calculated the Fourier spectrum (G) of the GF as where R and S are the spectra of the far field record and the source signal, respectively, and "*" represents the complex conjugation. A water level constant c is used to ensure the stability of the division operation (Langston, 1979). The water level is usually empirically set from 1 to 0.1% according to the SNR of the source signal. The records from the source station are used as the source signals (as in Maeda et al., 2015;Yang et al., 2018). Benefiting from the high SNR (~10 6 ) of the source signals, we used a very low water level c of 0.01%. The time domain GF is then obtained through inverse Fourier transform of G. The source signatures are effectively removed from the GFs with a low water level (Figures 3, 4, and 5).
The deconvolution process effectively reduces the source duration from~1.5 s ( During the experiment, we kept all working conditions for aigun fixed. However, slight fluctuation, especially the working pressure, is inevitable. After the continuous experiment we fired four additional shots with a working pressure 13 MPa, 2 MPa lower than used in the continuous experiment; the pressure difference is much larger than the possible working pressure fluctuation in the continuous experiments. The source wavelets are clearly visible at both the source station (a) and the distant station (b) but disappear in (c) the GFs. The raw waveforms are subject to trigger errors resulting from the finite sampling rate, which are corrected in the GFs. For clarity, we shifted the waveforms in (a) and (b) 0.5 s to the right. Subtle source changes induced by the working pressure fluctuations are efficiently eliminated via deconvolution ( Figure 5).

Measuring Traveltime Changes
All 138 vertical GFs were then linearly stacked to provide a reference GF for each station. Both individual and stacked GFs were further filtered with a fourth-order two-pass 2.5-5 Hz Butterworth band-pass filter ( Figure 6). For simplicity, we only used the vertical component signals herein. To minimize the possible phase interferences, we focused on the first arrival and maximum amplitude ( Figure 6). The first arrival and peak amplitude travel with apparent velocities of~4.5 and~2.3 km/s (Figure 4), which are regarded as P-wave and S-wave, respectively. The phase identification is also supported by the particle motion analysis ( Figure S1 in the supporting information). The changes in traveltime relative to the reference GF for P-and S-waves were then calculated using the sliding window cross-correlation technique (e.g., Niu et al., 2003;Wang et al., 2008).
Subsample time delays are usually required to monitor subtle velocity changes, which can be achieved using refined cross-correlating techniques (e.g., Silver et al., 2007;Wang et al., 2008). The cross-correlation refinement can be fulfilled by interpolating the waveforms to a very high sampling rate prior to cross-correlating or by interpolating the CC functions from the original or slightly interpolated waveforms. Both methods perform equally well (e.g., Céspedes et al., 1995). However, the computation cost is higher for the first method than that for the second. Therefore, we interpolated the waveform to 500 samples per second (5 times the original sampling rate) with a cubic spline function and then measured the subsample time delays from the cosine curve-fitted correlation functions (e.g., Céspedes et al., 1995;Wang et al., 2008). We also estimated the traveltime delays by cross-correlating the 100-time (i.e., 0.1-ms sampling  interval) interpolated waveforms. The results from the two methods are nearly identical (Figure 6c).
We calculated the cross-correlation between each GF and the corresponding reference trace with a 0.6-s-long moving window and running step of one sampling interval (i.e., 0.01 s). The average and standard deviation of dt within 0.3-s windows around the P-and S-waves (boxes in Figures 4 and 6) were taken as the traveltime delay and its measurement error, respectively.

Results
The traveltime changes (dt) of the P-and S-waves at station 53034 ( Figure 1) are consistent and oscillate with amplitudes of 2 and 1 ms, respectively (Figures 7a and 7b). Benefiting from the high signal repeatability (Figures 6c and 7c), the standard deviations of dt for the P-and S-waves are~0.2 and 0.1 ms, respectively.
Both dt values are dominated by daily (24 hr) variations (Figures 7a and 7b). The diurnal and semidiurnal modes are even notable in the frequency domain (Figure 7d). To further emphasize the diurnal and semidiurnal traveltime variations, we filtered dt from 2 to 80 hr. The filtered dt preserves the major features of the raw dt (Figures 7a and 7b).
Following the above procedure, we measured traveltime changes (dt) for P-and S-waves at all stations ( Figure 4). The filtered dt at all stations with epicentral distances <10 km (except at station 53278, which had serious clock drift) and the station 53281 (with an epicentral distance of 21 km) changed primarily within 1 ms and presented similar diurnal and semidiurnal patterns in both time and frequency domains (Figures 8a and 8b). The traveltime changes are divided by the absolute traveltimes dt/t (Figure 4), inverse to the relative velocity changes dv/v, are of the order 10 −4 (Figure 9a). These amplitudes of dv/v are of an order similar to the previous observations from active (e.g., Reasenberg & Aki, 1974) and passive source monitoring (e.g., Mao et al., 2019).

Precision of the Traveltime Measurements
The seismometer clock plays a primary role in seismic velocity measurements and may be affected by factors as air temperature (Silver et al., 2007). During this experiment, all seismometers were synchronized to the GPS clock every 2 hr. The clock drifts for all the synchronizations are no larger than 0.1 ms with more than 98% having drifts <10 μs (Figure 10), which is 2-3 orders smaller than the observed traveltime changes (several milliseconds; Figure 8a).
The continuous ground motion was digitized with a finite sampling interval, which is 0.01 s in our case. The waveforms trimmed to the origin time show subtle offsets (Figures 3a and 3b); these offsets primarily reflect the triggering error from the finite sampling interval (e.g., Niu et al., 2008;Silver et al., 2007). The triggering error can be corrected via the deconvolution process, and retrieved GFs (Figure 3c) show better alignment than the corresponding raw waveforms (Figure 3b) leaving only the traveltime changes of milliseconds ( Figure 8).
The precisions of the delay time measurements (ε) are subject to a theoretical boundary known as the Cramer-Rao lower bound ε CRLB (e.g., Céspedes et al., 1995): (2) where f 0 ,B,ρ,and T are the dominant frequency, fractional bandwidth, signal similarity, and window length used for the cross-correlation, respectively. For the airgun signal, we have f 0 ¼ 4Hz; B≅ 5−2 4 ¼ 0:75; ρ≅1; and T ¼ 0:6 s. Then the precision is dominated by the value of SNR. At the frequency of airgun signal (~4 Hz), seismic records are strongly affected culture activities, which show strong diurnal variation (McNamara & Buland, 2004). Most seismic signals used in this study have SNR above 200; then from equation (2) we obtain ε CRLB ≈ 1.5 × 10 −4 s, which agrees with the observed variation of~10 −4 s (Figures 7 and 9b) and with the boot-strap test (Sawazaki et al., 2018) result ( Figure S2). Since the P-and S-wave traveltimes are several seconds (Figure 4), the precisions of the relative traveltime or velocity measurements are at the order of 10 −5 .

Spatial Variation of Velocity Changes
All velocity changes show nearly in-phase diurnal and semidiurnal tendency with different amplitudes (Figures 8a and 9a). The amplitude of the relative velocity change (dv/v) decays rapidly from 10 −3 (P-wave) and 10 −4 (S-wave) at station 53277 (with 750-m distance) to~10 −5 (both P-and S-waves) at station 53274 (at 8.8-km distance) (Figures 1 and 9c). Under low stress perturbation, the seismic velocity change can be expressed as the product of the stress change Δσ and the velocity sensitivity η (¼ ∂v v ∂σ ). The velocity sensitivity has long been observed to decrease with the effective pressure or burying depth (e.g., Nur & Simmons, 1969). The stress changes may also decay exponentially with depth (Ben-Zion & Leary, 1986;Tsai, 2011). The observed distance-dependent velocity changes (Figure 9c) may reflect the depth decay of both η and Δσ since the ray path pierces deeper at larger propagation distances.
Note that the velocity changes at basin stations are stronger than velocity changes at mountain stations (Figures 9a and 9c). With lower stiffness, sedimentary basins are more susceptible to external loading than bedrock. Different material properties can also result in different stress changes in sediments and bedrock. Both factors contribute to the different velocity changes in the basin and mountainous regions.

Possible Mechanisms for Diurnal and Semidiurnal Velocity Changes
Diurnal and semidiurnal velocity changes are usually attributed to changes in solid Earth tides (Mao et al., 2019;Reasenberg & Aki, 1974), ocean tides (Yamamura et al., 2003), barometric pressure (Silver et al., 2007;Wang et al., 2008), and air temperature (Mao et al., 2019;Sens-Schönfelder & Larose, 2008). Herein, we compared the average velocity change (Figure 9b) with the tides, air temperature, and air pressure at different periods (Figures 11 and 12). The volumetric tidal strain was calculated with the Earth tide data processing package ETERNA (Wenzel, 1996), where we did not take the oceanic tides into consideration because our experiment site is more than 800 km from the nearest coastline (Figure 1 inset). The air pressure and temperature were measured at the source site with a sampling rate of twice an hour.

Solid Earth Tide
The strongest constituent of the solid Earth tide is the M2 tide (period 12.421 hr), which is believed to govern the semidiurnal seismic velocity changes (Mao et al., 2019). However, our velocity change correlate well to the tide in 1-day period (cc = −0.90; Figure 11b) but poorly correlates to the stronger half-day (cc=−0.21; Figure 11c) and broadband (cc = −0.49; Figure 11a) tides. Furthermore, Earth tides are believed to cause instant velocity changes (Mao et al., 2019;Reasenberg & Aki, 1974), which are not observed here. Therefore, the solid Earth tide is unlike the primary cause of the observed velocity changes.

Barometric Pressure
The barometric pressure variation may cause velocity changes by altering the bulk and pore pressures (Silver et al., 2007;Wang et al., 2008). And every kPa air pressure change will induce 10 −9 elastic strain (e.g., Rabbel & Zschau, 1985). The observed 0.5-kPa air pressure change will result in 5 × 10 −10 elastic strain at shallow part. However, the elastic response to the air pressure is believed to be instant (e.g., Hillers et al., 2015;Silver et al., 2007); the observed velocity changes and the barometric pressure changes are out-of-phase (Figure 11).

Journal of Geophysical Research: Solid Earth
Delayed velocity changes may occur in response to the poroelastic strains from changes in the underground water table (Tsai, 2011) and likely the air pressure. To explain the observed velocity changes with amplitude of 10 −3 -10 −4 , we would need several meters of water table change or tens of kilopascals of air pressure change (Tsai, 2011). Though we are lacking water table data, a daily water table change of several meters does not seem reasonable. And the observed air pressure change (hundreds of Pascal) is about 2 orders below the required level of tens of kilopascals.

Air Temperature
The air temperature can change the strain at depth via thermal expansion and the transmission of stress (Ben-Zion & Leary, 1986;Berger, 1975;Tsai, 2011). The strain changes (Ben-Zion & Leary, 1986;Berger, 1975), as well as the velocity changes (Tsai, 2011), at a certain depth can be predicted in the form of delayed, damped, and low-pass filtered air temperature variations. Regarding the air temperature as a standing wave with amplitude T 0 , frequency ω,wavelength λ, and wavenumber k = 2π/λ, the amplitude of the thermal strain A can be represented as a function of the linear thermal expansion α th , thermal diffusivity κ, incompetent layer thickness y b , and Poisson's ratio ν as (Tsai, 2011): The thickness of the incompetent layer y b determines the lag t between the thermal strain and the temperature Δt= y b 2 ffiffiffiffi 2 ωκ q (Ben-Zion & Leary, 1986).
The diurnal and semidiurnal velocity changes lag by 22 and 11 hr behind the air temperature with possible one-day and half-day cycle skipping, respectively (Figures 11b and 11c). Taking the representative value κ = 10 −6 m 2 /s, we estimated incompetent layer thicknesses of 0.95 and 0.67 m for the diurnal and semidiurnal changes, respectively. These thicknesses agree with previous estimations of~0.6 m by Ben-Zion and Leary (1986) and Tsai (2011).

Journal of Geophysical Research: Solid Earth
Berger's model is sufficient to explain the observed 10 −4 seismic velocity change (e.g., Richter et al., 2014) but only with a~3-hr delay. The model with an incompetent layer can predict the observed lag time but insufficient strain amplitude.
The semidiurnal air temperature change (~2°C) is about~1/6 of the diurnal temperature variation (~12°C), and the semidiurnal thermal strain is expected no larger than 1/6 of the diurnal thermal strain. And the thermal strain-induced semidiurnal velocity change is expected to be no larger than 1/6 of the diurnal velocity change, while the observed semidiurnal velocity change is~1/2 of the diurnal velocity change (Figures 11  and 12). Thus, the amplitude ratio between the semidiurnal and diurnal velocity change is hard to be explained by the simple thermal strain model. Richter et al. (2014) argued that material heterogeneities and topographic effects contribute to the time lag. These effects may also influence the amplitude of the thermal strain (Ben-Zion & Leary,1986) as well as the spatial variation of velocity change (section 5.2).
In summary, though we take the thermal strain as the most like cause of the velocity change, all possible factors may contribute to the observed velocity change, which cannot be perfectly explained by any single model. In addition, the measured traveltime changes are integrations along corresponding ray paths; thus, the depth-dependent velocity changes should be taken into consideration for a comprehensive understanding of the mechanism. A sophisticated model incorporating these effects, fine basin structure (Xu et al., 2018), and longer observations will help decipher the mechanism of these velocity changes.

Different Pand S-Wave Velocity Changes
Limited by the source characters, most previous studies focused on either P-wave (Reasenberg & Aki, 1974;Yamamura et al., 2003) or S-wave velocity changes (Mao et al., 2019;Meier et al., 2010). The airgun source used in this study can efficiently generate both P-and S-waves , which make it feasible for us to spontaneously monitor the P-wave velocity (Vp) and S-wave velocity (Vs) changes.
The Vp changes are stronger than the corresponding Vs changes at all stations ( Figure 9c). This observation is consistent with the results obtained from laboratory experiments, indicating that Vs is less sensitive to the effective stress than Vp, especially at shallow depths (low effective pressures; e.g., Nur & Simmons, 1969;Winkler & Liu, 2005). Amplitude difference in Vp and Vs changes may contribute to the observed different coseismic Vp (e.g., Pei et al., 2019) and Vs (e.g., Cheng et al., 2010) changes of the Wenchuan earthquake.
The seismic velocity changes of cracked rocks are generally attributed to the opening/closure of microcracks or pores. We estimate the P-and S-wave velocities ( Figure 13a) as functions of crack density (CD, number of cracks in certain volume) and saturation fraction (SF) using a crack model (O'Connell & Budiansky, 1974) based on self-consistent approximation. Both Vp and Vs increase with the increase of SF, and both decrease with the increase of CD. At low CD (CD<0.4) or high SF, Vp is more sensitive to SF than Vs as widely observed in laboratory (e.g., Nur & Simmons, 1969).
To reveal the microscale origin of our observed velocity changes (Figure 9a), we project the observations to the predicted Vs-Vp space (Figure 13s). The predicted velocities in Figure 13 are normalized to the velocities of intact (without crack) rocks. Taking the typical crust velocities Vp=6.5 km/s and Vs=4 km/s as intact velocities, we normalize the observed velocity changes (Figure 9a) to 0.6, in the Vs-Vp space (Figure 13b) for both P-and S-waves, which propagate at apparent velocities of 4.5 and 2.3 km/s, respectively (Figure 4). The observations cluster into two groups with different slopes. The first group (red points in Figure 13b) lies oblique to isolines of SF and CD. The slope of the first group is larger than the slope of any isoline of SF; it means the CD change itself is insufficient to explain the observed velocity changes of the first group. The second group (blue points in Figure 13b) distributes close to constant SF; the contour lines of SF are parallel with similar slope in the whole Vs-Vp space (Figure 13a). Therefore, velocity changes of the second group are mainly result from the stress-induced CD changes.
The different slopes observed in Figure 13b may originate from the difference in rock property along each ray path. The stations of the first group (red triangles in Figure 1) are located in sedimentary area and in the vicinity of a local fault (Figure 1), which is generally characterized with high crack/pore density (e.g., Mitra & Ismat, 2001) and high pore fluid pressure. The stations of the second group (blue triangles in Figure 1) distribute in the mountainous area, where igneous rocks are wide spreading (Luo et al., 2015), and are expecting to have lower CD and lower pore saturation than fault zone or basin.
Subjecting to external loadings, both CD and saturation ratio will change. The P-wave velocity of porous rock is very sensitive to the SF, especially approaching full saturation, as observed in laboratory (e.g., Naesgaard et al., 2007) and predicted with crack model (Figure 13a). The stress-induced SF changes in less porous mountain are likely weaker than in porous basin. Therefore, velocity changes in the mountain area (blue points in Figure 13b) are dominated by the CD changes and are weaker than the velocity changes in the basin (red points in Figure 13b), where both SF and CD contribute.

Conclusions
In this study, we measured seismic velocity changes using an active source with high precision (10 −5 ) and high-temporal resolution (~1 hr) at distances up to 21 km, which is longer than most previous studies. Airgun signals can be tracked to hundreds of kilometers with less than 1-day stacking (Chen et al., 2017). The long propagating distances and rich phases ) of airgun signals make it possible to monitor the crust with high temporal resolution and high precision (e.g., Leary & Malin, 1982). The airgun source has been operated for more than 7 years; more results on seismic velocity changes at different scales will be presented in the follow-up studies (e.g., Xiang et al., 2019).
Unambiguous diurnal and semidiurnal velocity changes were observed for both P-and S-waves, with the Pwave velocity changes being larger than the S-wave velocity changes. The amplitude of the velocity change decayed with the epicentral distance. The velocity changes in the sedimentary basin were stronger than velocity changes in the crystallized mountainous area. The velocity correlates better with the air temperature than with Earth tides; therefore, we suggest that the thermal strain is the main cause of diurnal and semidiurnal velocity changes. Different microstructural responses of basin and mountain materials to the external loadings explain the different velocity changes in basin and in mountain. This study emphasizes the importance of the synchronous Vp and Vs monitoring, which helps the understanding the physical mechanism of the seismic velocity change.  (Figure 1), respectively. The velocity changes from basin (or fault zone) and mountain stations have different slopes; this difference is even remarkable between red circles (stations 53264 and 53034) and blue circles (station 53274).