Impacts of Reservoir Water Level Fluctuation on Measuring Seasonal Seismic Travel Time Changes in the Binchuan Basin, Yunnan, China

An airgun source in a water reservoir has been developed in the past decade as a green active source that had been proven effective to derive short-term subsurface structural changes. However, seasonal water level fluctuation in the reservoir affects the airgun signal, and thus whether the airgun signals can be used to derive robust seasonal variation in subsurface structure remains unclear. We use the airgun data observed in the Binchuan basin to estimate the seasonal variation of seismic travel time and compare the results with those derived from ambient noise data in the same frequency band. Our main observation is that seasonal change δt/t from airgun is negatively correlated to the variation of dominant frequency and water table fluctuation in the reservoir. One possible explanation is that water table fluctuation in the reservoir affects the dominant frequency of the airgun signal and causes significant phase shift. We also compute the travel time changes in P-wave from the empirical Green’s function after deconvolving the waveforms from a reference station that is 50 m from the airgun source. The dominant frequency after deconvolution still shows seasonal variation and correlates inversely to the travel time changes, suggesting that deconvolution cannot completely eliminate the source effect on travel time changes. We also use ambient noise cross-correlation to retrieve coda waves and then derive travel time changes in monthly stacked cross-correlations relative to a yearly average cross-correlation. We observe that seismic travel time increases to its local maximum in the end of August. The travel time changes lag behind the precipitation for about one month. We apply a poroelastic physical model to explain seismic travel time changes and find that a combined effect from precipitation and evaporation might induce the seasonal changes as shown in the ambient noise data. However, the pattern of travel time changes from the airgun differs from that from ambient noise, reflecting the strong effects of airgun source property changes. Therefore, we should be cautious to derive long-term subsurface structural variation from the airgun source and put more attention on stabilizing the dominant frequency of each excitation in the future experiments.


