Temporal change in seismic wave attenuation using highly stable vibration sources

We developed a method to detect attenuation changes during seismic wave propagation excited by precisely controlled artificial seismic sources, namely Accurately Controlled Routinely Operated Signal System (ACROSS), and applied it to monitor the temporal changes for in situ data collected by previous studies. Our method, together with the use of the ACROSS sources, is less susceptible to noise level changes, from which conventional methods such as envelope calculation suffer. The method utilizes the noise level that is independently estimated in the frequency domain and eliminates the influence of the noise from the observed signal. For performance testing, we applied this method to a dataset that was obtained in an experiment at Awaji Island, Japan, from 2000 to 2001. We detected a change in amplitude caused by rainfall, variation in atmospheric temperature, and coseismic ground motions. Among them, coseismic changes are of particular interest because there are limited studies on coseismic attenuation change, in contrast to many studies on coseismic velocity decrease. At the 2000 Western Tottori earthquake (MW = 6.6, epicenter distance of 165 km), a sudden decrease in amplitude of up to 5% was observed. The coseismic amplitude reduction and its anisotropic characteristics, which showed a larger reduction in the direction of the major axis of velocity decrease, were consistent with the opening of fluid-filled cracks, as proposed by previous studies. The ΔQ-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {Q}^{-1}$$\end{document} corresponding to the amplitude change gives similar values to those reported in previous studies using natural earthquakes.


Introduction
Temporal changes in the propagation property of seismic waves have been studied using both seismic velocity and attenuation. In previous studies, coseismic changes, among others, have drawn the attention of researchers because they display clear synchronization to earthquakes and their relation to relevant phenomena can be discussed quantitatively. A decrease in seismic velocity is often reported as a typical phenomenon in earthquakes. Various methods have been used to detect seismic velocity changes at earthquakes in many locations. Seismic interferometry (Brenguier et al. 2008;Hobiger et al. 2016;Ikeda and Tsuji 2018), which does not require artificial sources, is generally used to measure the temporal variation in seismic velocity in recent years. In contrast, experiments with artificial seismic sources (Ikuta et al. 2002;Yang et al. 2014;Tsuji et al. 2018) were conducted to identify coseismic velocity changes with better resolution. In most cases the seismic velocity decreases during the earthquakes and recovers gradually over time. These phenomena are accounted for by the fracturing of subsurface rock (Sawazaki et al. 2009;Nakata and Snieder 2011) and/or crack opening by stress (Grêt et al. 2006;Silver et al. 2007) or pore pressure changes (Ikuta and Yamaoka 2004;Sawazaki and Snieder 2013) caused by the strong ground motions of earthquakes.
Fracturing, crack opening, and saturation increase of pore fluid cause attenuation, as shown in laboratory experiments (Lockner et al. 1977;Toksöz et al. 1979;Tao and King 1990;Best et al. 2007;Chichinina et al. 2009;Amalokwu et al. 2014;Zaima and Katayama 2018). These studies measured the changes in attenuation with fluid saturation and confining pressure, and reported an increase in attenuation with crack opening and saturation. Attenuation with cracks is interpreted as an increase in friction loss on the crack surface and/or energy loss due to fracture to generate cracks. The attenuation with saturation is interpreted as an increase in viscous losses in the cracks and/or with fluid flowing through the microchannels between cracks.
Coseismic changes in attenuation were also detected in the field. The spectral ratio of natural earthquakes (Chun et al. 2004;Titzschkau et al. 2010;Kelly et al. 2013) and coda-Q (Sato 1986;Fehler et al. 1988;Huang and Kisslinger 1992;Domínguez et al. 2003;Sugaya et al. 2009;Padhy et al. 2013) have been used. They are used to estimate the temporal variation in attenuation in the region surrounding earthquake faults or in volcanic regions. The analysis based on the spectral ratio assumed a model in which Q is independent of frequency, for example, because the amplitudes of the sources are difficult to estimate precisely. However, when we use artificial sources, the assumption may not be necessary. Yamamura et al. (2003) used a piezo-electric actuator, which produces seismic signals of a constant amplitude, as a seismic source to reveal the periodic variation in attenuation associated with earth tides. Despite the use of seismic sources with high resolution, no remarkable earthquakes occurred during their experiment, and they did not detect a coseismic change of attenuation. Nibe and Matsushima (2021) detected temporal change in attenuation in reservoir of a geothermal field associated with shut-in of production wells. Zhu et al. (2017) detected change in attenuation in a test field of geological CO 2 storage associated with injection of CO 2 . These studies used centroid frequency shift method (Quan and Harris 1997) to detect changes in attenuation.
We developed effective methods for high stability artificial sources to measure the temporal change in amplitude of seismic signals. We used the Accurately Controlled Routinely Operated Signal System (ACROSS) as a highly stable artificial source, which makes it possible to obtain accurate transfer functions. The stability of the ACROSS source is a principal advantage of measuring the temporal variation in attenuation, which is difficult to acquire when using natural seismic sources. The ACROSS was developed for precise monitoring of propagation properties through continuous seismic waves excited by a precisely controlled rotating eccentric weight (Kumazawa and Takei 1994), which enables us to continuously obtain precise transfer functions between sources and receivers (Kumazawa et al. 2007). With these characteristics, ACROSS is able to monitor small variations in the medium in various fields, such as the Sakurajima Volcano, Japan (Yamaoka et al. 2014;Maeda et al. 2015), or the Nankai Trough subduction zone .
Another advantage of estimating amplitude variation using ACROSS is the noise level estimation during the signal observation. Amplitude estimation of the signal can be affected by temporal variations in background noise. For example, the signal amplitude can be overestimated when we use the envelope of the signal in the time domain. In the usual operations of the ACROSS source, the noise level and the signal from the source are continuously monitored, which is achieved by a frequency modulation (FM) of the source frequency at a precisely constant time interval. With this operation, the signal has a discrete series in the frequency domain, and the background noise appears between the signal peaks.
In this study, we demonstrate a method to estimate the temporal variation in the amplitude by removing the noise variation effect, which is independently obtained during the observation of the ACROSS source. Additionally, we demonstrate the effectiveness of our method by applying it to the dataset obtained in a previous study by Ikuta et al. (2002), in which coseismic change in the seismic velocity was clearly visible. Ikuta et al. (2002) interpreted that the changes were caused by the self-opening of the preferred oriented, fluid-filled cracks due to pore pressure increase. Their interpretation is that the changes in seismic attenuation can be caused by viscous losses of fluid in the cracks.

