A difference in a transfer function of two horizontal components between the ground level and the borehole at KiK-net Mashiki Station Anisotropy of shear wave propagation in Kyushu area based on seismic interferometry

In this study, we evaluated the travel time of S-wave between the vertical array stations based on seismic interferometry, focusing on the difference in transfer function due to two horizontal components at the KiK-net Mashiki station (KMMH16). At that time, we surveyed the differences by back azimuth (BAZ) and the polarization direction of seismic waves. Furthermore, we expanded the survey to all KiK-net stations in the Kyushu district, to conrm whether the phenomena seen at KMMH16 is specic to this location. The result shows that the difference by the polarization direction in the travel time was larger than the difference by the BAZ. This result suggests that the difference in transfer function at KMMH16 were affected by the anisotropy of the S-wave velocity. We evaluated the leading S-wave polarization directions (LSPDs) and the strength of anisotropy (ΔV) for all KiK-net stations in the Kyushu district. The LSPDs roughly correspond to the results of previous studies. The LSPDs in the forearc area are nearly perpendicular to the crustal deformation whereas those in the back-arc area are nearly parallel to it. This characteristic is similar to one found by the previous research in the Tohoku district. We examined the change in anisotropy before and after the Kumamoto earthquake at two stations, KMMH16 and KMMH14 that are located near the source region. The changes in the LSPD and the ΔV before and after the earthquake were not notable. using


