Anisotropic Velocity Structure Beneath Shikoku, Japan: Insights From Receiver Function and Shear Wave Splitting Analyses

A combination of receiver function (RF) analysis and shear wave splitting can reveal the anisotropic properties of the Earth's velocity structure. The RF profile exhibits harmonic patterns partially generated from the discontinuities of anisotropic media near the receiver. The anisotropic properties are determined using splitting parameters, which include the fast polarization direction (FPD) and split time. We utilized the seismograms recorded by stations on Shikoku Island, Japan. Also, we applied the Bayesian information criterion to control the spike count of each RF, which consists of a series of spikes generated by time‐domain iterative deconvolution. Anisotropic differences arose between northern and southern regions along the 30‐km iso‐depth line of the top surface of the Philippine Sea Plate. In the southern part of Shikoku Island, the FPDs are sub‐perpendicular to the plunge of the subducting slab, while in the northern part, they are sub‐parallel to the plunge of the slab. Our results indicate that anisotropic strengths are weaker in the tectonic‐tremor band in northwestern Shikoku and stronger in the northern part of central and eastern Shikoku around the no tectonic‐tremor area. These anisotropic variations may characterize Shikoku Island’s geological structure.

• Fine-resolution anisotropic properties have been calculated above slab Moho in Shikoku Island, southwest Japan, using receiver function analysis • Weaker anisotropic strength was observed in tectonic-tremor-clustered areas with serpentinized mantle wedge, especially in northwestern Shikoku • Serpentinization of the mantle and metasomatic reactions lead to hydrofracturing and relatively low anisotropic strength

Supporting Information:
Supporting Information may be found in the online version of this article.
Plate subduction, particularly in the mantle wedge (MT), a wedge-shaped region between the subducting oceanic plate and the Mohorovičić discontinuity (Moho) of the overlying continental plate, can significantly influence seismic anisotropy.This region plays a crucial role in subduction zone dynamics, as well as in the generation of earthquakes and volcanoes.Olivine and serpentine are the primary minerals in the MT (Lassak et al., 2006;Shiomi & Park, 2008;Wang et al., 2019).Olivine is widely recognized as the primary driver of anisotropy (Silver & Chan, 1991).The contact between olivine and fluid under an appropriate thermal pressure condition can lead to serpentinization.The serpentine in the MT is commonly observed in the subduction system.The seismic azimuthal anisotropy is expected to be trench-parallel and trench-normal above steep and flatter subducting plates, respectively.This is because the mineral's geometry is projected onto the horizontal plane differently in each case (Katayama et al., 2009).Furthermore, the intricacies of the subduction zone, along the locally serpentinized MT on the plate interface, can give rise to spatial variations in the activity of slow earthquakes potentially influenced by pore pressure buildup (Katayama et al., 2013).
The Philippine Sea (PHS) Plate is subducted under the Eurasian Plate at a velocity of 4-5 cm yr −1 toward the northwestern direction (Seno et al., 1993).Previous studies have measured the geometry of the subducting PHS Plate beneath Shikoku Island (Hirose et al., 2008;Ito et al., 2009;Shiomi et al., 2006Shiomi et al., , 2020;;Ueno et al., 2008).In this study, the 30-km iso-depth line of the top surface of the PHS Plate is regarded as a rough dividing line between the northern and southern regions of Shikoku Island (Figure 1) (Hirose et al., 2008).Based on their results, as the plate is subducted, the PHS Plate dips gently at approximately 6° from the Nankai trough to central Shikoku, and the dip angle steepens toward the north (Ito et al., 2009).In western Shikoku, the PHS Plate maintains a dipping angle of approximately 8-10° and is in contact with the continental lower crust beneath the middle and southern portions of Shikoku.By contrast, the PHS Plate subducts beneath the MT in the northern area (e.g., Nakanishi et al., 2018) and the most active slow-earthquake areas in the world (e.g., Obara, 2002;Obara & Kato, 2016).Additionally, the Median Tectonic Line (MTL), the longest fault system in Japan, runs through northern Shikoku and displays distinctive structural and stress characteristics.Several previous studies have investigated the anisotropic properties beneath Shikoku Island.Anisotropy tomography of S-waves in the crust revealed that the fast polarization directions (FPDs) are sub-parallel to the trajectories of maximum horizontal compression in the NW-SE direction (Ishise & Oda, 2009).This suggests that the FPDs revealed from Ps phases on receiver functions (RFs) reflect the anisotropic properties in the upper crust and are aligned with the direction of S hMax (Nagaya et al., 2011).
The activity of tectonic tremors beneath the northern part of Shikoku Island exhibits lateral variation, particularly in the accumulated total energy over a decade (Annoura et al., 2016).The area with the highest energy release due to tectonic tremors is in the western section of Shikoku Island.In contrast, it is comparatively lower in the eastern end and middle of the island.Although tectonic tremors have shown the lateral variation on Shikoku Island, the seismic anisotropic structures of the entire island, especially in the MT, have yet to be thoroughly investigated and compared.As the fast axes of splitting shear waves reflect structural properties, such as EDA and fractured fault system in the crust and LPO in the mantle, we investigated the anisotropic properties and their spatial variation throughout Shikoku Island, considering the geological structure.In particular, we adopted the simultaneous analysis of RFs and shear wave splitting to estimate the seismic anisotropic structure and compare the lateral variation of slow earthquake activity on Shikoku Island.