Estimation of noise level
Our method fully utilized the unique characteristics of the signal generated by the ACROSS source. In general, ACROSS sources are operated with FM at a precise time period in repeating intervals. Thus, signal peaks of ACROSS appear with an interval of the reciprocal of the FM period in the frequency domain. In this operation, we can obtain spectral components that include only noise between all the signal components of ACROSS when Fourier transformation is applied to the data whose length is an integer multiple of the FM period. For example, when an identical modulation repeats every 10 s and Fourier transformation is applied on 50 s long data, the sampling interval in the frequency domain is 1/50 Hz, and the signal of the ACROSS source appears at every five samples (Fig. 1). We refer to the spectral components including the ACROSS signal as the "signal channel", and other components that include only background noise as the "noise channel".

Fig. 1
Characteristics of the ACROSS signal in frequency domain. Typical distributions of the signal channels and noise channels in the frequency domain are illustrated. The signal channels appear at a constant interval, and noise channels appear in between. The intervals for the signal and noise channels are reciprocal of the modulation period and the data length, respectively As the background noise is also included in the signal channels of the ACROSS, the noise causes a bias result in the estimation of the ACROSS signal amplitude. Therefore, it is necessary to make an adequate estimation of the noise that overlaps the signal, especially to estimate the temporal variation in the signal amplitude because the level of the background noise varies with time. The noise in a signal channel can be estimated from the neighboring noise channels under the assumption that the noise in both channels has the same stochastic nature.

Method to estimate the amplitude variation
We estimated a temporal change in the amplitudes of a target phase (e.g., P or S wave). The signal and noise variance for the target phases were estimated based on the method of Ikuta et al. (2002). In their study, the transfer function in the frequency domain G k corresponding to the target phase was calculated as the following convolution: where X j and S j are the spectral components of the signal channel and the source signal, respectively, H k−j is the Fourier transform of the Hanning window function, and N is the number of signal channels used in the analysis. Similarly, we can calculate the noise variance corresponding to G k as where ε 2 j is the variance of the noise in the signal channel, and calculated from the neighboring noise channels as where X m is the spectral component of the noise channels neighboring the corresponding signal channel, * denotes the complex conjugate, and M is the number of noise channels used. We express the temporal variations of the amplitude based on the total power ratio of the ACROSS signal at each time to those of reference transfer function. As the signal channel includes both signal and noise, the power of the signal can be estimated by subtracting the noise power from the power of the signal channel. Therefore, the ratio r can be calculated using the following equation: where G k0 is the reference transfer function obtained using Eq. (1). Summation over all available components of the ACROSS signal are calculated.
The uncertainty of the ratio is given by where α 2 and β 2 are the variances of the real and imaginary parts of the noise, respectively, and η k denotes the estimation error of σ 2 k . See the Appendix for the derivation of Eq. (5).