Introduction
Extensive field observations and lab experiments [1][2][3][4][5] have suggested that slowly evolving, low-amplitudes phases usually appear before the onset of natural hazards such as earthquakes, landslides, and volcanic eruptions. A better understanding of the related crustal deformation and responsible mechanisms of these phases are critical for effective early warnings of such hazards. Crust deforms under the loading of tectonic and climatological forces. In order to precisely measure the deformation from tectonic loading, we need Remote Sens. 2021, 13, 2421 2 of 21 to remove deformation related to climatological forces. Seismic velocity changes are usually used as a proxy to quantitatively understand crustal deformation. Temporal changes in seismic velocity can reflect fault zone co-seismic damage and post-seismic healing [6][7][8], volcanic eruption [6,9], groundwater level changes [10][11][12], temperature and atmospheric pressure variations [13][14][15][16][17][18], and solid earth and oceanic tidal deformation [17,[19][20][21]. Thus, high temporal resolution monitoring and high precision measurements of seismic velocity changes are necessary.
Seismic velocity changes have been estimated through measuring travel time or phase difference from active sources including explosion [21][22][23][24][25], electronic hammer [26], airguns [18,27,28], repeating earthquakes [29][30][31][32][33][34], and dephasing of ambient noise crosscorrelations (CCs) [6,35]. Among these methods, temporal changes in ambient noise amplitude and heterogeneous distribution of noise sources [12,17] may introduce a bias to the measurements of seismic velocity changes. The uncertainties in locations of repeating earthquakes might cause concern to the reliability of the measurements of seismic velocity changes. Explosions are not environmentally friendly and thus are not possible to operate in a populated region. In contrast, a large volume airgun array with stable signal excitation has been suggested as a better method to monitor seismic velocity changes [36].
Since 2011, a permanent seismic source [36][37][38] has been deployed in a water reservoir called Dayindian in Binchuan, Yunnan, China, to image velocity structure [39,40] and monitor seismic velocity changes [8,[41][42][43][44]. It has been in routine operation since September 2012. The source is composed of four 2000 in 3 LongLife airguns manufactured by Bolt Co., which can ignite sufficient seismic signal penetrating down to a depth of 10-20 km and propagating further to tens of kilometers away [36,38]. The four airguns were towed 2 m below the four corners of a 7 m × 7 m steel frame hanging under a tower crane [37]. The source causes negligible effects to the wildlife and local communities in terms of ground shaking and noise pollution [37].
We target the Binchuan region because it is located at the eastern margin of the collision zone between the Indian and Eurasia plates, accumulating strong shear strain [45], accommodating millions of people, and experiencing many devastating earthquakes. The largest one was the M7.75 Yongsheng earthquake which occurred on the Chenghai fault in 1515 [46,47]. This region has recorded numerous earthquakes, including eight events with magnitude greater than 7, and 70 with magnitudes greater than 6 [48].
Since the deployment of the airgun array, [18] used one-week hourly-shooting signals to observe diurnal and semidiurnal P-and S-wave velocity changes. They proposed that the thermal strains from air temperature changes are the primary cause. The reservoir water level did not fluctuate significantly during the one-week period. However, we must consider the effect of water level fluctuation when investigating seasonal changes in seismic velocity, as it may change a few meters from winter to summer. The authors of [49] observed that the dominant frequency of airgun data correlates to the changes in water level in the reservoir, which may potentially affect measurements of travel time changes. In this study, we use airgun data in the frequency ranges of 2-6 Hz to compute variation of the dominant frequency and amplitude, derive seasonal changes in seismic velocity, and explore how water level fluctuation of the reservoir affects seismic velocity changes through investigating the correlation between seismic velocity changes and variations of the dominant frequency and amplitude. We also use ambient noise in the same frequency range to obtain seasonal changes in seismic velocity for coda waves and compare them with velocity changes for body waves.

Airgun Data
The Binchuan airgun array shot a few tests in 2011 and began to operate routinely after September 2012. On average, a few shots were made every day ( Figure S1), accumulating hundreds to thousands of excitations per year. From June to August, the water level in the Dayindian reservoir is too low to excite airgun shot because of irrigation for farming ( Figure S1). However, in 2016, water level fluctuation in the reservoir did not affect the Remote Sens. 2021, 13, 2421 3 of 21 airgun operation and the average number of daily excitations is also higher than that in other years. We therefore use the airgun data in 2016 to investigate the seasonal variation of seismic arrival time of P-and S-waves.

Source Characteristics
Forty portable three-component short-period stations have been deployed since the operation of the airgun source, and we selected 14 stations with clear body wave arrivals to compute travel time change (Figure 1). Each station is equipped with a Güralp 40 T sensor (a flat response from 0.5-100 Hz) and a RefTek 130 data logger. We convert daily data from MiniSEED to sac and merge them if there is more than one segment in a day. We then extract all airgun signals from the original dataset, starting from the shooting time to 20 s after that.
We first inspect the waveform recorded at the nearest station CKT, which is~50 m from the airgun source. A near-field airgun wavelet is composed of primary pulse and bubble pulse ( Figure 2). Primary pulse is generated by high-pressure gas releasing into the water and is usually in high energy with a wide bandwidth, e.g., 7 to 30 Hz ( Figure 2). It was mostly used in shallow oil and gas exploration. In comparison, bubble pulse is generated by the oscillation of bubbles after high-pressure gas was released into the water. It is in a relatively low-frequency band from 2 to 6 Hz [36,38], so it was often used in deep seismic exploration ( Figure 2).
We compute power spectral density of the bubble pulse and define the frequency with the maximum power spectral density as the dominant frequency. We observe that the dominant frequency of the bubble pulse decreases in the summer at the reference station ( Figure 3a) and others (Figures 3b,c and S2), similar to the reports in [49]. In fact, the towing depth of the airgun array decreases in summer when water in the reservoir was used for farming. An experimental study from [50] also indicates that the dominant frequency decreases when the towing depth of the airgun array decreases. Such variation in frequency content is quite common. For instance, at a global scale, seasonal or multi-decadal variation of storm activities have been reported to cause long-term variation in noise frequency content [51,52]. It has been suggested that temporal variation of frequency content can cause apparent seismic velocity change [53]. In addition to the dominant frequency, we also observe that the maximum amplitude in the bubble pulse window decreases in the summer, which is similar to the reports in [50].   Original near-field waveform (a) recorded by a reference station CKT. The signal from 0 to 0.3 s is the primary pulse (gray area marked as primary). The signal from 0.6 to 2 s is bubble pulse. (b) The spectrogram of the waveform in (A). The primary pulse is in high energy and in a wide bandwidth, e.g., 7-30 Hz. The bubble pulse is in a relatively low frequency bandwidth, 2-6 Hz.

Airgun Data Quality
We interpolate the airgun signals to 0.001 s/sample and bandpass all waveforms with a 2-6 Hz filter, the dominant frequency band of the bubble pulse [38]. With the single shot data, we can observe clear P-wave and S-wave arrivals up to 30 km with apparent velocities of 5.5 km/s and 2.8 km/s, respectively (Figure 4), consistent with previous observations [38].
When inspecting the waveforms throughout the entire year, we observe strong clock drifts in the reference station CKT and other receivers, approximately from day 230 in 2016 (Figure 5a,c and Figure S3). Because the differential times are nearly the same among all stations ( Figure S4), it is mostly attributed to the timing log system of the airgun array. In additional to that, sensors in the 53273 and 53277 also have timing problem. To fix the clock drift, we determine P-wave arrival time at each station through a polarization analysis introduced by [54]. For a three-component airgun waveform in the passband of 2-30 Hz at station 53278 (Figure 6a), we first convert each component to an analytic signal.

Airgun Data Quality
We interpolate the airgun signals to 0.001 s/sample and bandpass all waveforms with a 2-6 Hz filter, the dominant frequency band of the bubble pulse [38]. With the single shot data, we can observe clear P-wave and S-wave arrivals up to 30 km with apparent velocities of 5.5 km/s and 2.8 km/s, respectively (Figure 4), consistent with previous observations [38].
When inspecting the waveforms throughout the entire year, we observe strong clock drifts in the reference station CKT and other receivers, approximately from day 230 in 2016 (Figure 5a,c and Figure S3). Because the differential times are nearly the same among all stations ( Figure S4), it is mostly attributed to the timing log system of the airgun array. In additional to that, sensors in the 53273 and 53277 also have timing problem. To fix the clock drift, we determine P-wave arrival time at each station through a polarization analysis introduced by [54]. For a three-component airgun waveform in the passband of 2-30 Hz at station 53278 (Figure 6a), we first convert each component to an analytic signal.

Airgun Data Quality
We interpolate the airgun signals to 0.001 s/sample and bandpass all waveforms with a 2-6 Hz filter, the dominant frequency band of the bubble pulse [38]. With the single shot data, we can observe clear P-wave and S-wave arrivals up to 30 km with apparent velocities of 5.5 km/s and 2.8 km/s, respectively (Figure 4), consistent with previous observations [38].
When inspecting the waveforms throughout the entire year, we observe strong clock drifts in the reference station CKT and other receivers, approximately from day 230 in 2016 (Figures 5a,c and S3). Because the differential times are nearly the same among all stations ( Figure S4), it is mostly attributed to the timing log system of the airgun array. In additional to that, sensors in the 53273 and 53277 also have timing problem. To fix the clock drift, we determine P-wave arrival time at each station through a polarization analysis introduced by [54]. For a three-component airgun waveform in the passband of 2-30 Hz at station 53278 (Figure 6a), we first convert each component to an analytic signal.

Airgun Data Quality
We interpolate the airgun signals to 0.001 s/sample and bandpass all waveforms with a 2-6 Hz filter, the dominant frequency band of the bubble pulse [38]. With the single shot data, we can observe clear P-wave and S-wave arrivals up to 30 km with apparent velocities of 5.5 km/s and 2.8 km/s, respectively (Figure 4), consistent with previous observations [38]. When inspecting the waveforms throughout the entire year, we observe strong clock drifts in the reference station CKT and other receivers, approximately from day 230 in 2016 (Figures 5a,c, S3). Because the differential times are nearly the same among all stations ( Figure S4), it is mostly attributed to the timing log system of the airgun array. In additional to that, sensors in the 53273 and 53277 also have timing problem. To fix the clock drift, we determine P-wave arrival time at each station through a polarization analysis introduced by [54].

A Modified Moving-Window Cross-Spectrum method
We follow the method of Moving-Window Cross-Spectral (MWCS) [30,57] to compute the phase shift between stacked reference and current airgun signals. For a reference and a current airgun signal at station CKT with the center time at 0.5 s and window length of 0.6 s (Figure 7a), two segments are mean-adjusted and cosine-tapered before transforming into the spectral domain. The cross-spectrum, ( ), between the two segments is constructed as: where ( ) and * ( ) are Fourier-transformed representations of the windowed reference and current signals. * denotes complex conjugation. is the interested frequency range. We can represent the cross-spectrum, ( ), as: where ( ) represents frequency-dependent phase shift of a current signal relative to the reference one. represents a constant phase shift. To quantify the reliability of the estimations of phase shift, ( ) and , we compute the cross-coherence ( ) (Figure 7b) between their energy densities as: in which v(t), e(t), and n(t) are vertical, east-west and north-south components. H represents Hilbert transform. i is √ −1. We then define a 0.02 s moving window and a 0.01 s moving step and compute the covariance matrix between windowed three-components v a1 , e a1 , and n a1 as: where the asterisks represent complex conjugation. We compute three eigenvalues λ 1 , λ 2 , and λ 3 from the covariance matrix and calculate the rectilinearity R [55,56] through: It has a value of 1 when the signal is linearly polarized and a value of 0 when it is unpolarized. Due to interference of noise in the airgun signals, we cannot directly determine the P-wave arrival time from the maximum of rectilinearity in different moving windows. We first select a P-wave template and compute moving window correlation coefficient between the template and each raw vertical waveform. We then select points with correlation coefficient greater than 50% of the maximum value ( Figure 6b). Among these points, we determine an arrival time with the maximum differential rectilinearity (Figure 6c). From the arrival time to 0.5 s after it, we finally determine the P-wave arrival time from the maximum value of the rectilinearity ( Figure 6d).
As subtle subsurface changes cannot be reflected by direct P wave arrival changes, we then align P-wave arrivals through correcting travel time difference from the first shot. The airgun waveforms are significantly improved after clock drift correction (Figure 5b,d). We still observe that S-wave coda at 3-6 s ( Figure 5d) varies in amplitude and duration, which is attributed to the effect of changing airgun array.
We then stack all the waveform at each station for the entire year as the reference waveform ( Figure 5b,d), and derive differential time by a modified moving-window crossspectrum method that will be introduced in the following section. We cannot use the entire time-series to compute delay time between one single shot and the reference because of the changes in the airgun array. We thus investigate how delay time changes for the time-series between P-wave initial arrival and S-coda wave, which is 0.7-3 s in Figure 5d.

A Modified Moving-Window Cross-Spectrum method
We follow the method of Moving-Window Cross-Spectral (MWCS) [30,57] to compute the phase shift between stacked reference and current airgun signals. For a reference and a current airgun signal at station CKT with the center time at 0.5 s and window length of 0.6 s (Figure 7a), two segments are mean-adjusted and cosine-tapered before transforming into the spectral domain. The cross-spectrum, X( f ), between the two segments is constructed as: where F re f ( f ) and F * cur ( f ) are Fourier-transformed representations of the windowed reference and current signals. * denotes complex conjugation. f is the interested frequency range. We can represent the cross-spectrum, X( f ), as: where φ d ( f ) represents frequency-dependent phase shift of a current signal relative to the reference one. φ c represents a constant phase shift.
To quantify the reliability of the estimations of phase shift, φ d ( f ) and φ c , we compute the cross-coherence C( f ) (Figure 7b) between their energy densities as: where the overlines represent 5-point smoothing.
The values of φ c and φ d ( f ) can be derived from the unwrapped phase of the crossspectrum (Figure 7c). The frequency-dependent phase shift φ( f ) between the reference and the current airgun signals is linearly proportional to frequency φ d ( f ) = 2π f δt so that the cross-spectrum can also be represented as: where the time delay δt is estimated from the slope m of a weighted linear regression of phase φ d ( f ) in frequency of interest. We estimate a constant phase shift φ c from the intercept of the linear regression. The slope m (Figure 7c) is estimated as: where n is number of samples within the frequency range. The weight function w i for linear regression at the ith sample in the frequency of interest is: The associated error e m for the slope is: in which σ 2 φ is the squared misfit of the data to the modelled phase shift.
All the above steps are used to compute one data point at the time 0.5 s in Figure 7d,e. If we obtain delay-time δt and constant phase shift φ c estimations for all moving windows, we can further estimate how they vary with the travel time. For each measurement with window centered at j second, where b corresponds to the relative time change δt j /t j . The estimation for b is: where p j = 1/e 2 m and t = ∑ p j t j / ∑ p j . The error e b of the slope is Linear regression of phase shift φ c over the travel time constrains φ c /t (Figure 7d). Linear regression of delay time δt over different moving window constrains δt/t (Figure 7e). In the following section, we investigate how δt/t changes seasonally. Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 21

/ Derived Directly from Corrected Airgun Waveforms
A critical point in the MWCS analysis is the selection of the length of the moving window. The window length needs to be short to satisfy a constant time shift and to obtain a relative high temporal resolution. On the other hand, the window length needs to be long enough to ensure an accurate measurement, as well as avoid cycle skipping. Thus, the tradeoff between the temporal resolution and the accuracy of seismic velocity measurements can help to constrain the window length. Previous studies [17,58,59] suggested that the window length needs to be equal or larger than the longest period of interest, which is 0.5 s in this study. We computed / of P-wave at station 53278 ( Figure S5

δt/t Derived Directly from Corrected Airgun Waveforms
A critical point in the MWCS analysis is the selection of the length of the moving window. The window length needs to be short to satisfy a constant time shift and to obtain a relative high temporal resolution. On the other hand, the window length needs to be long enough to ensure an accurate measurement, as well as avoid cycle skipping. Thus, the tradeoff between the temporal resolution and the accuracy of seismic velocity measurements can help to constrain the window length. Previous studies [17,58,59] suggested that the window length needs to be equal or larger than the longest period of interest, which is 0.5 s in this study. We computed δt/t of P-wave at station 53278 ( Figure S5) using different lengths (i.e., 0.5, 0.6, 0.7, 0.8 s) of moving windows and found that the estimation of δt/t is more stable when the moving window is longer than 0.5 s. We thus use a moving window of 0.6 s to compute the differential time for all stations.
The other important parameter in linear regression of delay time over moving window is the window of P-wave arrivals ( Figure S6). Rather than selecting the begin and end time of P-wave from the raw airgun signal, we handpick them from the diagram of delay time verse travel time. We first visually inspect and determine the rough arrival time of P-waves. We then select the begin and end time of P-wave based on the linearity of delay time in different moving windows.
Because the arrivals of S-wave are not as obvious as those of P-waves (Figure 4), we only investigate how the arrival time of P-and P-coda wave changes. We observe that δt/t fluctuates between −5% to 5% for all stations (Figure 8) The error bars represent standard deviation of δt/t in one day. For δt/t measurements of P-wave arrivals on uncorrected waveforms, we observe strong uncertainties after the day 220 ( Figure S7). After clock drift correction, δt/t measurements show smoother seasonal change. Considering the uncertainty of correction imposed through this procedure, we conclude that δt/t measurements are more reliable for days before 220 than after. Although scattering varies at different stations, nearly all of them show an abrupt increase near days 150 and decrease around days 220 ( Figure 8). For closer stations from the source, such as 53277 and 53278, measurements on signals shot within the same day are constrained in a narrow δt/t range. For farther stations, such as 53276, 53264 and 53258, measurements have relatively large uncertainty. We can barely observe clear P-wave arrivals for stations farther than 20 km except the station 53266, in which we observe clear travel time changes.

/ Derived after Deconvolution of the CKT Waveform
Deconvolution has usually been used to remove the source effect of the airgun. For instance, [18] and [28] deconvolve a receiver station from a reference station to retrieve the empirical Green's function, based on which they measure the relative time delay to compute diurnal and semidiurnal seismic velocity change in the Binchuan basin. However, when the airgun excitation conditions change significantly, including the towing depth and firing pressure, the deconvolution cannot eliminate the source effect [43]. To investigate whether deconvolution can eliminate the source effect, we apply a frequencydomain water level deconvolution method [18,60] to construct empirical Green's functions and measure seismic velocity changes in P-wave.
in which and are Fourier spectrum of the receiver and source signals, respectively. We treat the recording at the reference station CKT as the source signal. The asterisk represents the complex conjugation. ranges from 0.01% to 1%, depending on the signal-tonoise ratio. We use a same value 0.01% as of [18] because of high-quality of our datasets. We finally invert frequency-domain Green's function to time-domain time-series

δt/t Derived after Deconvolution of the CKT Waveform
Deconvolution has usually been used to remove the source effect of the airgun. For instance, [18] and [28] deconvolve a receiver station from a reference station to retrieve the empirical Green's function, based on which they measure the relative time delay to compute diurnal and semidiurnal seismic velocity change in the Binchuan basin. However, when the airgun excitation conditions change significantly, including the towing depth and firing pressure, the deconvolution cannot eliminate the source effect [43]. To investigate whether deconvolution can eliminate the source effect, we apply a frequency-domain water level deconvolution method [18,60] to construct empirical Green's functions and measure seismic velocity changes in P-wave.
in which R and S are Fourier spectrum of the receiver and source signals, respectively. We treat the recording at the reference station CKT as the source signal. The asterisk represents the complex conjugation. c ranges from 0.01% to 1%, depending on the signal-to-noise ratio. We use a same value 0.01% as of [18] because of high-quality of our datasets. We finally invert frequency-domain Green's function G to time-domain time-series through inverse Fourier transform. If deconvolution can completely remove source effect, we should observe constant dominant frequency within one-year cycles, and seasonal change should be in the same magnitude with these values measured from ambient noise, which is usually less than 1%. After deconvolution, we observe that the dominant frequency of P-wave in these empirical Green's functions correlates inversely to the delay time change (Figure 9). The dominant frequency still shows seasonal variation, which causes obvious phase shift and affects the differential time calculation. Thus, deconvolution cannot eliminate the source effect on computation of differential time.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 21 through inverse Fourier transform. If deconvolution can completely remove source effect, we should observe constant dominant frequency within one-year cycles, and seasonal change should be in the same magnitude with these values measured from ambient noise, which is usually less than 1%. After deconvolution, we observe that the dominant frequency of P-wave in these empirical Green's functions correlates inversely to the delay time change (Figure 9). The dominant frequency still shows seasonal variation, which causes obvious phase shift and affects the differential time calculation. Thus, deconvolution cannot eliminate the source effect on computation of differential time. To investigate the effect of deconvolution, we also compare / of P-wave in original airgun waveforms and the empirical Green's functions ( Figure 10). We do not correct airgun data through procedures described above before retrieving of the empirical Green's function. If all receivers and the reference station share exactly same clock drift, the deconvolution will eliminate the clock drift effect. If not, we might observe large uncertainties of δt/t measurements on Green's functions. For station 53268, we observe that / from airgun signals correlates inversely to that from the empirical Green's functions.
For other stations, they are positively correlated to each other. For the stations 53268 and 53274, we observe large uncertainty of / at day 280 from the Green's functions, which may reflect different clock drifts in these stations from the reference station CKT. To investigate the effect of deconvolution, we also compare δt/t of P-wave in original airgun waveforms and the empirical Green's functions ( Figure 10). We do not correct airgun data through procedures described above before retrieving of the empirical Green's function. If all receivers and the reference station share exactly same clock drift, the deconvolution will eliminate the clock drift effect. If not, we might observe large uncertainties of δt/t measurements on Green's functions. For station 53268, we observe that δt/t from airgun signals correlates inversely to that from the empirical Green's functions. For other stations, they are positively correlated to each other. For the stations 53268 and 53274, we observe large uncertainty of δt/t at day 280 from the Green's functions, which may reflect different clock drifts in these stations from the reference station CKT. Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 21 Figure 10. Comparison of / before and after the deconvolution.

The Magnitude of / Change from Airgun
Almost 5% travel time change in P-wave from airgun measurement is huge compared to ~0.2% S-wave velocity decrease after the 2004 Parkfield earthquake [7] and ~0.4% increase due to thermoelastic strain increase in the San Jacinto fault area [17]. However, it is in the same magnitude with observations from [61] who observed an abrupt 4% coseismic velocity decrease for P wave after the Wenchuan mainshock. The reason why travel time change for P-wave is a lot larger than that for S-wave is still not well known. Before attempting to elucidate possible mechanisms for such discrepancy, we first need to inspect that our obtained / reflects subsurface variation or not.
Seismic travel time changes computed from ambient noise cross-correlation is a typical method to monitor the magnitude and trend of property change in the shallow crust. In the following section, we use ambient noise cross-correlation to compute seasonal change in seismic coda wave and compare its magnitude and trend with these from airgun data. If the magnitude and trend of seismic travel time changes from airgun are similar to these from ambient noise, we may suggest that the seismic travel time changes in P-wave from airgun are related to crustal property change. Otherwise, they are related to the source.

The Magnitude of δt/t Change from Airgun
Almost 5% travel time change in P-wave from airgun measurement is huge compared to~0.2% S-wave velocity decrease after the 2004 Parkfield earthquake [7] and~0.4% increase due to thermoelastic strain increase in the San Jacinto fault area [17]. However, it is in the same magnitude with observations from [61] who observed an abrupt 4% co-seismic velocity decrease for P wave after the Wenchuan mainshock. The reason why travel time change for P-wave is a lot larger than that for S-wave is still not well known. Before attempting to elucidate possible mechanisms for such discrepancy, we first need to inspect that our obtained δt/t reflects subsurface variation or not.
Seismic travel time changes computed from ambient noise cross-correlation is a typical method to monitor the magnitude and trend of property change in the shallow crust. In the following section, we use ambient noise cross-correlation to compute seasonal change in seismic coda wave and compare its magnitude and trend with these from airgun data. If the magnitude and trend of seismic travel time changes from airgun are similar to these from ambient noise, we may suggest that the seismic travel time changes in P-wave from airgun are related to crustal property change. Otherwise, they are related to the source.

δt/t from Ambient Noise
If seismic travel time change is due to shallow crust property change, we should observe similar trend of travel time change for P-wave from airgun and S-wave from ambient noise. To verify the travel time changes derived from the airgun source, we compute seismic travel time changes from ambient noise for coda waves in the passband of 2-6 Hz for 91 pairs of stations with the interstation distance less than 30 km (Figures 11  and S6). We first scan MiniSEED data into MSNoise [62], downsample to 20 Hz, remove earthquakes and airgun signals by the root-mean-square (RMS) temporal normalization, and reduce the effect of heterogeneous distributions by spectral whitening. We then compute daily CCs and construct a monthly CC through stacking and averaging CCs in the past one month before the selected date. We stack and average CCs in the entire year to construct a reference CC. We finally compute the delay time between monthly CCs and the reference CC with a moving window of 3 s and a step of 1 s. We do not follow a similar monthly stacking procedure in the previous airgun computation because of high similarity of waveforms of airgun. If similar procedure is followed in the computation of travel time changes for airgun data, the final seasonal variation will be smoother, but the magnitude and trend will not change significantly. We define the starting time of the coda window as one time of interstation distance and specify the width of the window as 30 s.
We compare the trave time changes in three pairs among station CKT, 53265 and 53268 (Figure 11a,b). Similar travel time changes among all station pairs reflect the robustness of δt/t derived from ambient noise. Indeed, 29 out of 91 station pairs show similar magnitude of δt/t, from −0.1% to 0.1% ( Figure S8b). Station pairs with seasonal changes in seismic travel time distribute in a broad region, up to 30 km from the reservoir ( Figure S8a).

/ . from Ambient Noise
If seismic travel time change is due to shallow crust property change, we should observe similar trend of travel time change for P-wave from airgun and S-wave from ambient noise. To verify the travel time changes derived from the airgun source, we compute seismic travel time changes from ambient noise for coda waves in the passband of 2-6 Hz for 91 pairs of stations with the interstation distance less than 30 km (Figure 11 and Figure  S6). We first scan MiniSEED data into MSNoise [62], downsample to 20 Hz, remove earthquakes and airgun signals by the root-mean-square (RMS) temporal normalization, and reduce the effect of heterogeneous distributions by spectral whitening. We then compute daily CCs and construct a monthly CC through stacking and averaging CCs in the past one month before the selected date. We stack and average CCs in the entire year to construct a reference CC. We finally compute the delay time between monthly CCs and the reference CC with a moving window of 3 s and a step of 1 s. We do not follow a similar monthly stacking procedure in the previous airgun computation because of high similarity of waveforms of airgun. If similar procedure is followed in the computation of travel time changes for airgun data, the final seasonal variation will be smoother, but the magnitude and trend will not change significantly. We define the starting time of the coda window as one time of interstation distance and specify the width of the window as 30 s.
We compare the trave time changes in three pairs among station CKT, 53265 and 53268 (Figure 11a,b). Similar travel time changes among all station pairs reflect the robustness of / derived from ambient noise. Indeed, 29 out of 91 station pairs show similar magnitude of / , from −0.1% to 0.1% ( Figure S8b). Station pairs with seasonal changes in seismic travel time distribute in a broad region, up to 30 km from the reservoir ( Figure  S8a). . dt/t(%) measurements (b) from ambient noise among CKT, 53265 and 53268 (a) and dt/t(%) measurements from airgun at the station 53268 (c). Figure 11. dt/t(%) measurements (b) from ambient noise among CKT, 53265 and 53268 (a) and dt/t(%) measurements from airgun at the station 53268 (c).

Possible Mechanism of δt/t from Ambient Noise
Travel time changes from ambient noise are usually related to environmental parameters changes, such as precipitation or water table fluctuation, temperature, and air pressure. To explore the correlation between the travel time changes and each environmental parameter (Figure 12), we requested such environmental data from one meteorological station in Lijiang, 150 km from the airgun source, from China Meteorological Data Service Center (http://data.cma.cn, accessed on 8 May 2021). We also compared the travel time changes to the water table fluctuation in the reservoir (Figure 12a).
We first observe that air pressure increases and correlates to the travel time change in the long-term trend (Figure 12b), which is opposite to [15], showing that seismic velocity increased with air pressure. Temperature increases in the summer and usually causes seismic velocity increase with a time lag of months [17]. However, we do not observe clear short-term and long-term correlation between the temperature and travel time changes (Figure 12d). Precipitation affects travel time changes through water mass loading or pore pressure changes. Water mass loading affects the travel time instantaneously, but the pore pressure causes travel time change with months of time lag. We compute a 10-day moving average value for the precipitation (Figure 12c). We observe that major precipitation (Figure 12c) occurs at day 200 in the summer. Seismic travel time changes increase to a local maximum at day 230 and lag behind the maximum precipitation for about one month, which is similar to the observation in [12].
To investigate if the pore pressure affects seismic velocity change, we use a poroelastic model to estimate seasonal variation of seismic travel time in 2016. The authors of [63] proposed a solution for the one-dimension fully coupled diffusion equation, which was used for understanding seismic travel time changes [12,59,64]. Assuming water diffuses in the vertical direction with constant hydraulic diffusivity, we can compute the pore pressure changes P(h, T) at h kilometer depth from the surface on day T as: in which n is the number of time increments δT from the day of the rainfall to the time T; δp i is the ground water load changes (ρ·g·δh i ) at the sampled instant T i = iδT. We compute δh i by using precipitation subtracting evaporation in 2016 ( Figure S9). δT is set to 1 day in second. er f c[·] is the complementary error function. c is the hydraulic diffusivity. The depth h is set to 500 m, which is the maximum depth sensitivity of surface waves in the frequency ranges of 2-6 Hz. A mathematical transfer function from [12,59,64] to compute a synthetic seismic travel time change: δt t syn (T) = δt t (T) + cov δt t (T), P(T) var(P(T)) ·(P(T) − P(T) ), (21) where δt t syn is synthetic seismic travel time change; δt t (T) is the average over observed seismic travel time change. We find the optimal diffusivity c value of 0.01 m 2 /s (Figure 12f) through minimizing the residual σ 2 (c): where n is number of days. Overall, the synthetic seismic travel time change correlates to the observed one (Figure 12e). The climatological data are from a Lijiang station, around 150 km from the study region. The discrepancy between synthetic results and the observed ones might be due to the uncertainty of the climatological data. We also investigated the effect of noise amplitude change on seismic travel time measurements. To compute the noise amplitude of ambient noise in the passband of 2-6 Hz, we remove airgun signals from the original dataset and calculate root-mean-square amplitude of the noise ( Figure S10). We observed that ambient noise amplitude increases in the summer when the water level in the reservoir decreases. The increasing noise amplitude in the summer might be related to active anthropogenic activities. We do not observe an obvious correlation between seismic travel time changes and variation of noise amplitude, so we conclude that changes in noise amplitude introduce negligible effects on velocity changes. We also investigated the effect of noise amplitude change on seismic travel time measurements. To compute the noise amplitude of ambient noise in the passband of 2-6 Hz, we remove airgun signals from the original dataset and calculate root-mean-square amplitude of the noise ( Figure S10). We observed that ambient noise amplitude increases in the summer when the water level in the reservoir decreases. The increasing noise amplitude in the summer might be related to active anthropogenic activities. We do not observe an obvious correlation between seismic travel time changes and variation of noise amplitude, so we conclude that changes in noise amplitude introduce negligible effects on velocity changes.
We finally propose that pore pressure change from a combined effect of precipitation and evaporation causes seismic travel time changes in coda waves in the ambient noise cross-correlations.

Difference of / from Ambient Noise and Airgun
Changes of 5% in travel time from airgun are a lot larger than 0.1% changes from ambient noise.
Changing of noise source property and shallow crustal structures are two major factors causing such change. A change of 0.1% seismic travel time may better reflect structural change in shallow crust than that from airguns. First, [12,17] proposed that changes in ambient noise source property cause a bias on travel time measurements, but the effect is negligible. Second, when measuring seismic travel time change for coda wave from ambient noise, we usually use a longer time window to include all possible scattered waves and multiple reflections, which better sample the crustal structure than direct P-wave arrivals. We need to understand the property change in the airgun source before investigating the correlation between / and subsurface structural change. We finally propose that pore pressure change from a combined effect of precipitation and evaporation causes seismic travel time changes in coda waves in the ambient noise cross-correlations.

Difference of δt/t from Ambient Noise and Airgun
Changes of 5% in travel time from airgun are a lot larger than 0.1% changes from ambient noise.
Changing of noise source property and shallow crustal structures are two major factors causing such change. A change of 0.1% seismic travel time may better reflect structural change in shallow crust than that from airguns. First, [12,17] proposed that changes in ambient noise source property cause a bias on travel time measurements, but the effect is negligible. Second, when measuring seismic travel time change for coda wave from ambient noise, we usually use a longer time window to include all possible scattered waves and multiple reflections, which better sample the crustal structure than direct Pwave arrivals. We need to understand the property change in the airgun source before investigating the correlation between δt/t and subsurface structural change.

Mechanism of δt/t Change from Airgun
To understand the effect of the change in source property on travel time measurements, we compare seasonal change in travel time to the dominant frequency ( Figure 13). We observe that δt/t is inversely correlated to the dominant frequency of P-wave and water level fluctuation in the reservoir. One possible explanation is that water table fluctuation in the reservoir causes signal frequency content change so that we can observe phase shift. The authors of [65] observed that signal amplitude changes are positively correlated to significant travel time delay, which is opposite to our observation. Since the change in source property cause significant effect on deriving seasonal change in seismic travel time, we should put more attention on stabilizing dominant frequency of airgun signals in further experiments.
To understand the effect of the change in source property on travel time measurements, we compare seasonal change in travel time to the dominant frequency ( Figure 13). We observe that / is inversely correlated to the dominant frequency of P-wave and water level fluctuation in the reservoir. One possible explanation is that water table fluctuation in the reservoir causes signal frequency content change so that we can observe phase shift. The authors of [65] observed that signal amplitude changes are positively correlated to significant travel time delay, which is opposite to our observation. Since the change in source property cause significant effect on deriving seasonal change in seismic travel time, we should put more attention on stabilizing dominant frequency of airgun signals in further experiments.

Conclusions
We computed seismic travel time changes / in the Binchuan basin and adjacent area, estimated seasonal change in the dominant frequency and signal amplitude, and determined how they correlate with the water level changes in the reservoir. The airgun signals exhibit 5.5 km/s and 2.8 km/s in apparent P-and S-wave velocities, respectively. Because S-wave is strongly interfered from multiple reflections and reverberations, we thus only investigated how seismic velocity of P-wave changes. We observed that the / correlates to the variation of dominant frequency and water table fluctuation in the reservoir. One possible explanation is that water table fluctuation affects the dominant frequency of airgun signals so that we can measure obvious phase shift. We also investigated seismic travel time changes in P-wave in the empirical Green's function and found that the / correlates inversely to the frequency change.
We computed 91 pairs of ambient noise cross-correlation to derive how seismic travel time changes in the passband of 2-6 Hz. We observed that the travel time changes increase to its local maximum in the end of August, but the variation of seismic velocity does not correlate with the changes in the water table in the reservoir. The seismic velocity changes lag behind the precipitation for about one month. We applied a poroelastic physical model to explain seismic travel time changes and found that a combined effect from precipitation and evaporation matches the observations. The optimal diffusivity to explain seismic travel time changes is 0.01 m 2 /s. Changes in ambient noise amplitude introduce negligible effects on seismic travel time changes.
Supplementary Materials: The supplementary materials are available online at www.mdpi.com/xxx/s1. Figure S1: Airgun excitation time from 2011 to 2020 (Up) and the number of average daily excitation (Down). From June to August, the water level in the reservoir is too low

Conclusions
We computed seismic travel time changes δt/t in the Binchuan basin and adjacent area, estimated seasonal change in the dominant frequency and signal amplitude, and determined how they correlate with the water level changes in the reservoir. The airgun signals exhibit 5.5 km/s and 2.8 km/s in apparent P-and S-wave velocities, respectively. Because S-wave is strongly interfered from multiple reflections and reverberations, we thus only investigated how seismic velocity of P-wave changes. We observed that the δt/t correlates to the variation of dominant frequency and water table fluctuation in the reservoir. One possible explanation is that water table fluctuation affects the dominant frequency of airgun signals so that we can measure obvious phase shift. We also investigated seismic travel time changes in P-wave in the empirical Green's function and found that the δt/t correlates inversely to the frequency change.
We computed 91 pairs of ambient noise cross-correlation to derive how seismic travel time changes in the passband of 2-6 Hz. We observed that the travel time changes increase to its local maximum in the end of August, but the variation of seismic velocity does not correlate with the changes in the water table in the reservoir. The seismic velocity changes lag behind the precipitation for about one month. We applied a poroelastic physical model to explain seismic travel time changes and found that a combined effect from precipitation and evaporation matches the observations. The optimal diffusivity to explain seismic travel time changes is 0.01 m 2 /s. Changes in ambient noise amplitude introduce negligible effects on seismic travel time changes.
Supplementary Materials: The supplementary materials are available online at https://www.mdpi. com/article/10.3390/rs13122421/s1. Figure S1: Airgun excitation time from 2011 to 2020 (Up) and the number of average daily excitation (Down). From June to August, the water level in the reservoir is too low so it is not suitable for the airgun to operate. Data is more complete in 2016 than other. The number of average daily excitation is 7 times/day in 2016, Figure S2: Continued, Figure S3: Different Airgun shots in 2016 at different stations. We observe strong clock drifts, Figure S4: Delay time of P-wave arrivals relative to the time in the first day of 2016, Figure S5: Seismic travel time change for station 53278 using different moving windows. We observe that the estimation of seismic travel time change is more stable when the length of the window is longer than 0.5 s, Figure S6: Determination of window for P-wave arrivals. We inspect the waveform to get a rough estimate of the begin and end time of P-wave arrivals. We then handpick the exact arrival time of P-wave based on the linearity of data points in the diagram of delay time verse travel time in different days, Figure S7: δt/t measurements on P-wave arrivals at four stations with or without clock drift correction, Figure S8: (a) Station pairs with similar seasonal changes of seismic travel time distribute in a broad region, at most 30 km from the reservoir. (b) 29 out of 91 station pairs show similar seismic travel time change δt/t. Figure S9: Precipitation substracts evaporation. Figure S10: Ambient noise amplitude increases in the summer when the water level in the reservoir decreases. The increasing noise amplitude in the summer might be related to active anthropogenic activities.