Seismic Record Data Set
We utilized seismic records obtained from a high-sensitivity seismograph network (Hi-net) and a broadband seismograph network (F-net) (National Research Institute for Earth Science and Disaster Resilience, 2019a, 2019b) (Figure 2).An average of 166 seismograms at each station were selected from the U.S. Geological Survey earthquake catalog (USGS, 2022) and filtered at 0.1-20 Hz.Only high signal-to-noise ratio (≥10) teleseismic records with magnitudes greater than M w 5.8 for epicenters between 30° and 90° were used for RF calculation.We corrected the horizontal sensor orientations of the Hi-net seismograms (Shiomi, 2013) and rotated the seismograms into the LQT-coordinate system.Based on the power spectral density of each cut component calculated for both the signal and noise windows (seismic record before the P-phase arrival), we applied a Gaussian low-pass filter with a 3-dB cutoff at 1 Hz using a Gaussian factor of 5.5.

Incorporation of the Bayesian Information Criterion (BIC) Into the Time-Domain Iterative Deconvolution Method
According to a study by Ligorría and Ammon (1999), RFs calculated using the time-domain iterative deconvolution method typically exhibit numerous spikes, with each spike representing a converted Ps phase.However, the iterative procedure typically persists in generating spikes until a predefined condition is satisfied, which relies on the variance between adjacent squared errors (δε).To better control the size of the resulting model (i.e., the number of spikes), we propose a new conditional judgment based on the Bayesian information criterion (BIC).Schwarz (1978) derived BIC as where L indicates the maximized value of the likelihood function for the estimated model.We also assume that the error or disturbance of data follows a normal distribution; K is the number of estimable parameters in the approximating model, which means spike number in the iterative deconvolution RF method; n is the number of samples.The maximized value of the likelihood function L can be expressed as the averaged summation of the squared residuals,    2 .Then, Equation 1 can be written as where ε i is the estimated residual from the fitted model at each point.Then the "suitable" RF will have the minimum BIC value under the assumption that the error or disturbance of data follows normal distribution.
Silver and Chan (1991) also proposed a method to estimate the uncertainty of the grid search results.It is assumed that the signal is combined with Gaussian white noise and obeys the F distribution, as follows: Here, E t indicates the summation of the squared T-component energy and k = 2 indicates two free parameters.A quarter of the width and length of the confidence interval indicates the standard error.The confidence interval for the FPD (φ) and the delay times (δt) of the split shear waves is defined as Here, α = 0.05 decides the confidence level (95% confidence interval).In particular, v indicates the degree of freedom, and we adopted the corrected formulae rederived by Walsh et al. (2013) for its calculation of grid search uncertainty in our results.