Introduction
During the 2016 Kumamoto earthquake, the damage was concentrated around the center of Mashiki town. Some studies suggest that the site ampli cation had affected to the damage of Mashiki town (e.g., Motoki et al. 2016, Cho et al. 2017. Strong ground motions were recorded at the KiK-net Mashiki station (KMMH16) not only at the ground level (GL) but also in 252-m deep borehole (BH) and were used analyses as base data to analyze the cause of the heterogeneous distribution of the wooden houses during the main shock (e.g., Yamada et al., 2017). However, at KMMH16, a transfer function of the NS component from the BH to the GL differs from that of the EW component (e.g., Motoki et al., 2016). The subsurface model obtained by PS logging and microtremor does not depend on the horizontal direction of polarization direction. Since the subsurface model cannot explain the different transfer function due to the polarization direction, it is an important issue in evaluating site ampli cation to clarify the cause of the different transfer functions due to the polarization direction.
There are two possible contributing factors as to why the transfer function relies on the horizontal direction. One is the irregularity of the boundaries of sediment layers, and the other is the anisotropy of the S-wave velocity. When there is an irregular layer boundary such as a step structure, the site ampli cation changes depending on the incident angle and the BAZ. Then, consequent transfer function of the horizontal component relies on the direction of the seismic waves. The anisotropy of the S-wave velocity of the propagating medium has been reported to be caused by micro cracks arranged in the propagation media (Crampin 1994). We should consider the mechanism of the different transfer functions due to the horizontal direction when exploring the subsurface model to evaluate the ampli cation factors.
In recent years, some researchers have used seismic interferometry to evaluate the travel time between vertical arrays (e.g. Yamada et al., 2010). It is also used to monitor the recovery process of stiffness loss in subsurface soils due to the strong nonlinearity of the surface ground (Yamada et al., 2010, Sawazaki 2017). Takagi and Okada (2012) and Nakata and Snieder (2012) used seismic interferometry to evaluate the anisotropy of S-wave velocity, but their target was limited to northeastern Japan.
In this paper, we evaluated the travel time of the S-waves between the vertical array stations based on the seismic interferometry to clarify the mechanism of the propagation characteristics of the two horizontal components of KMMH16. We focused on the differences in seismic waves by BAZ and polarization direction. Furthermore, we expanded the evaluation of the difference between two horizontal components to all KiK-net stations in the Kyushu district, to con rm whether the phenomena seen at KiK-net Mashiki is unique.

Method And Target Sites
In this paper, seismic interferometry was used to investigate the travel time of an S-wave between two stations in a vertical array. We calculated the deconvolved waveform of the S-wave records at the GL and BH, and extracted the time lapse when the amplitude of the deconvolved waveform is peak. The deconvolution is evaluated by the following equation.
where Wε (f) is the Fourier spectrum of the deconvolved waveform, u b (f) is the Fourier spectrum of the record at the BH, u s (f) is the Fourier spectrum at the GL, and ε is a coe cient for preventing the in nity to Wε (f) by dividing by 0; an asterisk means a complex conjugate. We set ε to 1% of the average of the denominator in the frequency range from 1 to 13 Hz as per the work by Nakata and Snieder (2012). In this paper, since we focused on S-wave propagation, we set the analysis time to 5 s to include the arrival of direct S-wave. The analysis start time was set to 1 second before the arrival of the rst S-wave obtained from the travel timetable provided by the Japan Meteorological Agency (Ueno et al., 2002). When we found that the selected time did not correspond to the direct S-wave from visual inspection, we adjusted the time to t the direct S-wave arrival manually.
We picked a peak time of the deconvolved waveform with a precision of 0.001 s. We tted a quadratic function to the peak point and points next to the peak, and selected the peak time of the tted function (Nakata and Snieder, 2012). We regarded the peak time as the travel time of S-wave from the BH to the GL as the work by e.g. Yamada et al. (2010). We con rmed that this evaluated time corresponded to the S-wave travel time from the PS logging model in the later section, and Supplementary File 1 shows that the accuracy of the estimation of travel time is more accurate than the sampling frequency. The deconvolved waveform was applied to a bandpass lter of 1-20 Hz, and we took the ensemble average for each earthquake record.
We targeted 77 KiK-net stations in the Kyushu district except one station, OITH07. The transfer function of the NS component is larger than that of the EW component at OITH07 regardless of frequency and there was a possibility of some unknown noises being included. We analyzed the strong motion records from the observation start time of each station to the end of March 2017. We discarded records in which the rst P wave was not observed. The sampling frequency of KiK-net changed from 200 to 100 Hz in several years after the observation start. In this paper, we resampled the records by 200 Hz with 100 Hz.
Before the deconvolution analysis, we estimated the orientation of the seismometer at the GL by correlation analysis with the seismometer at the BH (Kato et al., 2001). For the orientation of the seismometer at the BH, we referred to the estimated results published on the Hi-net website (Shiomi et al., 2003). Supplementary le 2 shows the estimated installation orientation of the seismometer at the GL. Figure 1 shows the velocity waveforms, the trajectory for each 0.12-s and the deconvolved waveform at the KMMH16. Looking at the time difference between the peak phases corresponding to the GL and BH in Figure 1 (A), the travel time of the S-wave is 0.50 s in the NS direction and 0.38 s in the EW direction, suggesting that the propagation velocities differ depending on the direction. We can also con rm in the particle motion in Figure 1 (B) that the prominent phase of the EW component propagated from (I) at the BH to (III) at the GL and that of the NS component propagated from (II) to (IV). Figure 1 (C) shows the deconvolved waveform of the GL record against the BH record for the same earthquake. The time difference, 0.38 and 0.50 s, seen in the waveform in Figure 1 (A), approximately appears as the difference between the travel times from the deconvolved waveforms. Figure 2 shows the average deconvolved waveforms among records of less than 100 Gal at the GL of KMMH16. The difference between the NS and EW components is more than 0.1 s. Figure 2 also shows the vertical S-wave travel time of the subsurface model based on microtremor exploration (Motoki et al., 2016) and the PS logging model provided by the National Research Institute for Earth Science and Disaster Resilience (http://www.kyoshin.bosai.go.jp/cgi-bin/kyoshin/db/sitedat.cgi?0+KMMH16+kik). The travel time from the microtremor model corresponded to that from the deconvolved waveform of the NS component, whereas the travel time from the PS logging model was faster than that from the deconvolved waveform of the EW component. In this paper, we performed the examination by using different earthquake data set to investigate the cause of the different travel times between the NS and EW components. We also examined the effect of the arrival and polarization directions on the different travel times shown in Figure  The difference in the travel time due to the difference in the wave's arrival direction is smaller than the difference between the NS and EW components. Although the travel time by BAZ between 90° and 135° is slightly shorter than that of other directions, the difference in the travel times between the NS and EW components is maintained. The difference in the travel times of the NS and EW components cannot be explained by the arrival direction of the wave, suggesting that the difference in the travel times along the NS and EW components is not mainly due to the irregularity of the sedimental layers.

The in uence of the polarization direction
Given that the irregularity of the sediment layers barely affects the difference in the propagation characteristics of the NS and EW components, the anisotropy of the S-wave velocity can be considered to be the main factor in the cause of the different travel times. We investigated the change in propagation velocity when rotating the polarization direction to determine whether the difference in propagation characteristics is due to anisotropy. As in the study by Nakata and Snieder (2012), deconvolution was performed by rotating the polarization direction at the GL and BH by 10°. Figure 4 (A) shows the deconvolved waveform of each polarization direction, and Figure 4 (B) shows the travel time. The deconvolved waveform changes with the polarization direction, and the travel time varies discontinuously at 50° to 60° and 140° to 150°. The discontinuous change due to the change in polarization direction suggests that the propagation medium between two points has two different velocities. The second peak is indicated by a blue dot in the deconvolved waveform near the azimuth where the travel time changes discontinuously. Looking at these waveforms in Figure 4 (A) (for example, 140° and 150°), there are two peaks indicated by red and blue dots, and the red ones were selected according to the amplitude of the phase and plotted in Figure 4 (B). This suggests that S-wave splitting occurred around that direction, and the difference in the NS and EW propagation characteristics is due to the anisotropy of the S-wave velocity.

Examination using records from the Kyushu district
We evaluated the travel time from deconvolution for all stations in the Kyushu district to investigate whether there was a difference in the propagation velocities due to the polarization direction as shown in KMMH16 in section 3.2. Figure 5 (A) shows examples of the travel time at these stations. As with KMMH16, at KMMH03 the travel time changed discontinuously along the polarization direction. At FKOH01, the travel time changed continuously. At FKOH07, the travel time changed discontinuously at one side. Figure 5 (B) shows the relationship between the travel time difference and classi cation by the transition type of the travel time difference. In the case of continuous change, it is limited to where the travel time difference is within 0.04 s. The peaks due to the two different velocities appeared continuously because the time difference shortened, the phases overlapped, and only one peak was seen. We performed numerical experiments simulating anisotropy with two velocities, and con rmed that the transition type depended on the length of the travel time difference of two horizontal components. The numerical tests are shown in Supplementary le 1.
The ratio of the shear-wave velocity anisotropy ΔV was calculated using equation (2) in the same manner as noted in the study by Crampin (1994).
(2) V fast and t fast represent the propagation velocity and the travel time of the deconvolved waveform along the leading S-wave polarization direction (LSPD) respectively. V slow and t slow are the same physical quantities as V fast and t fast but along the axis perpendicular to the LSPD. We adopted two different methods for obtaining the LSPD. If the travel time changed continuously along the polarization direction, as shown for FKOH01 in Figure 5 (A), we set the LSPD as angle θ at which the difference in travel time between θ and θ+π/2 becomes the maximum. If the travel time changes discontinuously, we divided two clusters cluster 1 and cluster 2 as shown in Figure 5 (A), and then calculated the average direction of the cluster. We set the LSPD as the average direction including the shortest travel time, indicated as cluster 2 in Figure 5 (A). Figure 6 shows the deconvolved waveforms along the LSPD and the perpendicular axis to the LSPD at each observation station. The dots plotted in the waveform represent the travel times, and the time on the horizontal axis is adjusted to that the average of the travel times in the two directions is located at the center of each waveform. The time between the two dots in the horizontal axis indicates the strength of the anisotropy. There are no clear peaks in both components for KMMH05, but clear peaks for other observation stations. Figure 7 shows the results of the S-wave velocity anisotropy at each observation station. The direction of the colored line indicates the direction of the LSPD, and its length and color indicate the ratio of the anisotropy represented by equation (2). The threshold value of the size classi cation refers to the classi cation by Crampin (1994): ΔV of 0.045 or lower indicates intact rock, and ΔV of 0.10 or higher indicates heavily fractured rock. We displayed ΔV larger than 0.2 as an additional threshold to emphasize the stations with strong anisotropy.
As shown in Figure 7, the anisotropy of the S-wave velocity appears at not only KMMH16 but also other stations. The anisotropy ratio of KMMH16 is larger than those of other stations except for NGSH03. Also, some areas of the Nagasaki, Miyazaki and Oita prefectures show strong anisotropy. The stations with strong anisotropy (except NGSH03) are mainly in the area where the normal fault occurs, as shown by Matsumoto et al. (2015). The distribution of LSPD seems to roughly correspond to the results of the Kyushu region by Ishise and Oda (2009), although the depth corresponding to the results of Ishise and Oda (2009) is thought to be several tens of kilometers.
Similar to Nakata and Snieder (2012), we compared the direction of crustal deformation with the direction of anisotropy. The black arrows in Figure 7 indicate the direction and horizontal velocity of the crustal deformation by the GNSS (Headquarters for Earthquake Research Prevention [HERP], 2013). Since it is thought that the northern part of Kyushu has a complicated plate deformation, we focused on the southern part of Kyushu for comparison. The area indicated by the gray dashed line in Figure 7 (almost the forearc side) and the area indicated by the solid gray line (back arc side) have different characteristics. Although there are some exceptions, there are many stations on the forearc side where the LSPD is almost perpendicular to the direction of the plate deformation. On the back-arc side, there are many stations where the LSPD is nearly parallel to the direction of the plate deformation. We measured the direction angle of the deformation of the GNSS station with Figure 7 in order to quantitatively compare with the angles of the LSPDs. The comparisons between the back-arc side and forearc side are shown in Figure 8. It is con rmed that the LSPDs were almost perpendicular to the direction of the plate deformation on the forearc area, and parallel to those on the back-arc area. This feature seen in the LSPD was pointed out by Nakajima and Hasegawa (2008) in the Tohoku area. Although the target depth in this paper is less than 500 m, which does not correspond to the target depth in Nakajima and Hasegawa (2008), the anisotropic characteristics show similar results.
Regarding the anisotropy of the subsurface structure deeper than the BH, the waveforms at the KMMH16 BH in Figure 1 show that the S-wave of the EW component arrived earlier than that of the NS component. The LSPD from the hypocenter to the BH station shows the same trend as that from the BH to the GL. Also, Crampin and Chastin (2003) indicated that the APE model of the propagation medium due to stress can be applied to below the critical depth of 500-1000 m that is vertical stress equal to minimum horizontal stress. Above this depth, crack distribution was controlled by stress release and lithologic phenomena (Crampin and Chastin, 2003). We considered that this disturbance near the GL might be not so much as to rearrange the original crack distribution because the LSPD near the GL corresponds to the LSPD deeper than the BH.
The ΔV at KMMH16 is 0.24, which is the second largest among the 77 observation stations targeted here. We compared the ΔV at KMMH16 with those reported by previous researches. The anisotropy of the Swave velocity in northeast Japan calculated by Nakata and Snieder (2012) was 0.07 for NIGH13.
According to the results from a similar analysis of all KiK-net stations in the Tohoku area, the maximum anisotropy is 0.20 at MYGH20 (Takagi and Okada, 2012;Mizuno and Nakai, 2005). The anisotropy around the Tanna fault, based on VSP by Nakao et al. (1994), is estimated to be about 20%, and they reported the anisotropy was strong. Of the 54 stations reported on by Crampin (1994), only one site had a greater anisotropy than KMMH16. Compared with these values, it can be said that KMMH16 exhibits strong anisotropy. Crampin (1994) con rmed that anisotropy increases near the fault zone and in the area of volcanic rocks associated with high heat ow. KMMH16 is near the fault and close to Mt. Aso., and KMMH16 meets several conditions for strong anisotropy.
We con rmed whether the travel time obtained by the seismic interferometry corresponds to the travel time of the S-wave velocity, comparing the average velocity from the observed travel time with that from the PS logging models. Figure 9 shows the comparisons between the observed average S-wave velocity and the calculated one by the PS logging models and the vertical axes of the left and right panels show the propagation velocity along the LSPD and the LSPD+π/2, respectively. The observed average propagation velocity of them corresponds well with the average S-wave velocity of the PS logging. We can conclude that the travel time obtained from the deconvolution represents the propagation of the Swave in average. Figure 10 shows the relationship between the velocity from the BH to the GL along the LSPD and the ratio of the anisotropy. Although there are ve exceptionally large stations, it can be said from these results that as the average velocity along the LSPD increases, the anisotropy ratio also increases. We considered that in sediment layers with low velocity, vertical cracks seem to be only a small few cracks formed in deposition process. This indicates that anisotropy tends to appear in a layer where the S-wave velocity is high.

Characteristics Of S-wave Anisotropy
We referred to the rock classi cation based on the geological map (Wakita et al., 2009) to con rm whether the difference in anisotropy is due to the difference in geology. For example, the rock divisions of KMMH16 and KMMH03, which have strong anisotropy, are non-alkaline pyroclastic ows. Here, the average of the anisotropy ratio, the mean ΔV, is calculated for 25 stations of igneous rock, and we obtained mean a ΔV of 0.094 with a standard deviation of 0.06. The average ΔV of all the stations is 0.086, and the standard deviation is 0.06. We used the t-test for statistical processing. Because of the tvalue to be 0.58, the difference in anisotropy from igneous rocks was judged to be insigni cant according to the used threshold of 0.05 used in engineering. We considered that the difference in the stress depends on the location and is larger than the rock classi cation.
Next, we investigated the anisotropy change before and after the 2016 Kumamoto earthquake. Figure 11 shows the transition of the travel time, the LSPD and the ΔV (the ratio of the anisotropy) for KMMH16 and KMMH14, which are located near the fault of the 2016 Kumamoto earthquake. The length of each term was set so that the term after the earthquake was approximately logarithmic interval. This expression comes from the ndings that the recovery process in S-wave is proportional to the logarithmic axis of time by Sawazaki (2017). At both stations, the travel time before the foreshock was shorter than before the main shock, and also longer after the main shock. Although the travel time gradually approached the state before the earthquake, it still remained longer than the travel time before the 2016 Kumamoto earthquake. The direction of the LSPD was stable before and after the earthquake, and the ratio of anisotropy did not change noticeably. The pre-and post-earthquake stress changes around these stations (e.g., Goto et al., 2017) are estimated to be large.
Although the travel time was delayed after the earthquake, the anisotropy hardly changed in both direction and strength. Cao et al. (2019) reported that the anisotropy changed near the epicenter of the 2004 Niigata-ken Chuetsu earthquake. Takagi and Okada (2012) also reported that the anisotropy change of FKSH14 after the 2011 off the Paci c coast of Tohoku Earthquake had occurred. The difference between our target stations (KMMH16 and KMMH14) and the stations noted above was the anisotropy strength. The station with an anisotropy change has weak anisotropy with almost less than 0.045 that means pre-fracturing (Crampin, 1994) at FKSH14 with a ΔV of 0.0176 and N.YNTH with ΔV of 0.0538 (Mizuno and Nakai, 2005). Since there might be many large cracks at the observation stations with strong anisotropy, such as KMMH16 and KMMH14, the effect of changing the stress eld at a depth of several kilometers or deeper seems to be small, if the changing the stress eld does not affect the arrange of the cracks. Therefore, the effects of the earthquakes in changing anisotropy are small compared with the observation stations with weak anisotropy indicated by Cao et al. (2019) and Takagi et al. (2012).
Finally, we con rmed the effect of the anisotropy on the difference in the transfer function of the horizontal directions. Figure 12 shows examples of the comparison of the transfer function between the LSPD and the axis perpendicular to the LSPD. At stations with weak anisotropy shown in Figure 12 (A), both transfer functions are similar. On the contrary, at stations with strong anisotropy shown in Figure 12 (B), the shape of the transfer function greatly differs depending on the horizontal direction. When evaluating ground motions at observation stations with strong anisotropy, it is necessary to pay attention to the reproducibility of the ground motions. The reproducibility may be reduced unless the transfer function is considered separately in relation to the horizontal direction.

Conclusion
We focused on the difference in the horizontal component of the transfer function of KMMH16 concerning the horizontal direction. We investigated the travel time of the S-wave between the stations in the vertical array based on seismic interferometry, from three viewpoints: the effects of irregularity of the sediment layer boundaries, the anisotropy of S-wave velocity, and the ground nonlinearity during strong earthquakes.
In the deconvolved waveforms at KMMH16, the difference in the travel times among the epicenter locations was smaller than that between the NS and EW components. Based on the transition of the travel times of the deconvolved waveforms at KMMH16, the difference in the travel times depend on the polarization direction. It can be considered that the medium has two different S-wave velocities. From the results noted above, it was inferred that the difference in the propagation characteristics at KMMH16 was mainly due to the anisotropy of the S-wave velocity.
Regarding the change before and after the earthquake at KMMH16 and KMMH14 near the source region of the 2016 Kumamoto earthquake, a change in travel times was observed, but the change in anisotropy was small. This suggests that the change in the travel times before and after the earthquake was mainly due to the nonlinear behavior of the subsurface media during strong motions.
According to the observations of 77 stations in the Kyushu district, there were some stations showing two velocities. The LSPDs observed in southern Kyushu were almost perpendicular to the direction of crustal deformation on the forearc side and parallel to those on the back-arc side. These characteristics were similar to those of northeast Japan, indicated by Nakajima and Hasegawa (2008).
The transfer function between the stations in the vertical array differs depending on the orientation of the horizontal component at stations with strong anisotropy. It suggests that an evaluation of site ampli cation using a single velocity model may reduce the reproducibility of the observed ground motions. Figure 1 (A) Velocity waveforms integrated from the records at KMMH16, (B) particle motions in velocity, (C) deconvolved waveforms of the surface records with the borehole records The shaded boxes in (B) mean the time swinging started in the NS or EW direction. The dots represent travel times of deconvolved waveforms.

Figure 2
Deconvolved waveforms averaged with records less than 100 Gal at the GL of KMMH16. Comparison of travel times from deconvolved waveforms with the vertical S-wave travel time from the PS logging model from the KiK-net website and the microtremor model by Motoki et al. (2016).    Deconvolved waveforms of the LSPD with black lines and the axis perpendicular to the LSPD with gray lines. Dots represent the travel times of each waveform. The time on the horizontal axis is adjusted to that the average of the travel times in the two directions is located at the center of each waveform.    The transition of the travel times, LSPDs, and ΔV (strength of anisotropy) before and after the 2016 Kumamoto earthquake at KMMH16 and KMHH14.