Performance test with a synthetic data
To demonstrate the advantage of our method, we conducted a performance test with a synthetic dataset of (3) Fig. 2 Result of the performance test. Amplitudes estimated using our method and the envelope calculation are shown by the dots in the left and right panels, respectively. Horizontal axes indicate the amplitude of noise relative to the root mean square of the signal. Vertical axes indicate the amplitude estimated with each method. The bold line shows the average of the calculation for each noise amplitude. The gray dashed line in the left panel indicates 95% confidence intervals, which are theoretically estimated from the noise amplitude signal and noise. By using the synthetic data, we estimated the amplitude of the signal using our method and compared the result with that of the envelope calculation. The signal was created as a series of spectral peaks with constant amplitudes and random phases in a certain frequency range. The signal was converted to a time series with inverse Fourier transform, and the synthetic data were obtained by applying a Hanning window. We produced frequency series with an increment of 0.01 Hz, and inserted spectral peaks between 10 and 20 Hz with a 0.1 Hz interval. It means 9 noise channels between every neighboring pair of signal channels. A Hanning window between 2 and 4 s was applied on the time series. Noises that follow a normal distribution with zero mean were added to the time series.
We compared the amplitudes of the signal in the synthetic data estimated by our method to that of the envelope calculation. Standard deviations of the noise were given with a step of 0.01 from 0 to 5, which were normalized by the root mean square (RMS) of the signal. We computed 1000 times for each step to compare the theoretical estimation of error with the scattering of the data.
The results of the performance test are shown in Fig. 2. The amplitudes estimated with our method did not show systematic bias resulting from the noise. The estimated amplitude average was almost constant with variable noise amplitudes. The scattering of the amplitude estimation coincides with the theoretical error estimation. In contrast, the estimation with the envelope calculation clearly has a bias result that increases with noise amplitude.