Shear Wave Polarization Anisotropy From Receiver Function
The RF can extract the Ps-waves converted at the impedance discontinuities beneath a seismic station.The Ps-waves can constrain anisotropic properties only above the discontinuities, whereas the splitting properties of teleseismic S-waves or SKS-waves also reflect deeper structures, such as the mantle.The splitting patterns from Ps-waves are more reliable for constraining anisotropic structures above the impedance discontinuities, such as the MT and continental crust (CC).
S-wave anisotropy on RFs affects the Q-and T-components of the RFs due to the azimuthal dependence of incidence waves, making a different pattern.Regarding the back-azimuthal variation, mainly the four-and two-lobed patterns are shown in the RF profile.These are caused by the horizontal anisotropic symmetry and tilt anisotropic symmetry axes or the dipping impedance discontinuity, respectively.
Due to the intricate nature of the converted S-wave pattern in RFs, determining both the FPD and δt from a single RF trace may introduce bias.Here, we neglected the difference in slowness of each trace and assumed all RFs to be generated from the same anisotropic media.Consequently, we gathered all RF traces over the back-azimuth and analyzed them simultaneously, considering the combined pattern resulting from all factors.Like previous studies, the time window was selected through a manual visual check of Q-component RFs to encompass as many specific phases as possible (e.g., Silver & Chan, 1991).We formulated an objective function similar to the traditional grid search employed in previous studies to seek the "best" parameter pair of (φ, δt) that minimizes the total T-component energy of all traces.
To validate the performance of the combined method, we conducted a straightforward synthetic test that simulates an anisotropic medium with a dipping interface (Table 1 and Figures 3 and 4 Information S1).Following the splitting analysis, the anisotropic patterns generated by a tilted interface were successfully eliminated (Figure 3b).The synthetic results validate the effectiveness of this approach in eliminating the energy pattern resulting from dipping interfaces.Specifically, we found this approach results in standard errors of less than 5.25° for FPD and 0.1 s for delay time.Additionally, the synthetic test yields slightly biased results in φ and δt, but well recovers the modeled anisotropic parameters with a limited back-azimuth coverage in real data scenarios under the dipping angle of 8-12° (see Figure S2 in Supporting Information S1).

BIC in Receiver Function Calculation
The original time-domain iterative deconvolution method required users to define a maximum iteration number and δε as a conditional criterion.In this study, we applied BIC as a criterion of equal significance, along with δε. Figure 5 illustrates the variation curves of the sum of square error (SSE) and BIC values during the deconvolution process using observed data.The SSE value decreased as the iteration number increased, eventually satisfying the conditional criterion and stopping the iteration.In contrast, the BIC value reached its minimum at an earlier iteration stage.Incorporating BIC reduced the number of spikes in both Q-and T-components by one-third to one-half compared to the traditional conditional criterion, SSE.These results demonstrate the effectiveness of our approach in controlling the resulting model size while maintaining its reliability.We also observed phases in T-component RFs; however, these phases did not show a rigorous four-lobed pattern like an anisotropic horizontal layer.We found polarity reversals in the T-component RF phases, and the corresponding change in the polarization was also identified in the Q-component (Figure 6a).The polarity reversals in the T-component RF appeared with a back-azimuth interval of about 90°, which suggests a subtle four-lobed pattern.This back-azimuthal pattern reflects the combined effects of anisotropy and the dipping slab Moho.Clear positive phases from the slab Moho discontinuity in the Q-component of the RFs appeared at around 4-5 s (Figures 6b-6d; see Table S1 for the synthetic velocity models of those stations).

Receiver Functions
We observed a series of coherent significant negative phases slightly earlier than the positive converted phases from the slab Moho in some stations above the tectonic-tremor-clustered area (Figure 7).This may suggest that the structure beneath the stations north of the 30-km iso-depth line is considered to be the area where the oceanic crust (OC) transitions from contact to separation with the CC (Kimura et al., 2019).Consequently, we expected that a MT would result in a tight positive-negative-positive phase pattern in the RF, as shown in Figure 7.

Anisotropic Features Revealed by Splitting Analysis
Based on the anisotropic multi-layer assumption, we identified distinct, consistently converted phases at the interface, encompassing all seismic anisotropic effects above it.We focused on the anisotropic effects above the slab Moho: CC, MT, and OC, as shown in Figure 8.The calculation of FPDs was conducted individually for each station.These FPDs and delay times were estimated with an acceptable standard error of 2σ ≤ 20° (Shiomi et al., 2020).
The results of the FPDs showed different spatial characteristics between the southern and northern regions.In the south, the FPDs show a NE-SW direction; in the north, they tend to rotate in the direction of the slab plunge.We found that the azimuthal anisotropic fast axes changed from a trench-parallel to a trench-normal direction while moving from the southern to the northern area.This change exhibited a strong correlation with the dehydration and mineral evolution resulting from the subduction of the young and hot PHS Plate (Horn et al., 2020;Ji & Yoshioka, 2017;Lee & Kim, 2021).The anisotropic strength was relatively weaker along the tectonic-tremor-clustered area, particularly in the northwestern part of Shikoku Island, with an average split time of 0.057 s.

