Ground water-induced changes in velocities of P and S waves (Vp and Vs) measured using an accurately controlled seismic source

We analyze temporal variations in the travel times of both P and S waves (Vp and Vs) for 14 months at Toyohashi (central Japan) with a continuously operating vibration source that could efficiently produce both P and S waves. Seismic waves produced by the source, which is called the accurately controlled routinely operated signal system (ACROSS), are recorded by three nearby seismic stations, and the travel time variation at each station is estimated using the transfer function calculated from the recorded data. Long-term variations in Vp and Vs are observed and can be interpreted by the change in fluid saturation and crack density of the subsurface rocks. The variation in fluid saturation and crack density are consistent with that in the groundwater level, which is measured at the station nearest to the ACROSS. Short-term responses to rainfalls are observed at the station nearest to the ACROSS system; the interpretation of the changes in crack density and saturation is inconsistent with the ground water observation, partly owing to the initial response to rainfall. This can be interpreted as an air–water mixture within pores or cracks on a fine scale.


Introduction
Changes in atmospheric conditions, such as rainfall and air temperature, result in changes in seismic velocity, and many studies have been conducted on the effects of other factors on seismic velocity. Evaluation and quantitative estimation of such effects is also necessary for detecting the velocity change in the deeper part of the crust, which is also a subject of geoscience. Rainfall is the most effective source of change in seismic velocity because it results in the direct infiltration of water to fill the pores of the subsurface soil or rock (Martini et al. 2009). Changes in air temperature have recently been found to influence the subsurface seismic velocity (Richter et al. 2014;Wang et al. 2020). This may be difficult to understand considering the short penetration depth of atmospheric temperature variation, although the effect of thermal strain (Wang et al. 2020) should be studied by further observations. Seasonal velocity changes are generally observed in relation to groundwater level variation (Clements and Denolle 2018;Sens-Schönfelder and Wegler 2006). In addition to atmospheric conditions, strong earthquakes also affect seismic velocity (Nakata and Snieder 2011;Niu et al. 2008). This was initially observed by Ikuta et al. (2002) and is now a well-known phenomenon.
The seismic velocity of the subsurface medium changes with variations in the rock conditions, such as increase or decrease in cracks or fluids within the rocks. Nur and Simmons (1969) found the velocity change of seismic waves to be a function of the pressure of dry and saturated rocks through laboratory experiments. King et al. (2000) also investigated the velocity change with increasing saturation via laboratory experiments. The results of both studies showed that the effect of fluids on seismic velocity differs for P and S waves. O'Connell and Budiansky (1974) theoretically proved that seismic velocity depends on the spatial density and fluid saturation ratio of the cracks in subsurface rocks. Their results showed that the P-wave velocity (Vp) and S-wave velocity (Vs) differ in terms of their dependence on the cracks and fluids in subsurface rocks; therefore, the state of cracks and fluid within the subsurface medium can be evaluated using Vp and Vs variations.
Both natural and artificial sources were used to observe the temporal variation in the seismic velocity. Natural sources, such as ambient noise, are now commonly used in combination with seismic interferometry (Brenguier et al. 2008;Liu et al. 2014;Hobiger et al. 2016). Seismic interferometry has also been applied to coda waves (Grêt et al. 2006) or shear waves (Miao et al. 2018) of natural earthquakes to detect seismic wave velocity changes. Artificial sources, such as vibroseis (Karageorgi et al. 1992), piezoelectric actuators (Yamamura et al. 2003;Silver et al. 2007), and seismic airguns Wang et al. 2020), have also been used to study the temporal variation in seismic velocity with high resolution. However, it is not easy to extract P waves from the ambient noise (Nakata et al. 2019), making it difficult to estimate the fluid effect that appears mainly on Vp.
In our study, an artificial seismic source, called Accurately Controlled Routinely Operated Signal System (ACROSS), is used to monitor temporal variations in seismic velocities of P and S waves. It generates vibrations by precisely controlling the rotation of eccentric weights. Because ACROSS source has high temporal stability and can therefore obtain a high signal-to-noise ratio (SNR) by stacking, it has the advantage of continuous and highresolution monitoring of velocity variations (Yamaoka et al. 2001). Ikuta et al. (2002) and Ikuta and Yamaoka (2004) monitored travel time variation for 15 months near a fault that caused the Kobe earthquake in 1995 and detected changes due to earthquake-induced strong ground motion. Saiga et al. (2006) analyzed 1-monthlong records to show that the signal was sensitive to temperature variations and rainfall. Additionally, Maeda et al. (2015) found a difference in transfer functions between the active and inactive periods of Sakurajima volcano, Japan. Furthermore, Tsuji et al. (2018) detected coseismic and secular changes in the Tokai region by examining the travel time variations for 10 years.
We measured travel time variations of both P and S waves by taking advantage of a new type of ACROSS source and interpreted the variations by obtaining responses to the fluid and cracks within the subsurface medium. This new type of ACROSS source (see Fig. 1) has a horizontal rotation axis to produce both vertical and horizontal forces, in contrast to conventional ACROSS sources, which have a vertical rotation axis. Therefore, the new type of ACROSS source effectively excites both P and S waves, and can extract P waves from the signal better than seismic interferometry does from the ambient noise. The source was deployed at an observatory of Nagoya University in March 2014 and began continuous operation at the end of 2017 after performance tests. We analyzed the signals from the continuous operation of ACROSS to monitor travel time variations of both P and S waves and discussed the influence of rainfall and groundwater level on the observed variations.