Application to the Awaji data
The dataset used in this study was obtained in an experiment at the Awaji site in western Japan (Ikuta et al. 2002;Fig. 3 Location of ACROSS vibrators and borehole seismometers. a Location of the site where ACROSS was deployed (yellow star) and epicenter of two earthquakes, which are mentioned in the text (red stars). The area of the map (a) is indicated by a black square in map (b). c Location of the ACROSS vibrator and borehole seismometers on a vertical cross section. The yellow star indicates the ACROSS vibrator. The red and blue circles show the 800-m and 1700-m sensors, respectively Ikuta and Yamaoka 2004), where a surface fault appeared at the 1995 Hyogo-ken Nanbu (Kobe) Earthquake (M jma = 7.3) (Nakata and Yomogida 1995). At this site, an experiment with two ACROSS vibrators was conducted for 15 months from January 2000 to April 2001. They monitored changes in travel time for P and S waves and found coseismic delays in two earthquakes that caused relatively large ground motion at the site (Ikuta et al. 2002;Ikuta and Yamaoka 2004).
The two ACROSS vibrators, which were deployed at the surface of granite ground, generated elastic waves with FM in different frequency ranges covering between 10 to 23 Hz. As the FM period was 5 s, the sources produced discrete signals with an interval of 0.2 Hz in the frequency domain. Unlike recent experiments with ACROSS sources, no rotation switching had been implemented at that time, that is, the weight rotated only clockwise from the top throughout the experiment.
Seismic waves excited by the ACROSS vibrators were recorded at the bottoms of the two boreholes which were drilled into granite basement in the experiment site. The seismometers in the boreholes are three-component velocity-type sensors with a natural frequency of 3 Hz and 2 Hz for vertical and horizontal components in the 800-m-deep borehole, and 4.5 Hz for the 1700-m-deep borehole. Thus, frequency response of the seismometers can be assumed to be flat between 10 to 23 Hz, which is the frequency range used by the ACROSS sources. Both sensors were located nearly vertically downward from the sources with small horizontal distances of 50 m and 120 m from the vibrators (Fig. 3). The water head in the 800-m borehole is higher than the surface (Kitagawa and Kano 2016).
We estimated a temporal change in attenuation using a series of transfer functions that were obtained in Ikuta et al. (2002). In this study, we performed an independent data screening to remove the observation period in which ACROSS vibrators are not in proper operation. The first step of the data screening was conducted based on the operation logs of the ACROSS source. The rotation frequency and mass position were logged in at 1-s intervals and were used to eliminate the operation stop period. The cross covariance of the transfer functions for each 1-h interval with the reference time is also used for data screening in the second step. We removed the data for which the cross covariances were less than 0.75. The low cross covariance may imply flaws in the seismic observation systems.
After screening, the continuous wave records were divided into 100-s segments that were stacked together in hourly intervals. Stacking was performed with a weighted average adopting the inverse maximum amplitude in each segment as the weight. The stacked waveforms were deconvolved with the source function of the ACROSS vibrators' force in the frequency domain to obtain transfer functions. As we used 100-s-long data for stacking, its Fourier transform contained spectral components with an interval of 0.01 Hz. As the period of FM was 5 s, the source signal provided spectral components with an interval of 0.2 Hz. This caused the ACROSS signal channel to appear with every 20 components in the frequency series for each vibration source.
We selected P and S waves in the transfer function in the time domain, which is calculated by the inverse Fourier transformation of the signal channel. P and S waves were selected for the 800-m and 1700-m boreholes in the same way as Ikuta et al. (2002) and Ikuta and Yamaoka (2004) for comparison. The waveforms that included either P or S waves were selected with a 20% Hanning window of 0.3 s centered at the travel time of P or S waves The area with the red area shows the P and S phases selected for calculating the change in amplitude and travel time. We selected P waves from the UD component (P-UD) and S waves from two horizontal components (S-H1 and S-H2) estimated from the typical velocity of the country rock. P waves were selected from the vertical components, and S waves were selected from the horizontal components (Fig. 4). The transfer functions in the time domain, shown in Fig. 4, were obtained by inverse Fourier transformation of those in the frequency domain. P, S, and several later phases were identified in the transfer functions for the seismometers in both the 800-m and 1700-m boreholes.
We applied our method to P and S waves. For S waves, we analyzed the signal of each component of the seismometers instead of synthetizing the radial and transverse components. Thus, we name two horizontal components, H1 and H2, for both sensors in the 800-m and 1700-m boreholes. The H1 and H2 components in the 800-m borehole were directed N126 • E and N216 • E, and those of the 1700-m borehole were directed N90 • E and N180 • E. The S phase that appears in each component is denoted as S-H1 and S-H2.
The estimation errors of the amplitude variation were calculated based on Eq. (5). The dataset of the experiments obtained by Ikuta et al. (2002) only includes the mean amplitude of the noise channels for each of the signal channels. Thus, we assume that the variances of the real and imaginary parts of the noise ( α 2 and β 2 in Eq. 5) are identical, and η 2 k is negligible. We calculated the estimation error of the amplitude variation using the following equation:

Temporal changes in travel time and amplitude
Temporal variations in the amplitudes and travel times estimated for P, S-H1, and S-H2 observed by the sensors in the 800-m and 1700-m boreholes are shown in Fig. 5. The hourly data show the daily and long-term, or eventrelated, variations. This study calculated the variation in travel time for each component and illustrated the hourly variation instead of the plots in Ikuta et al. (2002), which showed radial and transverse components with a moving average over 24 h. We confirmed that the variations in travel time are similar to those in Ikuta et al. (2002) (see the top three panels in Fig. 5). Long-term variation in P and S waves, such as gradual delay from January to March and gradual advance from April to May in 2000, were observed. A sudden delay in extremely heavy rain on September 11 was observed in all the components. The delay was preceded by a 2-month-long advance of travel time, which corresponds to the dry summer season at this site. This change has been interpreted by Ikuta et al. (2002) as the change in the source motion with reference to the weight rotation caused by the change in energy dissipation immediately local to the source. Coseismic delay and gradual recovery were observed at the teleseismic event on October 6, 2000. In addition, changes that may be associated with the second and fourth water injection experiments were observed for the 800-m-borehole, although Ikuta et al. (2002) did not mention this event.
Several characteristics of the amplitude variation were seen in the experimental period (see the bottom three panels in Fig. 5). The variation measured with our method reflects both intrinsic and scattering attenuation.
The overall pattern of the amplitude variation in P waves was similar for both the 800-m and 1700-m boreholes and gradually decreased from January to February, soon after the onset of the ACROSS experiment. The extremely heavy rainfall on September 11 also changed the amplitudes of the S wave after a gradual change from mid-August. The amplitude of S-H1 decreases for the 800-m borehole after a gradual increase but increases for the 1700-m borehole after a gradual decrease. However, only a few changes were observed in P and S-H2. In contrast to the travel time change, it is difficult to interpret the cause of the amplitude variation with only the source-ground interaction. A possible cause is the interference of multipath waves that suffer radiation variation near the source.
Diurnal variations that correlated well with temperature changes (Fig. 6) were observed. For the 800-m borehole, the amplitude had a positive correlation with temperature for all the phases. The P and S-H1 showed a larger variation than that of S-H2. The maximum variation in P and S-H1 reached approximately 5%, while S-H2 reached approximately 3%. For the 1700-m borehole, P and S-H1 showed a positive correlation with the temperature, whereas S-H2 showed a negative correlation. The daily variation in P (approximately 10%) is larger than that of S-H1 and S-H2 (approximately 3%). Good correlations with temperature were discussed by previous studies (Ben-Zion and Leary 1986; Richter et al. 2014;Wang et al., 2020) and were attributed to the thermal strains. Our results show no increase in delay or decay with depth, indicating the region of stiffness change is constrained near the surface. This shows that the thermal strain is a possible cause of the diurnal variation.
On October 31, 2000, a step-like decrease in the amplitude of all the phases for the 1700-m boreholes was observed. This coincides with the time when the top of the borehole was opened for temperature measurement (Yamano and Goto 2005). Interestingly, no change was observed in the travel time. This may be caused by a change in the coupling between the sensors and the surrounding rock, which may reduce only the amplitude. Maximum acceleration and maximum velocity at four earthquakes with relatively large ground motion. The maximum acceleration and velocity were calculated by averaging of maximum amplitude of acceleration and velocity waveforms at four accelerometers in the Awaji Island, respectively. The velocity waveform was calculated from acceleration using a high-pass filter with a cut-off frequency of 0.1 Hz. and D in the last lines in captions of each earthquake indicates epicenter distance and depth of the earthquake, respectively The change in the coupling of the sensor might be caused by the change in the water flow in the borehole induced by the opening of the well top.

Coseismic change
During the ACROSS experiment, four earthquakes with relatively large ground motions were observed at the site. Among them, two earthquakes, the Western Tottori Earthquake (WT) on October 6, 2000, and the Geiyo Earthquake (GY) on March 24, 2001, showed coseismic change in amplitude as well as seismic velocity. The magnitudes and epicenter distances of WT and GY were 6.6 and 165 km, and 6.4 and 215 km, respectively. Figure 6 shows the temporal changes in amplitude during the 10 days before and after the earthquakes. The amplitude decreases for both earthquakes, and the change at the WT was greater than that at the GY, which followed the same tendency as travel times. For the WT, a sudden drop of approximately 3% and 5% erwe observed in S-H2 for the 800-m borehole and S-H1 for the 1700-m borehole, respectively. This may imply that the coseismic change in attenuation is homogeneous in the part shallower than the 1700 m in this site. Gradual recoveries lasted approximately 1 week after the earthquakes were observed in P and S for both boreholes (Fig. 6a). The amplitude recovered to a similar level as before the earthquake. In contrast, for the GY, only S-H2 of the 800-m borehole showed a clear step-like decrease of 1% (Fig. 6b). Gradual recoveries were barely observed in any of the phases of the sensor.
To discuss the relationship between the ground motion at the earthquakes and the coseismic changes in the During the water injection experiment, the transfer function of the 1700-m sensor could not be analyzed because of the noise caused by the injection. The gap between horizontal lines with a pair of triangles on the right side of the panels show the mean error of the amplitude change amplitude, we calculated both peak ground acceleration (PGA) and peak ground velocity (PGV). Figure 7 shows the PGA and PGV for the four earthquakes calculated from the acceleration seismogram observed by K-net (National Research Institute for Earth Science and Disaster Resilience 2019). The PGAs in this figure are the average of the PGAs at four stations, HYG024, HYG025, HYG026, and HYG027. The PGVs are also calculated in the same way with velocity seismogram which is calculated by numerical integration of the acceleration record using a high-pass filter with a cut-off frequency of 0.1 Hz.
The effect of PGA on coseismic velocity change was already discussed by Ikuta and Yamaoka (2004). They claimed that only larger two of the four earthquakes showed delay in travel time and the earthquake with larger PGA showed larger delay. They concluded that the coseismic delay and the PGA correlate positively, even though a threshold exists. On the other hand, previous studies (Sawazaki and Snieder 2013;Brenguier et al. 2014) showed that coseismic PGV and derived maximum dynamic strain (MDS) are more effective on seismic velocity than PGA. Sawazaki and Snieder (2013) also reported a correlation between coseismic delays and MDSs in the range where the MDSs exceed a threshold.
We compare the PGA and PGV with the coseismic change in amplitude. Similar to the previous studies on seismic velocity, a larger change is observed at larger coseismic ground motion even for the amplitude change. The PGA and PGV at WT were about 1.5 and 15 times as large as those at GY, respectively. The PGA and PGV at GY were about 2 and 4 times as large as those at the other two earthquakes, respectively. The coseismic changes in amplitude at WT were very clear, but very small changes were observed at GY. This may reflect the large difference in PGV between the two earthquakes. A threshold may exist between the ground motions of GY and other two earthquakes because no coseismic change in amplitude are observed for smaller two earthquakes.