Characteristics of Receiver Function
Our observations of the slab Moho-converted phases in eastern Shikoku Island indicate that these phases occur, on average, at a delay time of approximately 4-5 s.Specifically, in the area north of the 30-km iso-depth line in eastern Shikoku, the arrival time of phases is earlier by approximately 0.5-1 s compared to that in western Shikoku.Ito et al. (2009) reported that the slab in eastern Shikoku changes its dip angle from 6° to a steeper angle further north of the MTL.We calculated the average arrival time of slab Moho-converted phases at each station in Shikoku Island (Figure S4 in Supporting Information S1).Since each Hi-net sensor is installed in boreholes at depths of approximately 100-200 m, we have excluded the effects of the low-velocity layer (average V s < 1.4 km/s) from our discussion with considering the quarter-wavelength resolution.
We hypothesize that the relatively near arrival time of the slab Moho-converted phases in this area may arise from the variations in the thickness of sedimentary layers (V p < 4.8 km/s) and metamorphic rock layers (V p ≈ 5.7 km/s), which become thinner toward the north.This hypothesis is consistent with the findings of a previous survey conducted by the Cabinet Office of Japan in 2003.The opposing influences of a later time lag due to a deeper slab Moho and an earlier time lag caused by a thinner sedimentary layer may counterbalance one another, resulting in nearly identical arrival times for the slab Moho Ps phases.
We assumed several synthetic velocity models for typical stations and corresponding synthetic RF profiles (see Table S1 and Figure S5 in Supporting Information S1) by referring to previous studies (Cabinet Office Japan, Disaster Management, 2003;Ito et al., 2009;Kimura, 1979;Shiomi et al., 2020).Our observations of the Moho slab in western Shikoku are consistent with previous studies (Nagaya et al., 2011;Shiomi et al., 2020).Regarding the observations in eastern Shikoku at a depth of approximately 30 km, we considered that the converted Ps phases from slab Moho might exhibit relatively consistent characters in the Q-component.Notably, the delay time of phases in RF depends on several factors, such as interface depth and V p /V s .However, in eastern Shikoku, when compared to western Shikoku, the variation in sedimentary thickness makes it challenging to determine the actual depth of the Moho slab.
Seismic multiples can significantly influence RFs, impacting methods such as the H-κ stacking used to determine crustal thickness (Zhu & Kanamori, 2000).Furthermore, reverberations from sedimentary layers can also affect RF analysis (Cunningham & Lekic, 2019;Yu et al., 2015;Zhang & Olugboji, 2021).A previous survey in Shikoku Island indicated that the top interface of the layer with a shear wave velocity of ∼3 km/s is located at depths of 1.2-1.6 km and 1.6-2.5 km in the northern and southern parts, respectively (Cabinet Office Japan, Disaster Management, 2003).When reverberations from sedimentary layers are taken into account, it becomes apparent that previous studies might have potentially overestimated the anisotropic properties of Shikoku Island (Watanabe & Oda, 2015).This is probably because these studies treated the reverberations as converted phases RUAN ET AL.
10.1029/2023JB027178 9 of 13 generated from a relatively deeper interface within the CC beneath the F-net stations on Shikoku Island, consequently leading to an overestimation of the anisotropic strength in the upper crust.

Anisotropic Properties in the Shikoku Island
EDA remains a valuable tool for understanding crustal anisotropic features.Uchide et al. (2022b) investigated crustal stress orientations in Japan using focal mechanism solutions from small earthquakes.Their results for Shikoku Island show that the local crustal S hMax is oriented in the E-W to NE-SW direction to the south of the MTL.Our results show that FPDs exhibit a shift in an iso-depth direction in most of the southern part of Shikoku Island (Figure 9), potentially indicating that oriented fractures caused by bending-related faulting of the oceanic slab dominate the anisotropy direction (Faccenda et al., 2008).In addition, the FPD results are consistent with the orientations of S hMax around the MTL in eastern Shikoku Island.This observation can be attributed to the high degree of alignment in fractures and microcracks due to local stress.To the north of the MTL, the orientations of the S hMax are in the WNW-ESE direction.Meanwhile, our FPD results in this region showed a different plunging direction of the slab.The difference in direction in the northern part of Shikoku Island further suggests that the medium of OC and/or MT mainly contribute to the shear wave splitting.