Source and data
We monitored travel time variations of P and S waves from December 8, 2017 to February 7, 2019 using the signals from ACROSS observed at three seismic stations located around the source. We also compared the change in travel time with the rainfall and groundwater level data to interpret the origin of the temporal change in the seismic velocity. In this section, we briefly overview the ACROSS source and experimental sites used in our study.

Locations of the ACROSS source and observation stations
The locations of the ACROSS source and observation stations are shown in Fig. 2. The ACROSS source was installed at the Mikawa Observatory of Nagoya University in Toyohashi City (central part of Japan). The seismic signals generated by the ACROSS source were monitored by three observation stations: NU.MIK (Mikawa), N.THNH (Toyohashikita), and N.MKBH (Mikkabi). NU.MIK is located in a vault of the Mikawa Observatory and equipped with a Streckeisen STS-2 seismometer, which is a three-component broadband velocity-type sensor; the distance of the station from source is 170 m. N.THNH and N.MKBH belong to Hi-net stations Fig. 1 ACROSS source used in this study. Eccentric weights rotate around a horizontal axis to produce a sinusoidal signal. Both horizontal and vertical forces generated by the rotation provide efficient excitations of P and S waves (National Research Institute for Earth Science and Disaster Resilience National Research Institute for Earth Science and Disaster Resilience (2019) and are equipped with short-period velocity sensors in boreholes at depths of 203 and 103 m, respectively. The distances from the source to N.THNH and N.MKBH stations are 3.4 and 6.0 km, respectively. We used the rainfall data from the Automated Meteorological Data Acquisition System (AMeDAS) of the Toyohashi and Mikkabi stations operated by the Japan Meteorological Agency (JMA). Ground water levels were observed at three sites, TYH1, TYH2, and TYE, which were located adjacent to the source (see Fig. 2) by the National Institute of Advanced Industrial Science and Technology (AIST). The TYH2 station sampled the water pressure at a depth of 130-150 m, and the THY1 and TYE stations sampled the water pressure at a depth of between 180 and 200 m. The details of this observation are summarized in Matsumoto et al. (2007).

Signal generation using the ACROSS source
ACROSS sources generate sinusoidal waves of various frequencies by rotating the eccentric weights. Frequency modulations (FM) are applied to the rotation to generate signals in a certain width of the frequency band. Because the frequency modulations are repeated in a precise constant time interval (FM period), the generated signal is expressed by the sum of sinusoids, whose frequencies are integer multiples of the reciprocal of the FM period.
We used an FM period of 50 s, with an upsweep of 37.5 s and a downsweep of 12.5 s. The modulation frequency was centered at 10.005 Hz and ranged from 5 to 15 Hz. Therefore, this modulation produced sinusoidal signals at 10.005 Hz ± 0.02, 0.04, 0.06 … Hz between approximately 5 and 15 Hz. The source switched its rotation direction every 2 h to synthesize two perpendicular linear forces after the completion of data acquisition. Because the eccentric weight of the ACROSS source rotates around the horizontal axis, we can synthesize both vertical and horizontal linear forces that are oriented perpendicular to the rotation axis. Thus, P and S waves can be efficiently excited using the vertical and horizontal linear forces that can be separated by this operation. Technical details of the analysis method can be found in the textbook by Kasahara and Hasada (2016).

Availability of data
Before analysis, we removed data that were unusable owing to troubles with the source or seismic stations. We first picked out the days when the ACROSS vibrator was not working correctly by looking at the operation logs of the source acquired every second. Furthermore, we compared the logs of the rotation frequencies and angular positions of the eccentric weights with those in the period in which we considered that ACROSS stably operated. We chose the day of reference as June 17, 2018, and removed the days that included mean deviations of more than 0.001% for frequency and 0.12% for angular position. Consequently, we could remove the days when the source was under maintenance or stopped because of power outages. In addition to the problems in the source The map data were provided by the Geospatial Information Authority of Japan operation, we selected the days that indicated disorders in the seismographs. Sensor check signals included in the seismograph of the Hi-net stations (N.THNH and N.MKBH) were used to identify sensor disorders using the method described by Kunitomo (2014). We calculated the cross-correlation of the sensor-check signal by referring to the data obtained on December 1, 2017. The reference signal was chosen from the period in which stable cross-correlation was observed. Data with a crosscorrelation coefficient of less than 0.8 were excluded. For the NU.MIK station, which did not have a sensor-check signal, we assumed that no trouble originated in the sensor because a feedback type sensor (STS-2 seismometer) was deployed in the vault. The usable days of the UD (up-down), NS (north-south), and EW (east-west) components of each observation station are listed in Table 1. Finally, 87% of the data obtained during the monitoring period could be used in our study, excluding the NS component of the N.MKBH station, which mostly suffered from sensor trouble.

Data processing
We calculated the travel time variations from the transfer functions between the observed waveforms at the seismic stations and the source function of the ACROSS vibrator. The transfer functions were calculated using the standard procedure described by Kasahara and Hasada (2016). The travel time variations were calculated using the method described by Ikuta et al. (2002). We now briefly explain the data processing to obtain travel time variations for both P and S waves, along with data stacking and transfer function calculations.

Stacking of observation data
We made the stacking of the observed waves to enhance the SNR and improve the detectability of the temporal variation in travel time. In the stacking, we used the weighted stack method proposed by Nagao et al. (2010), as shown in the following equations: for the data of the kth day [ X f , k ] at frequency f . K represents the number of usable dates for the analysis. The weight is given by the noise variance [ ε f , k 2 ] , as shown in Eq.
(2). Stacking is performed on the data of every 400-s-long section. Because each section includes eight cycles of frequency modulation, the signal of the ACROSS source appears every eight components in the frequency series; therefore, other components include only noise. The noise variance is assumed to be equivalent to the squared amplitude of the frequency component adjacent to the component of the ACROSS signal. ε s f in Eq.
(3), which is the noise level after stacking, is used in the following analysis to obtain the estimation errors in travel time variation.

Calculation of transfer functions
Transfer functions, which represent the subsurface structure, are obtained by deconvolution of the daily stacked data using the force generated by the ACROSS source. We used the force theoretically calculated employing the operation parameters of the ACROSS source. Subsequently, we calculated the daily transfer functions for the clockwise and counter-clockwise rotations and synthesized the transfer functions for both vertical and horizontal source excitations. (1) the calculation of travel time variation, as described in the following section. The transfer functions are labeled with the components of the receiver and source excitation. For example, the vertical motion excited by vertical and horizontal forces is labeled as Uv and Uh, respectively. Nv and Ev, and Nh and Eh represent the northward and eastward motions excited by the vertical and horizontal forces, respectively. The direction of the horizontal excitation is indicated by the black arrow in Fig. 2. In our study, two horizontal components at the seismic stations were linearly combined to produce a horizontal component, whose orientation was parallel to that of the horizontal excitation at the source. The component is labeled as Hh in the following section.

Travel time variation of P and S waves
We calculated the daily variation in travel time using the method described by Ikuta et al. (2002). We calculated a cross-spectrum for both P and S waves between the transfer function recorded for each day and the reference transfer function, which is shown in Fig. 3. The travel time difference ΔT was calculated from the phase difference as shown in Eq. (4), in which a weighted average over all spectral components is adopted. The uncertainty δt N is estimated from the noise level in the frequency domain using Eq. (5), which can be interpreted as the measurement error of this experiment. The uncertainty arising from Eq. (4) is Eq. (6), which is the standard error in estimating the T as the weighted average: where A n represents the square root of the amplitude of the cross-spectral density, θ n the phase difference from the reference, and ω n the angular frequency in the nth frequency component. X n and ∼ σ n represent the Fourier transforms of the transfer function and noise at the nth frequency, respectively. The summations are applied over the frequency components of the ACROSS signal used. Equation (5) is corrected from the equation proposed by Ikuta et al. (2002).
We selected P and S waves by applying Hanning windows to the transfer functions, as shown in Fig. 3 and Table 2. The positions and lengths of the windows were determined such that both the initial arrival and maximum amplitude of the P or S waves were included. We used the Uv component for the P waves and the Hh component for the S waves. The horizontal direction of the Hh component was determined to have the largest amplitude in the horizontal particle motion, as shown in Fig. 3d-f. For the N.MKBH station, however, we used the Eh component without synthesizing of Hh component because the sensor trouble in the NS component seriously reduced the number of available days. We decided to use the Eh component rather than Hh because the direction of the maximum amplitude was almost in the east-west direction.

Results and discussion
We now present the results of our analysis and discuss the influence of rainfall on seismic velocity variations. We examined the dependence of V P and V S on the crack density and fluid saturation using the theoretical relation provided by O'Connell and Budiansky 1974). Considering SNR, which decays with distance, we stacked the data for 60 days for the N.THNH and N.MKBH stations rather than using 1-day stacked data, which were used for the NU.MIK station.

Temporal variation in travel time
The temporal variation and uncertainties in the travel time obtained in our calculations are shown in Fig. 4ac. In temporal variations, both short-term changes, which continue for several days after rainfall events, and long-term changes, which appear throughout the observation period, can be recognized. First, we describe the variations for NU.MIK, which is located closest to the source and ground water measuring facilities. In Fig. 4, the high precipitation period from April to September in 2018 is shown by orange lines, and the low precipitation periods from December 2017 to March 2018 and from October 2018 to February 2019 are shown by blue lines. At NU.MIK, many short-term changes in the travel time are seen as step-like delays, with a maximum step of approximately 0.3 ms, followed by a gradual recovery. Comparing these changes with the precipitation data, we found that they were preceded by rainfall, and significant changes were seen in the high precipitation season in the given year. The magnitude of this change was larger in the P waves than in the S waves. The long-term variation was also higher for the P waves than for the S waves. A large advance of P wave was observed in late June 2018, whereas a small long-term change was observed in the S wave during this period. However, in May 2018, a large delay in the S wave and a small change in the P wave were observed. The magnitudes of these changes were approximately 0.9 and 0.6 ms for the P and S waves, respectively. This is notable considering the overall change during the entire period of approximately 1.2 and 0.7 ms for the P and S waves, respectively.
Figures 4b, c shows temporal variations in the travel times at the N.THNH and N.MKBH stations, for which the source-receiver distance is considerably larger than that in the case of NU.MIK. Because of the low SNR, a moving average of 60 days was applied for these stations. A 60-day-long window was shifted at a 1-day step, and averaged data were plotted at the middle of the period when the window was applied. The data were not plotted if the data deficit within the window exceeded 20 days. In these plots, the short-term variations were suppressed but larger variations could be seen for the P wave (more than that seen for the S wave), even upon considering the difference in uncertainties between the P and S waves.

Temporal relationship between velocities of P and S waves (Vp and Vs, respectively)
Figure 5a-d shows the relationship between the changes in Vp and Vs. Figure 5a-c shows the plots for all the periods in Fig. 4a-c. Figure 5d shows the representative variation during the period of rainfall indicated by the orange bar in Fig. 4a, d, e. The velocity change rates δV /V were calculated from the travel time variation as follows: where δT represents the travel time variation, and T was chosen as the travel time at the middle of each time window shown in Fig. 3, which almost corresponds to the part of the maximum amplitude. The orange and blue lines in Fig. 5 correspond to those in Fig. 4. The velocity variation at NU.MIK (see Fig. 5a) mainly comprises three trends: a decrease in Vp with a small δV /V = −δT /T , increase in Vs during the low precipitation period, a decrease in Vs without a major change in Vp at the beginning of the high precipitation period, and an increase in both Vp and Vs during the high precipitation period. The amplitude of variation was larger in Vp than in Vs for both long-and short-term variations. The N.THNH station (see Fig. 5b) shows a trend in which only Vp variation predominates. A small positive correlation between Vp and Vs is observed in this plot. This is also the case for N.MKBH, except for the small negative correlation between Vp and Vs. Figure 5d shows the typical response of Vp and Vs to rainfall at NU.MIK, corresponding to Period 1 (June 19-July 4, 2018) shown in Fig. 4. The travel times of both the P and S waves changed in response to rainfall, as shown in Fig. 4, although their responses were different. In the P wave, a sharp bend appears in the rainfall response, whereas the S wave shows a blunt bend. The rainfall-relevant change can be divided into three parts, as seen in Fig. 5d: Vp decreased immediately after the rainfall with a small decrease in Vs; in the first half of the recovery response, Vp increases with a small change in Vs; in the second half, both Vp and Vs increase. This type of temporal pattern is common for most responses to rainfall observed in NU.MIK.

Influence of crack density and fluid saturation
O'Connell and Budiansky (1974) theoretically calculated the Vp and Vs of a cracked medium. They formulated the dependence of seismic velocity on crack density and fluid saturation and showed the difference between Vp and Vs. Figure 6 shows the temporal variation of Vp and Vs obtained in our study, plotted on the diagram using the equations derived from O'Connell and Budiansky (1974). The base diagram shown in Fig. 6h indicates the change in Vp by green lines and that in Vs by yellow lines, as a function of crack density and fluid saturation (see Appendix 1). This figure provides an opportunity to understand the effect of saturation and crack density on Vp and Vs. Vp and Vs in these graphs represent the velocity ratios with respect to an intact (crackless) medium. The dependence on fluid saturation and crack density is different for Vp and Vs.
We interpreted the changes in Vp and Vs, which are shown in Fig. 5, by converting them into the fluid saturation and crack density, as shown in Fig. 6. Although it is difficult to determine the seismic velocity for the intact condition of the country rock through which the seismic wave propagates, we tried to interpret it by assuming the initial ratios of Vp and Vs for those of the intact medium to be 0.6. Once the initial ratio is assumed, we can convert the velocity variation on the Vp-Vs plane onto the crack density-fluid saturation plane. In this process, we used the following equation: where ε and ξ denote the crack density and fluid saturation, respectively. J is generally called Jacobian. In the conversion, J can be calculated as the slope of Vp and Vs in the direction of ε and ξ for the given initial ratio of Vp and Vs. From this equation, we can obtain each (dε, dξ) on each day from (dV P , dV S ) using J −1 .
In the following sections, we show the changes in fluid saturation and crack density calculated from the changes in Vp and Vs, and we also examine whether we can infer the changes in the fluid saturation and crack density by monitoring Vp and Vs.

Long-term variation
For all the stations, the observed variations in Vp and Vs almost follow a linear trend, in which the fluid saturation and crack density have a positive correlation. This trend results from the observation that the variation in Vp is considerably larger than that in Vs.
In NU.MIK, the saturation and crack density decreased during the dry season. This may be interpreted as a crack closure owing to decrease in fluids in the cracks. The crack density slightly increased at the beginning of the wet season, and both the fluid saturation and crack density increased in the following period. This may reflect fluid infiltration by rainfall and the resultant crack openings. However, N.THNH and N.MKBH showed small difference between the dry and wet seasons. Because these two stations were located farther than NU.MIK from the ACROSS source, the P and S waves passed through  Fig. 4. Furthermore, d shows short-term velocity change associated with the rainfall indicated in Fig. 4 for NU.MIK. The circular and triangular marks indicate the beginning and end of the period, respectively the deeper part of the crust. This may indicate that the contribution of the variation in transfer function in the deeper part for N.THNH and N.MKBH is larger than that for NU.MIK because the seasonal difference in velocity might be smaller in the deeper part than in the shallower part.

Effects of short-term rainfall
As seen in the precipitation data observed by AMeDAS, the selected location experiences heavy rainfall every year. Therefore, it is interesting to discuss the temporal responses of Vp and Vs to a single rainfall event. Figure 6d shows a typical example of a changing cycle associated with a rainfall event. The cycle is plotted on the diagram with the saturation and crack density for NU.MIK, corresponding to Fig. 5d. The changing cycle comprises three stages, as explained in Sect. 4.2. In the first stage, immediately after rainfall, both saturation and crack density decreased, corresponding to the period in which Vp decreased (with a slight change in Vs). In the second stage, both saturation and crack density increased corresponding to the period in which Vp increased (with a small decrease in Vs). In the third stage, both saturation and crack density continued to increase, with a larger increase in saturation. These changes in saturation and crack density are the result of sudden delay and recovery of both P and S waves. Similar changes in travel times of P and S waves can be seen in Fig. 4 frequently throughout the observation period. Other examples for the rainfall events, which have high precipitation and certain period before next rainfall, are also analyzed for confirmation. We analyzed the data in Periods 2-4 of Fig. 4 and the plots are shown in Fig. 6. The variations for all the periods show a decrease in both saturation and crack density immediately after the rainfall, similar to that in Period 1. These interpretations seem strange because the fluid saturation may increase as a response to rainfall and decrease subsequently. This incompatibility is discussed in the next section.

Relationship with change in ground water level
We examined the relationship between variations in Vp and Vs, and ground water level, measured by observational wells located near the ACROSS site, as shown in Fig. 2. We attempted to interpret the ground water level variation using the fluid saturation and crack density, which are presumed by velocity changes, as stated in the previous sections. As shown in Fig. 4d, TYH2 records short-term changes, such as rainfall response, whereas TYH1 and TYE are insensitive to rainfall and may reflect the pressure of deeper aquifers. Therefore, we used TYH2 to compare the short-term changes. For comparison of long-term variations, only TYH1 was used because the records of TYH1 and TYE were almost the same, and the  (1974), which is shown in h as a base diagram. The green and yellow solid lines show the ratios of Vp and Vs to those of the intact medium, respectively. a-g Shows the variations obtained in our experiment. The plotted region in these figures corresponds to the small rectangle of h, in which we assume the initial velocity ratios for both P and S waves to be 0.6. a-c Show the long-term variations at NU.MIK, N.THNH, and N.MKBH, respectively. d-g Show the short-term changes associated with the rainfall events at NU.MIK. The thin and dotted lines in a-g represent the contours for the value in every 0.01 and 0.005, respectively long-term trend of TYH2 was similar to that of THY1. The comparison was made only for NU.MIK because the observation well was located close to the ACROSS source and NU.MIK. The long-term variations in Vp and Vs can be compared with the ground water level of TYH1 (see Fig. 4d). Generally, the wet season corresponds to a period of high groundwater levels. When the wet season begins, the groundwater level increases and remains high until the dry season begins. In the period of high groundwater level, Vp gradually increases (travel time decreases for P) with small variations in Vs. The variation in seismic velocities during this period result in an increase in both the crack density and fluid saturation of the subsurface rocks, and this may indicate a continual diffusion of water in the ground. In the dry season, the ground water level decreases and Vp continually decreases with a small change in Vs. This indicates that the crack density and fluid saturation decrease with the water level. A clear decrease in Vs (travel time increase of the S wave) is observed in the period around the beginning of the wet season, i.e., from March to May 2018. Considering the decrease in Vp during this period, this change can be interpreted as an increase in the crack density and a small change in the saturation.
The short-term velocity change (Period 1) can be compared with the data for the shallower aquifer at TYH2. Although TYH2 sampled the water pressure at a depth of 130-150 m, the aquifer was efficiently connected to the shallow part where the rainwater could easily penetrate. In contrast to the long-term variation, the ground water level increased when the saturation and crack density decreased, as indicated by the decrease in Vp.
However, the interpretation in terms of crack density and water saturation, especially during the first stage of rainfall, seems inconsistent with the groundwater observations. One possible explanation is the limited application of the study conducted by O' Connell and Budiansky (1974) with respect to the velocity response to rainfall. Mavko and Nolen-Hoeksema (1994) discussed the effect of a gas-liquid mixture of pore fluid on the bulk Vp using Gassmann's relation. On a fine scale and within each crack, the bulk modulus of the gas-liquid mixture was almost equivalent to that of the gas, even for a small gas content, whereas the density increased with the liquid content. Therefore, the bulk Vp decreased because of the increase in density without any change in the bulk modulus. The relation in O'Connell and Budiansky (1974) does not consider this, but treats the fluid saturation as a ratio of saturated cracks to the total number of cracks.
At the initial stage of rainfall, the water infiltrates the ground, thereby increasing the water content of each crack in the shallow subsurface. This may increase the density of the medium without changing the bulk modulus, decreasing both Vp and Vs. At the later stages, Vp increased owing to the decreasing water, as indicated by the ground water level.

Conclusion
We analyzed travel time variations of both P and S waves using a new type of ACROSS vibration source that efficiently excited both P and S waves. From the analysis, the long-and short-term variations resulting from water infiltration in the ground could be recognized. The variation in Vp was larger than that in Vs at all the observation stations. The observed variations in Vp and Vs were then converted into variations in the crack density and fluid saturation of subsurface rocks using the relation described by O'Connell and Budiansky (1974). By comparing the velocity variations with the ground water level observed near the source, we found that the long-term variations in the ground water were consistent with the interpretation of variations in the fluid saturation and crack density; however, the short-term variations associated with rainfall events were inconsistent. This may be interpreted as an air-liquid mixture of fluids in individual cracks. The velocity change at shallow depths, estimated by converting the observed Vp and Vs into fluid saturation and crack density in our method, will be used to eliminate the effect of the shallow part, to extract the change in the deeper part.

Appendix 1
Drawing diagram on velocity change from O' Connell and Budiansky (1974).