Change at water injections
In this experiment site, water injection experiments were conducted repeatedly from 1997 to monitor healing process in the Nojima Fault zone after the 1995 Kobe Earthquake (Ando 2001). Four injection operations were conducted during the ACROSS operation period. The first operation was performed on January 4, 2000 but was aborted due to the leakage of water at the well top. The second, third, and fourth operations were successfully performed on January 22-26, January 31 to February 5, and March 3-11, respectively. During these three operations, the water was injected at constant pressures of 2.9, 4.0, and 4.5 MPa, respectively. A total of 457 kL of water was injected during the four operations (Nishigami 2001). Water was pumped into the 1700-m borehole to feed the water into the bedrock. It is assumed that the water leaked to the bedrock at a depth of 580 m, which was revealed by temperature monitoring using an optical fiber thermometer (Yamano and Goto 2001).
Significant changes in the transfer function in both the amplitude and travel time were observed during two of the four injection operations (Fig. 8). Increases in amplitude were observed in P and S-H1 of the 800-m borehole in the second and fourth periods of water injection. The timing of the increases in travel times were similar to those of the increase in amplitude. However, no significant change was observed in the third period. The similarity of the variation between the amplitude and travel time with the variation in the atmospheric temperature in the second and fourth periods indicates that the changes during the periods of the water injection experiment may be due to changes in atmospheric temperature. Variation in the atmospheric temperature might change the elastic property of the surface part of the ground around the ACROSS source and affect the transfer function.
We tried to estimate the effect of atmospheric temperature on the amplitude variation (green dots in Fig. 8) using the method as explained later for Fig. 12. We assumed that the variation follows a linear function of temperature and calculated its coefficients using the data for 100 h before the beginning of each period of the injection operation. As a result, the changes were reduced especially in the second periods, suggesting the effect of atmospheric temperature. However, the changes in the fourth experiment remained, which may indicate the effect of water injection to increase the efficiency of wave transmission by pore pressure increase.