Thermal and Hydrous Condition in the Northern Part
The depth of the top surface of the PHS Plate corresponded to approximately 30-40 km in the northern area of Shikoku Island (Figure 8).Within the OC, the thermal gradient becomes larger at depths greater than 30 km, and the temperature is ∼400-500°C at the target depth of the PHS Plate (Yoshioka et al., 2008).Notably, the temperature distribution approximately 2.5 km above the PHS Plate surface ranged from 450 to 600°C (Tatsumi et al., 2020).Commonly, chrysotile and lizardite tend to form at 50-300°C, transforming into antigorite at ∼300-400°C, and then stabilizing at ∼400-600°C (Evans, 2010;Lafay et al., 2013).Tectonic tremors may be associated with the punctuated dehydration of the subducting OC, where brittle deformation occurs otherwise a ductile creep regime (Condit et al., 2020).According to Mizukami et al. (2014), serpentinization and metasomatic reactions can result in a significant amount of residual water, elevating pore fluid pressure within SiO 2 -rich fluid systems.

Tectonic Tremor and Serpentinite Layer Above the Plate Interface
Differences in tectonic-tremor activity may be caused by lateral variations in pore fluid pressure and the mechanical strength of the hanging wall mantle.
Our observations reveal a strong anisotropy with a delay time of 0.2 s in the N-S direction beneath the stations where the subducting slab plunges deeper and without tectonic tremors.This sub-parallel slab plunge is commonly detected in subduction zones (Castellanos et al., 2020;Karato et al., 2008;Zhao et al., 2021).The lineation fabric (LPO) of olivine is likely generated by the flow-induced strain within the MT (Nakajima & Hasegawa, 2004) and/or a serpentinized layer on the plate interface with low-angle subduction (Shiomi et al., 2020).
In northeastern Shikoku, the strong anisotropy results from the converted phase of the slab Moho, which encompasses anisotropy contributions from the CC, OC, and MT.The layered serpentinite above the plate interface within the MT exhibits the same strength of anisotropy as the olivine in the MT, which is 10 times thicker than the serpentinite layer (Katayama et al., 2009).Generally, it seems unlikely that the relatively large delay time is solely  RUAN ET AL.
10.1029/2023JB027178 10 of 13 attributed to olivine LPO.However, it is possible that the accumulative effects of all anisotropic media could account for such a delay time, with a serpentinized layer with a thickness of approximately 1 km (Figure 10).
The observed FPDs along the iso-depth contours between a depth of 30 and 40 km, mainly around the sources of tectonic tremors, had an average delay time of 0.057 s, possibly reflecting the degree of serpentinization of the MT (Figure 8).The observed strength and direction of FPDs are consistent with those in the western Kii Peninsula (Saiga et al., 2013).Also, the change in FPDs is consistent with that of Shiomi et al. (2020).The dehydration of the OC results in the serpentinite layer on the plate interface generally growing above the plate interface around the MT (Katayama et al., 2012;Lee & Kim, 2021).Furthermore, the serpentinite layer or the serpentinized MT corner can host tectonic tremors and regulate slow-earthquake activity and plate interface coupling (e.g., Katayama et al., 2012;Tarling et al., 2019) (Figure 10).
We observed weaker anisotropy along the iso-depth contours across the tectonic-tremor band (Figure 8) in comparison to other subduction zones with serpentinized layers, such as the Ryukyu subduction zone (Katayama et al., 2009).Nakajima and Hasegawa (2016) proposed the difference in strength of anisotropy due to the drainage ability of the overlying plate in megathrust areas.One possible explanation for this weaker anisotropy within the MT is the high fluid pressure conditions around the tectonic-tremor sources.Fault weakening and shear failure caused by fluid overpressures from metasomatic reactions or serpentinization may disrupt the anisotropic features of foliate serpentine at the plate interface.This phenomenon, known as hydrofracturing, mainly occurs in regions with high tectonic-tremor excitation efficiency (Tarling et al., 2019).Furthermore, according to Tang et al. (2015), increasing pore fluid saturation causes a decrease in anisotropic parameters.The alternative perspective is the existence of a narrow MT under low-angle subduction.The tectonic-tremor sources can be concentrated near the top of the MT with a thickness of a few kilometers (e.g., Sawaki et al., 2021;Shelly et al., 2006).In this case, we expect a delay time <0.05 s in shear wave splitting due to layered serpentinite anisotropy (Katayama et al., 2009).

