# Temporal change in seismic wave attenuation using highly stable vibration sources – Earth, Planets and Space

Jan 14, 2022

### 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:

$$begin{array}{c}Mleft(tright)=at+Hleft(t-Tright)bullet left{bleft(t-Tright)+cright},end{array}$$

(7)

where (a) represents the change rate before the earthquake, (mathrm{H}left(tright)) denotes the Heaviside step function, (b) represents the recovery rate after the earthquake, (c) is the coseismic step, (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 (pi /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

$$begin{array}{c}left{begin{array}{c}{S}_{1}left(omega right)=mathrm{exp}left{iphi left(omega right)right}\ {S}_{2}left(omega right)=Amathrm{exp}left[ileft{phi left(omega right)+frac{pi }{2}+omega delta {t}_{0}right}right]end{array}right.,end{array}$$

(8)

where (omega) is the angular frequency of the signal, (phi (omega )) is the initial phase for each angular frequency, and (delta {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 (omega), is constant and provided (phi left(omega right)) as a uniformly distributed random number between 0 and (2pi).

We assumed that the principal axes of anisotropy do not change, we got the transfer functions after the coseismic step as follows:

$$begin{array}{c}left{begin{array}{c}stackrel{sim }{{S}_{1}}=mathrm{exp}left[ileft{phi left(omega right)right}+omega delta {t}_{1}right]\ stackrel{sim }{{S}_{2}}=Amathrm{exp}left[ileft{phi left(omega right)+frac{pi }{2}+omega left(delta {t}_{0}+delta {t}_{2}right)right}right]end{array}right.,end{array}$$

(9)

where (delta {t}_{1}) and (delta {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 (theta) as follows:

$$begin{array}{c}{S}_{theta }={S}_{1}mathrm{cos}theta +{S}_{2}mathrm{sin}theta ,end{array}$$

(10)

$$begin{array}{c}stackrel{sim }{{S}_{theta }}=stackrel{sim }{{S}_{1}}mathrm{cos}theta +stackrel{sim }{{S}_{2}}mathrm{sin}theta .end{array}$$

(10a)

We calculated the delay in travel time as a function of (theta) 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 (delta {t}_{0}), (delta {t}_{1}), (delta {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 (delta {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.

The results are shown in Fig. 11. For WT, the directions of the major axis were N (90^circ) 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 (theta). 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:

$$begin{array}{c}stackrel{sim }{{r}^{2}}=frac{{sum }_{k}{({{G}_{1}}_{k}{{G}_{1}}_{k}^{*}+{{G}_{2}}_{k}{{G}_{2}}_{k}^{*})-({{sigma }_{1}^{2}}_{k}{+{sigma }_{2}^{2}}_{k})}}{{sum }_{k}{({G}_{1}}_{k0}{{G}_{1}}_{k0}^{*}{+{G}_{2}}_{k0}{{G}_{2}}_{k0}^{*})},end{array}$$

(11)

where subscripts 1 and 2 indicate the components of H1 and H2, respectively, and ({{G}_{1}}_{k0}) and ({{G}_{2}}_{k0}) are the reference transfer functions. Figure 13 shows the change in the invariant amplitude for the entire period. As a result, 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 (gamma) can be written as

$$begin{array}{c}gamma left(fright)=frac{frac{1}{d}Amathrm{exp}left(-frac{2pi fd}{2c}{Q}^{-1}right)}{frac{1}{d}Amathrm{exp}left(-frac{2pi fd}{2c}{Q}_{0}^{-1}right)}=mathrm{exp}left(-frac{pi fd}{c}Delta {Q}^{-1}right),end{array}$$

(12)

where (d) is the distance of the sensors from the source, (A) is the amplitude at the source, (omega) is the angular frequency of an elastic wave, c is the velocity of an elastic wave, and ({Q}_{0}^{-1}) and ({Q}^{-1}) are the inverse of the quality factor before and after the earthquake, respectively. We provided (Delta {Q}^{-1}={Q}^{-1}-{Q}_{0}^{-1}).

We may assume that the amplitude ratio obtained in this study, before and after the earthquake, is equivalent to the amplitude ratio (gamma) at the mean frequency used in our experiment (overline{f }). Then, the change in (Delta {Q}^{-1}) can be written as follows:

$$begin{array}{c}Delta {Q}^{-1}=-frac{c}{pi overline{f}d}mathrm{ln }left(kappa right),end{array}$$

(13)

where (kappa) 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 (Delta {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 left[mathrm{km}/mathrm{s}right], {V}_{S}=2.5 left[mathrm{km}/mathrm{s}right]) for (c), and 16 Hz for  (overline{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 (Delta {Q}^{-1}) calculated by Eq. (13).

(Delta {Q}^{-1}) estimated using Eq. (13) with the amplitude ratio shown in Table 1, ({V}_{P}=4.0 left[mathrm{km}/mathrm{s}right], {V}_{S}=2.5 [mathrm{km}/mathrm{s}]), and  (overline{f }=16) Hz.

These (Delta {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 (Delta {Q}^{-1}) in the fault region of the 2004 Parkfield Earthquake (MW = 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 (MW = 7.6). The ({Q}_{S}) was changed from 238 to 157, which corresponds to 2.2 × 10–3 of (Delta {Q}^{-1}).