Anisotropy of the coseismic changes
Anisotropic changes in travel time were a key aspect in revealing the mechanism of coseismic change in Ikuta and Yamaoka (2004). In this chapter, we aim to analyze the anisotropic nature of coseismic change in attenuation. We reinvestigated anisotropic change in travel time and conducted a method to estimate anisotropy of the amplitude change with the ACROSS monitoring data. Ikuta and Yamaoka (2004) estimated the orientation of anisotropy by calculating the coseismic change of S wave splitting. They calculated the magnitude of the coseismic delay of S waves in six directions by rotating the oscillation direction with an interval of 30°. They estimated the direction of the principal axes of anisotropy by fitting an ellipse to the azimuthal distribution of the delay. We reanalyzed the anisotropy of coseismic changes in the same way, except for the azimuthal interval. We estimated the travel time delays with an azimuthal interval of one degree using the following procedure.
We used the transfer functions between 100 h before and 24 h after the earthquakes and synthesized transfer functions for all directions. We calculated the changes in travel time of the transfer functions using cross spectrum with reference to the average of the transfer functions in the first 24 h of the selected period. Because large variations that are quite similar to the changes in atmospheric temperature are observed, we reduced the effect of the atmospheric temperature by assuming that the variation follows a linear function of temperature. We calculated the coefficients of the linear function using the data obtained from the 100 h before the earthquakes.
The coseismic delay in each direction is modeled by the following equation: where a represents the change rate before the earthquake, H(t) denotes the Heaviside step function, b represents the recovery rate after the earthquake, c is the coseismic step,   8) and (9) for this calculation t is the time from the beginning of the modeling period, and T is the time of the earthquake. We estimated the parameters a, b and, c with a least square method. Figure 9 shows the azimuthal variation in the coseismic delay with a one degree interval. Our result coincides with the result of Ikuta and Yamaoka (2004) at their calculation azimuths. However, the azimuthal pattern of coseismic delay may not be approximated by an orbit as they did for fewer data points.
To reveal the meaning of the pattern, we performed a simulation with synthetic transfer functions. We modeled the wave that traveled vertically downward in a uniformly anisotropic media from the ACROSS source. As the S wave from the ACROSS source is expressed as two components of sinusoidal vibration with a right angle and a phase difference of π/2 , we used the direction of the two components of the transfer function as the principal axes of anisotropy. The transfer function in the principal axes in the frequency domain can be written as where ω is the angular frequency of the signal, φ(ω) is the initial phase for each angular frequency, and δt 0 is the difference in travel time due to the anisotropy of the surrounding media. We assumed that the amplitude ratio between the two functions, A and ω , is constant and provided φ(ω) as a uniformly distributed random number between 0 and 2π.
We assumed that the principal axes of anisotropy do not change, we got the transfer functions after the coseismic step as follows: where δt 1 and δt 2 are the steps in travel time in the directions of S 1 and S 2 . Next, we synthesized a transfer function with a linear combination for the direction angle θ as follows: We calculated the delay in travel time as a function of θ and found the azimuthal patterns similar to those shown in Fig. 9. We can reproduce a similar pattern by giving proper parameters in the calculation, as shown in Fig. 10. By using this model, we updated the azimuthal direction of the principal axes of anisotropy. We searched for the parameters that provide the best fit to the observed angular pattern shown in Fig. 9 using the least squares method. The parameters we searched for were δt 0 , δt 1 , δt 2 , A , and azimuthal direction of the major axis of anisotropy. The major axis denotes the principal axis showing larger delay. As the calculation has a trade-off between δt 0 and A , A is given independently with the observed transfer function. A is given as the ratio of the amplitude of the transfer functions in the azimuthal direction of two principal axes of anisotropy.
(10) S θ = S 1 cosθ + S 2 sinθ, The results are shown in Fig. 11. For WT, the directions of the major axis were N 90 • E for both the 800-m and 1700-m borehole sensors. For GY, the directions were N85° E and N95° E for the 800-m and 1700-m borehole sensors, respectively. These directions coincide well with the results of Ikuta and Yamaoka (2004), in which the direction was N100°E for WT and GY at the 800-m sensors. Now, we can estimate the anisotropic change in amplitude using the directions of the principal axes of velocity anisotropy. In our simulation, we found that the azimuthal variation in the amplitude was affected by the change in the travel time even if no change in amplitude was given to the synthetic data. However, the changes in amplitude are not affected by the travel time change for the direction of the principal axes, which is apparent from Eqs. (8)-(10). Therefore, the anisotropy of coseismic change in amplitude is evaluated using the variation in the direction of the principal axes of anisotropy.
Changes in the amplitude and result of the model fitting are shown in Fig. 12. The coseismic changes in the two principal axes are estimated using the same procedure as for the travel time analysis. We also corrected the daily variation correlated with the atmospheric temperature. To check the consistency of the correction, we showed a variation that is synthesized from the atmospheric temperature with a proper coefficient as Appendix Fig. 14. For WT, a greater decrease in amplitude is obtained in the major axis of velocity anisotropy for both depths. For GY, no significant difference was observed for either depths.
Our results show that the direction of greater attenuation corresponds to that of larger delay, which are consistent with laboratory experiments with cracked and saturated media (Tao and King 1990;Chichinina et al. 2009). Thus, the attenuation anisotropy is consistent with the mechanism of coseismic change proposed by Ikuta and Yamaoka (2004).

Rotation-invariant value of change in amplitude
The result of our simulation indicates that the estimation of amplitude change based on a single component may give artificial results, especially for S waves. As indicated by Eqs. (10) and (10a), the sum of the square of the two components with a right angle gives a solution that is independent of the azimuthal angle θ . Therefore, we may use the RMS of two horizontal components of transfer functions as azimuth-invariant values. An equation to obtain the invariant amplitude for the S wave can be derived from Eq. (4) as follows: where subscripts 1 and 2 indicate the components of H1 and H2, respectively, and G 1k0 and G 2k0 are the reference transfer functions. Figure 13 shows the change in the invariant amplitude for the entire period. As a result,  Table 1 The amplitude ratio of coseismic step The amplitude ratio of the transfer functions before and after the earthquakes in each sensor and each phase. "-" indicates that no significant change was observed Step at WT Step  we confirmed the coseismic changes at WT and GY, and some other changes associated with the heavy rain and opening of the well top.