Conclusions
We conducted a joint analysis of RF and shear wave splitting to unveil the anisotropic features beneath Shikoku Island in southwestern Japan.In southern Shikoku, FPDs of the anisotropy are typically sub-parallel to the iso-depth contours of the top boundary of the PHS Plate, possibly reflecting the faults and the oceanic slab bending.In northern Shikoku, FPDs from the northeastern region rotate toward the slab plunge, indicating the existence of a serpentinized MT.However, in the region of the tectonic-tremor-clustered, the anisotropic strength was lower, especially in the northwestern region of the island.We did not observe a strong anisotropic strength, although serpentine is present in this region under suitable thermal-hydrous conditions.The reduced anisotropic strength in the northwestern region of Shikoku Island may be due to hydrofracturing and high pore fluid saturation, which can weaken the anisotropic strength of the antigorite alone at the slab interface.To better understand these factors and investigate the specific structure, implementing the Markov chain Monte Carlo-based inversion will help us better evaluate the model variability.

Figure 1 .
Figure 1.Epicenter distribution of teleseismic events and stations used in this study.The triangles depict the Hi-net and F-net stations in the northern and southern Shikoku Island (red and navy blue triangles, respectively).The dark red lines indicate Median Tectonic Line (MTL).The dashed lines indicate the top interface of the Philippine Sea (PHS) Plate by Hirose et al. (2008).

Figure 2 .
Figure 2. Hi-net and F-net stations in northern and southern Shikoku Island, marked using red and navy blue, respectively.

Q
-component RF traces exhibited primary phases with positive amplitudes, indicating the presence of the Moho of the subducting slab.Here, we show four example stations located approximately on the 20 and 35 km iso-depth lines in eastern and western Shikoku (Figure 6; the station's locations are shown in Figure 1; all RF profiles are shown in Figure S3 in Supporting Information S1).At station N.OOZH, located in the western part along the MTL, clear positive phases appeared at ∼5-6 s in the Q-component of the RFs, which were inferred as converted Ps phases from the slab Moho of

Figure 3 .
Figure 3. Receiver functions profile of synthetic test, (a) before and (b) after shear wave splitting correction.The phases used for the analysis were emphasized based on the time window.A two-layer model with a tilted interface and detailed assumed parameter values are listed inTable 1 with the dipping angle of 10°.

Figure 4 .
Figure 4. Grid search contour obtained from the synthetic test.The top of the panel shows the best-fit parameter pair with standard error and the time window used for the analysis.The yellow line represents the 95% confidence interval.

Figure 5 .
Figure 5. (a) Variations of the BIC and SSE values versus the iteration number.Yellow dots indicate the minimum value in BIC curves.Green dots indicate where the iteration fit the SSE condition and stopped.(b) Top and bottom panels show the Q-component receiver function traces by Bayesian information criterion (BIC) and sum of square error (SSE), respectively.

Figure 6 .
Figure 6.Examples of receiver functions (RFs) obtained at the N.OOZH, N.NAKH, N.AYKH, and N.HWSH stations of NIED Hi-net.Station locations are shown in Figure 1.The vertical and horizontal axes represent the back-azimuth and time lag after P arrival, respectively.The RFs were grouped into 5° back-azimuthal bins.The blue area indicates the positive amplitude of each RF, and the red area is the negative amplitude.The yellow lines indicate the arrival of the slab Moho-converted phases estimated from the Q-component RFs at each station.The star symbols indicate the positions of changes in phase polarization.

Figure 7 .
Figure 7. Receiver function profiles from the station above the tectonic-tremor-clustered area listed from west to east; (a) N.MISH, (b) N.UWAH, (c) N.SJOH, and (d) N.AYKH.The emphasized consistent negative phases just before the phases were converted from slab Moho indicate the existence of a relatively high-velocity layer.

Figure 8 .
Figure 8. Distribution of fast polarization directions in Shikoku district.The legend shows the splitting analysis results from different anisotropic layers in the different color bars with light-color bars indicating uncertainty.The sky-blue-colored circles indicate the distribution of tectonic tremors(Maeda & Obara, 2009;Obara et al., 2010).

Figure 9 .
Figure 9.The distributions of fast polarization direction are same as in Figure 8.The black bars indicate the direction of normalized S hMax (Uchide et al., 2022b).

Figure 10 .
Figure 10.Schematic diagram of our anisotropic observations in Shikoku Island.The anisotropic areas considered in this study were plotted in different colors.Stress-induced microcracks, bending-related faults, and deformed serpentine layers may cause anisotropy in the upper continental crust, upper oceanic crust, and mantle wedge (MT), respectively.The tectonic tremors occur at the corner of the MT.

Table 1
Assumed Model for the Receiver Function Synthetic Test , Figures S1 and S2 in Supporting