Comparison of coseismic change in terms of Q
We tried to convert the changes in amplitude to that of Q −1 for comparison with the coseismic changes observed in other studies. When we assumed that Q is independent of frequency, the amplitude ratio γ can be written as where d is the distance of the sensors from the source, A is the amplitude at the source, ω is the angular frequency of an elastic wave, c is the velocity of an elastic wave, and Q −1 0 and Q −1 are the inverse of the quality factor before and after the earthquake, respectively. We provided Q −1 = Q −1 − Q −1 0 . We may assume that the amplitude ratio obtained in this study, before and after the earthquake, is equivalent to the amplitude ratio γ at the mean frequency used in our experiment f . Then, the change in Q −1 can be written as follows: where κ is the ratio of the amplitudes before and after the earthquake in this study. The amplitude ratio in this study is calculated by (1 + aT + c)/(1 + aT ) using a , T , and c in Eq. (7). Herein, we assumed no velocity change in the equation because the coseismic velocity change was sufficiently small (< 0.5%).
We calculated Q −1 for the P-UD component and the change in the invariant amplitude of the S wave using Eq. (13). We used V P = 4.0[km/s], V S = 2.5[km/s] for c , and 16 Hz for f as the typical velocity of the surrounding rocks and the center of the source frequency of the ACROSS operation. Table 1 shows the amplitude ratio, and Table 2 shows Q −1 calculated by Eq. (13).
These Q −1 values have a similar order of magnitude as reported in previous studies using spectral ratio methods. For example, Kelly et al. (2013) reported that Q −1 in the fault region of the 2004 Parkfield Earthquake (M W = 6.0) was of the order of 1.0 × 10 -3 . Wang and Ma (2015) found a decrease in Q S associated with the 1999 Chi-Chi Earthquake (M W = 7.6). The Q S was changed from 238 to 157, which corresponds to 2.2 × 10 -3 of Q −1 .

Conclusions
We developed a method to detect amplitude changes by utilizing the stable artificial seismic source, ACROSS, and applied it to in situ data. The amplitude change was obtained by calculating the ratio of the power of the transfer function of each period to that of the reference time. The noise that is estimated independent of the signal is used to avoid an apparent change in amplitude caused by the drift of the noise level. We applied this method to the data acquired by Ikuta et al. (2002). A coseismic decrease in amplitude was detected for two earthquakes with a maximum decrease of approximately 5%. The anisotropy of the coseismic changes showed that the direction of larger delay in travel time shows a greater decrease in amplitude. This result is consistent with the results of laboratory experiments with cracked media. We calculated Q −1 of the coseismic change corresponding to the amplitude change, and consistent values with previous studies were obtained.

Appendix 1: Derivation of error formula
In the routine procedure of the analysis of ACROSS data, signal and noise channels are separated in the frequency domain, as explained in the text. We may assume that the noise in the signal channels has the same statistical property as the noise in the noise channels. Therefore, a signal channel is composed of signal and noise as where G 0 is the signal and ε is the noise. Note that G 0 and ε are complex numbers. The noise ε follows the statistical property as where <> represents the expectation and " * " denotes complex conjugate and σ 2 0 is the variance of noise. Because the noise variance is estimated from the noise channels around the signal channel, the variance of the noise in the signal channel has an error δ . Thus, we express the noise variance as: The error of the variance δ follows the statistical property as: We calculate the expectation and variance of a component of the numerator of Eq. (4), which is the square of the signal subtracted by the variance of noise. The expectation is expressed as follows: and the variance is expressed as where α 2 and β 2 denote the variance of the real and imaginary parts of ε , respectively. (21) Because Eq. (4) gives the amplitude by calculating the sum of the signal channels, the error of the changes in the amplitude is expressed as where k denotes the index of the signal channels. We may replace G 0k and σ 2 0 k with G k and σ 2 k , respectively. The error is given by the following equation: See Fig. 14. (23) . Fig. 14 The temperature effect in the amplitude variation. The amplitude variation (blue and red dots) in major and minor axes of principal axes and fitting with variation in atmospheric temperature (green lines) are shown from 100 h before the occurrence of WT and GY. The top two rows show the results for the 800-m sensors and the next two rows show those for the 1700-m sensors. The bottom row shows variation in atmospheric temperature. The left and right panels show the results for the WT and GY earthquakes, respectively