Temporal Variations in Frequency‐Dependent Shear‐Wave Anisotropy Above a Plate Interface Following Episodic Slow‐Slip Events

Recent observations beneath Kanto, Japan have shown that seismic activity and seismic attenuation within the overlying continental plate change with time due to drainage caused by slow‐slip events (SSEs) along the upper boundary of the Philippine Sea plate. However, associated changes in rock properties have not been investigated. In this study, we estimate frequency‐dependent shear‐wave anisotropy to provide a detailed insight into the structural change associated with drainage. We perform shear‐wave splitting analysis in frequency ranges of 1–4, 2–6, and 4–8 Hz for 306 earthquakes that occur during September 2009–August 2021 and recorded at the Metropolitan Seismic Observation network. Obtained time differences between fast and slow S waves (delay time) range from almost zero to 0.16–0.18 s, exhibiting spatio‐temporal variation and frequency dependence. The fast S‐wave polarization directions are generally consistent with the direction of the maximum horizontal compressional axis in the study region, which suggests that the observed anisotropy is probably caused by the NE–SW‐oriented fractures developed under the regional stress field. The temporal variation in delay times is correlated with SSEs activity with a lag time of 0.0–0.1 year. Furthermore, comparisons between observed frequency‐dependent delay times and numerical calculation of fracture‐induced anisotropy suggest that the average fracture radius is almost constant (0.30–0.35 m) over time but fracture density temporally varies from 0.025 to 0.035. We infer that the fracture density is probably enhanced by opening of the NE–SW‐oriented fractures during the upward migration of fluids that are expelled from the plate interface.


Introduction
Recent studies have shown that fluid supply from subducting plates into the overlying continental plate (upper plate) could trigger earthquakes in the crust and mantle wedge (e.g., Halpaap et al., 2019;Nakajima & Hasegawa, 2016).In addition, structural changes in the upper plate have also been observed in some subduction zones.For example, in the northern part of the Hikurangi subduction zone, seismic activity in the upper plate increased after the occurrence of shallow slow-slip events (SSEs) (e.g., Shaddox & Schwartz, 2019), which was concurrent with seismic velocity and anisotropy changes in the upper plate (e.g., Wang et al., 2022;Zal et al., 2020).One of the possible causes to explain the coupled process between enhanced seismicity and structural changes in the upper plate is drainage due to SSEs by breaking impermeable seals along the plate interface (e.g., Wang et al., 2022).Nakajima and Uchida (2018) reported long-term spatiotemporal correlations between seismicity in the upper plate (hereafter called supra-slab earthquakes) and repeating earthquake activity along the plate interface of the subducting Philippine Sea (PHS) plate beneath Kanto region, Japan (Figure 1).Associated temporal changes in Pand S-wave attenuation in the upper plate have been investigated, and it is suggested that drainage from the plate interface repeatedly occurs in response to SSEs in approximately 1 year intervals (e.g., Ito & Nakajima, 2023;Nakajima & Uchida, 2018).However, it remains unclear how the rock property (e.g., size and shape of fractures and cracks) changes within the upper plate when fluids are supplied.
Seismic anisotropy in a medium with fluid-filled fractures is considered to be caused by the rapid movement of fluids between fractures and pores in the rock matrix (Thomsen, 1995).The frequency dependence of anisotropy generated by this mechanism was theoretically investigated by Hudson et al. (1996) and Pointer et al. (2000).Since their model is composed of micro cracks on a grain-size scale, the theoretical model can be applied to anisotropy of small rock samples measured for sonic and ultrasonic laboratory measurements.However, the model was not comparable to anisotropy of rocks that contain fractures larger than grain sizes as observed in a frequency range of seismic wave.Chapman (2003) established a theory that describes the frequency dependence of seismic anisotropy by incorporating meso-scale fractures, which are much larger than the grain size.This theory was subsequently applied to the shear-wave splitting observed in the vertical seismic profiling data by Maultzsch et al. (2003), and their calculations produced anisotropy compatible with the observations of anisotropy induced by rocks containing meso-scale fractures.
In this study, we investigate temporal variations in frequency-dependent S-wave polarization anisotropy in the upper plate beneath Kanto, Japan (Figure 1) and discuss its correlation with SSE activity along the plate interface.Seismic anisotropy provides clues to infer the density and orientation of fractures, and its frequency dependence can potentially constrain the size of the predominant fractures (e.g., Chapman, 2003;Maultzsch et al., 2003).Therefore, if seismic anisotropy and its frequency dependence change associated with inferred drainage during SSEs, then we can evaluate the temporal variations in the orientation, intensity, and size of fractures above the plate boundary.

Data
We carried out shear-wave splitting analysis using waveforms of 306 earthquakes that occurred during September 2009-August 2021, whose P-and S-wave arrival times were manually picked in this study.Hypocenters of these earthquakes were determined by the Japan Meteorological Agency (JMA) using a one-dimensional velocity model (JMA2001) (Ueno et al., 2002).The analyzed earthquakes consist of 262 plate-interface earthquakes (2.5 ≤ M ≤ 5.5) at depths of 40-50 km and 44 supra-slab earthquakes (1.5 ≤ M ≤ 3.0) at depths of 30-35 km (Figure 1c).Waveforms used for the splitting analysis are the same as those used for the analysis of seismic attenuation by Ito and Nakajima (2023), which enables direct comparisons of the results in this study with the results by Ito and Nakajima (2023).Three seismograph stations (E.YSPM, E.NKGM, and E.NSUM) used in this analysis are a part of the Metropolitan Seismic Observation network (MeSO-net) stations (triangles in Figures 1b  and 1c).The stations consist of three-component acceleration seismographs that are installed in ∼20 m-deep boreholes and are recorded at a 200 Hz sampling frequency (Aoi et al., 2021;Sakai & Hirata, 2009).Seismic waveforms were filtered in three frequency ranges of 1-4, 2-6, and 4-8 Hz because high signal-to-noise ratios in these frequencies ensure stable analyses.Incident angels of seismic waves to the stations are all within 15°so that the contamination by converted waves can be negligible (Figure S1 in Supporting Information S1).
To estimate the slip rate along the upper surface of the PHS plate, which is used to compare with the shear-wave splitting result, we selected repeating earthquakes beneath the supra-slab earthquakes that occurred during the 2009-2019 period from the repeating earthquake catalog by Igarashi (2020) (stars in Figure 1c).The criteria to select repeating earthquakes were the same as in Nakajima and Uchida (2018) who used repeating earthquakes that occurred within a ±3 standard deviation distance from the epicentral centroid of the supra-slab earthquakes (gray stars in Figure 1c and gray bins in Figure 1d).We first estimated the seismic moment (M 0 in dyne•cm) of the selected repeating earthquakes from the JMA magnitude (M JMA ) with logM 0 = 1.5MJMA + 16.1.We then calculated slip (d in cm) for each repeating earthquake with the relationship of d = 10 2.36 × M 0.17 0 (Nadeau & Johnson, 1998), which is used as a proxy for representing the slip along the plate interface.The slip estimated from each repeating earthquake is shown in Table S1.

Estimation of Delay Time and Fast Direction
In an isotropic medium, an S wave that is radiated from the source arrives at a station preserving the original polarization direction.However, in an anisotropic medium, the S-wave velocity depends both on the polarization direction of the wave and on inclination and azimuth of the ray path to stations (e.g., Nur & Simmons, 1969).Thus, an S wave vibrating or traveling in a certain direction splits into two waves with mutually orthogonal polarizations except for the case in which the original S-wave polarization is parallel or perpendicular to the symmetry plane of the anisotropic medium.The split waves arrive at a station with a lag time.This phenomenon is known as "shear-wave splitting" (e.g., Savage, 1999).We can characterize the strength and direction of anisotropy by measuring two parameters, the time difference between fast and slow S waves (delay time) and the leading S-wave polarization direction (fast direction).We used the cross-correlation method (e.g., Nakajima & Hasegawa, 2004;Silver & Chan, 1991) to estimate the delay time and fast direction of shear-wave splitting.We can calculate the splitting parameters by rotating two horizontal (NS and EW) components of the observed seismic waveform to the fast and slow S-wave components (red and blue lines in Figure 2a).Since the fast and slow S waves originate from a single S wave that was radiated from the source, the two waveforms after the rotation are expected to be similar with each other when the lag time is corrected (black dashed line in Figure 2a).We rotated NS and EW components of the observed seismograms from 0°to 175°with an interval of 5°, and one of the rotated components is shifted by a time lag, varying from 0 to 1 s with an interval of 0.005 s.When the value of the cross-correlation coefficient becomes maximum, the rotation angle was regarded as the fast direction, and the amount of the time lag was regarded as the delay time.A horizontal particle motion of the original waveforms exhibits an ellipse shape for an anisotropic media, while that corrected for the anisotropic effect becomes nearly linear (Figure 2b).
We applied the multiple window analysis (Savage et al., 2010;Teanby et al., 2004) to avoid the bias of the length and the onset of time windows on the splitting parameters.We set 250 time windows for each waveform that start from 0.1 to 0.0 s from the manually picked S arrival with a 0.02 s interval and end from +0.4 to +0.65 s from the S arrival with a 0.005 s interval (see gray shaded region in Figure 2a).We visually ensured that these time windows include one or two cycles of the S waves and excluded waveforms that do not contain the first one cycle of the S wave within a given time window.We then classified the obtained 250-splitting parameter pairs into several dominant clusters by the k-means method (MacQueen, 1967) based on the scattered distribution of the delay times and fast directions.When a given cluster had the largest number of parameter pairs (N ≥ 100), we considered that the cluster as the optimal sets of the delay time and fast direction, and the parameter pair with the minimum estimation error within that cluster was used in the subsequent analysis (Teanby et al., 2004).Figure 2c shows examples of the estimates of splitting parameters for two earthquakes.We obtained stable measurements of splitting parameters for the 250 windows for an earthquake (2009-11-24), and the splitting parameter pairs were classified into a single cluster (left panel in Figure 2d).We regarded the delay time and fast direction at a time window at which we observed the minimum estimation error as the optimal delay time and fast direction (red cross).In contrast, splitting parameters estimated for an earthquake (2012-12-16) were somewhat scattered (right panel in Figure 2c), and the parameter pairs were classified into three clusters with the number of observations of 151, 77, and 22 (right panel in Figure 2d).We regarded a pair of splitting parameters with the minimum estimation error in the cluster with the maximum number of parameter pairs (N = 151 in this case) as the optimal pair of delay time and fast direction (red cross in Figure 2d).
We performed this analysis for 262 plate-boundary and 44 supra-slab earthquakes and obtained 225, 226, and 187 splitting parameters for plate-interface earthquakes and 17, 28, and 12 parameters for supra-slab earthquakes, respectively, at 1-4 Hz, 2-6 Hz, and 4-8 Hz.We confirmed that there are no systematic biases in the estimates of the delay times and fast directions with respect to the number of clusters classified by the k-means method (Figure S2 in Supporting Information S1).We visually checked each waveform and excluded the measurement when cycle skipping occurred.In the result section, we show the results of E.YSPM station (red triangle in Figure 1), where we observed the highest signal-to-noise ratio among the three stations.The analysis results at other two stations (E.NKGM and E.NSUM; black triangle in Figure 1) are shown in Figures S3 and S4 in Supporting Information S1.

Estimated Delay Times and Fast Directions
The obtained delay times for plate-interface earthquakes range from almost zero to 0.16-0.18s for all frequency bands (Figures 3b-3d).We found that the delay times gradually decrease for all frequency bands at the horizontal distance of ≥∼14 km (the horizontal distance is calculated from point B in Figure 1b), even though the lengths of ray paths to stations become longer at larger horizontal distances (Figure 3a).The delay times of the supra-slab earthquakes range from 0.08 to 0.12 s, and these values are generally smaller than most of the delay times for the plate-interface earthquakes.Therefore, we infer that an anisotropic medium exists between the plate interface and the supra-slab earthquakes at horizontal distances of 8-14 km.It is noted that the delay times observed for the supra-slab earthquakes suggest that a medium above the supra-slab earthquakes also possesses some degree of anisotropy.The splitting parameters estimated in this study are listed in Tables S2-S7.
Figure 4 shows rose diagrams of the fast directions inferred from the plate-interface and supra-slab earthquakes.Both the plate-interface and supra-slab earthquakes show the fast directions in N30°-70°E in all frequency bands, and their average directions are 53°(1-4 Hz), 47°(2-6 Hz), and 46°(4-8 Hz) for plate-interface earthquakes, and 64 (1-4 Hz), 48°(2-6 Hz), and 56°(4-8 Hz) for supra-slab earthquakes.These directions are largely comparable to the direction of the maximum horizontal compressive axis (∼60°) inferred from active folds in the Kanto region (purple lines in Figure 4) (Sugiyama et al., 1997).This suggests that anisotropy observed in this study is probably caused by the preferable alignments of fractures or cracks that have been developed under the NE-SW-oriented regional stress field (e.g., Crampin, 1994;Crampin & Peacock, 2008).

Temporal Variations in Delay Times and Relationships With SSEs
To investigate temporal relationships between SSEs occurring along the plate interface and delay times, we plot the individual delay times of plate-interface earthquakes for July 2009-June 2021 (colored circles in Figure 5).We here adopted the delay times of plate-interface earthquakes at the horizontal distance of 8-14 km to extract earthquakes whose ray paths are considered to be affected by the strong anisotropic medium above the plate interface.The delay times of ≥0.08 s were selected so that the delay times for plate-interface earthquakes are higher than or comparable to those for supra-slab earthquakes within 95% confidence intervals (see Figures 3b-3d).The total numbers of delay times shown in Figure 5 are 91, 90, and 76 at 1-4, 2-6, and 4-8 Hz, respectively.
We evaluated the temporal variations in the plate-interface slip rates, which can serve as a proxy for the occurrence of SSEs, from the individual sequences of repeating earthquakes by calculating 0.4 year centered average slip rates at every 0.1 year interval (gray region in Figure 5) (e.g., Nakajima & Uchida, 2018; Uchida Journal of Geophysical Research: Solid Earth 10. 1029/2023JB028526 et al., 2016)).Figure 6 shows temporal variations in the delay times for plate-interface earthquakes and the plateinterface slip rates that are both averaged in 0.4 year time windows at every 0.1 year interval.
Figure 7 shows the results of the cross-correlation analysis between the plate-interface slip rates and delay times within 0.4 year sliding windows for the 2009-2019 period.We find that the delay times have the strongest temporal correlation with the plate-interface slip rates, which is used as a proxy for SSEs, with a statistical significance of P value < 0.05 when the lag time of delay times against the slip rates is 0.0-0.1 years (color-filled circles in Figure 7a).Of note, similar correlation between the two independent phenomena are coherently observed at the three different frequency bands and at other two stations (Figure S4 in Supporting Information S1).Therefore, we consider that an anisotropic structure above the plate interface changes with time in response to the occurrence of SSEs.The lag times of 0.0-0.1 years between the occurrence of SSEs and delay times are almost comparable to those observed for the temporal variations in the strength of P-and S-wave attenuation above the plate interface (Ito & Nakajima, 2023;Nakajima & Uchida, 2018).
In the cross-correlation analysis, we used only the delay times of ≥0.08 s for the plate-interface earthquakes that occurred in the horizontal distance of 8-14 km.To assess the effects of these threshold values on the temporal relationship between the slip rates and delay times, we performed the same analysis for two delay-time values (0.00 and 0.04 s) and three values of horizontal distance ranges for plate-interface earthquakes (8-13 , 8-15 , and 8-16 km).The results show that temporal correlations calculated by the different threshold values are almost identical to those observed for the plate-interface earthquakes in the horizontal distance of 8-14 km with the delay times of ≥0.08 s (Figures S5 and S6 in Supporting Information S1).We also confirmed that different lengths of time windows for the calculation of the average slip rates and delay times yield similar cross-correlation patterns.
The time plots of delay times derived from supra-slab earthquakes are shown in Figure S7 in Supporting Information S1.In contrast to plate-interface earthquakes, the delay times for supra-slab earthquakes are limited to a small number of data, and we cannot reliably assess a temporal variation in anisotropy above the supra-slab earthquakes.However, the delay times for supra-slab earthquakes are generally smaller than those of plateinterface earthquakes over the entire analysis period.Therefore, we infer that an anisotropic region between the plate interface and supra-slab earthquakes shows a temporal variation concurrent with SSEs along the plate interface.

Fracture Property Estimation From Frequency-Dependent Anisotropy
The strength of anisotropy was calculated using the formula Φ = dt L × V S × 100 (Hudson, 1981) where dt and L are the individual observed delay time (≥0.08 s) and ray-path length to E.YSPM station from plate-interface  (Sugiyama et al., 1997).

Journal of Geophysical Research: Solid Earth
10.1029/2023JB028526 earthquakes, respectively.We assumed an average isotropic shear-wave velocity V S of 3,636 m/s at depths of 0-40 km, which was calculated from the JMA 2001 velocity model (Ueno et al., 2002).Figure 8a shows temporal variations in the strength of anisotropy at 1-4, 2-6, and 4-8 Hz, and Figure 8b shows the average the strength of anisotropy at each frequency over the entire period.
As shown in Figure 8b, it is found that the average anisotropy strength becomes weaker in higher frequency bands.Such frequency-dependent anisotropy has been observed in a highly fractured geologic layer in eastern Utah, the United States, and is discussed in terms of the effective size and density of fluid-filled fractures (e.g., Chapman, 2003;Maultzsch et al., 2003).It is noted that frequency-dependent anisotropy could be apparently caused by the difficulty to determine delay times that are too small or large compared to a period of waveform (e.g., Silver & Chan, 1991;Wolfe & Silver, 1998).However, periods of waveforms used in the splitting analysis (0.125 ∼ 1 s) are on the same order of magnitude as the observed delay times (0.08-0.18 s).Therefore, we infer that the frequency-dependent anisotropy observed in this study is not apparent and is caused by drainage-related physical process above the subducting PHS plate.We used the theoretical model of frequency-dependent anisotropy derived by Chapman (2003) and evaluated the temporal change in the average radius and density of fractures using the temporal variations in the strength of the observed anisotropy.We calculated the theoretical strength of anisotropy as a function of frequency of seismic wave using equations (S1)-(S7) in Text S1 in Supporting Information S1.Then, we estimated the average fracture radius and fracture density by fitting the theoretically calculated strength of anisotropy to the observed one.In the theory of Chapman (2003), two types of inclusions with different characteristic dimensions are defined: one is "micro crack" on a grain-size scale and the other is "fracture" on a scale much larger than grain size.Since micro crack is generally considered to have no influence on the frequency dependence of anisotropy for the frequency band of seismic waves, the micro-crack density was assumed to be 0 in our calculation, and only the fracture density and fracture radius were treated as free parameters.We fixed other parameters except for the fracture radius and density to typical values proposed in previous studies (see Table S8 in Supporting Information S1).We classified the observed strength of anisotropy into "strong-anisotropy period" and "weak-anisotropy period" when the strength of anisotropy is larger and smaller than 95% confidence intervals of the mean value for the entire period, respectively (Figure 8a).We then calculated the mean values of anisotropy strength for the strong and weak anisotropy periods, respectively, for the three frequency bands (Figures 8c and 8d), and compared the mean values with the theoretically calculated strength of anisotropy to estimate the fracture density and fracture radius at each period.
Figures 9a and 9b show the results of fitting between the observed and theoretical strength of anisotropy at the strong and weak anisotropy periods, respectively.The fracture radius and density were estimated from ranges of the theoretical anisotropy strength that can cover the observed anisotropy strength at 1-4, 2-6, and 4-8 Hz (gray bands in Figures 9a and 9b).The obtained results show that the fracture density ranges from 0.019 to 0.048 with the mean of 0.035 for the strong anisotropy period and from 0.014 to 0.035 with the mean of 0.025 for the weak anisotropy period (Figure 9c).In contrast, the fracture radius ranges from 0.05 to 2.34 m with the mean of 0.35 for the strong anisotropy period and from 0.05 to 1.51 with the mean of 0.30 for the weak anisotropy period.These results suggest that the average fracture density increases with increasing anisotropy strength but the average fracture radius changes little with time.

Possible Mechanism for Seismic Anisotropy Change
Our observations show that the temporal variation in delay times is correlated with SSEs activity with a lag time of 0.0-0.1 year (Figure 7).Seismic structure changes within the overlying continental plate associated with SSEs have been suggested from the analyses of P-and S-wave attenuation structure in the same area (Ito & Figure 7. (a) Cross-correlation coefficients at E.YSPM station between the plate-interface slip rates and the delay times averaged in 0.4 year time windows, which were calculated for the 2009-2019 period.Colored-filled circles indicate the cross-correlation coefficients where the null hypothesis that the delay time is randomly distributed, irrespective of plate-interface slip, can be rejected with a statistical significance P ≤ 0.05.The maximum cross-correlation and lag time values are shown.(b) Scatter plots of the slip rate and delay time that are shifted by the lag time that yields the maximum cross-correlation coefficients.Nakajima, 2023;Nakajima & Uchida, 2018).Similar structural changes concurrent with SSEs have been identified through the investigation of seismic velocity and anisotropic structures in the northern part of the Hikurangi subduction zone (e.g., Wang et al., 2022;Zal et al., 2020).These structural changes above the plate interface concurrent with SSEs are interpreted as being caused by drainage from the plate interface.This study further suggests that the average fracture density increases when the anisotropy strength is enhanced (Figure 9).
Here we proposed a simple process that can explain the temporal changes in seismic anisotropy observed in this study.First, a low-permeability seal along the overpressurized plate boundary is broken by periodic occurrences of SSEs (Ito & Nakajima, 2023;Nakajima & Uchida, 2018;Wang et al., 2022;Zal et al., 2020).Then, fluids liberated from the plate boundary permeate gradually into the upper plate, which probably forces pre-existing, preferably oriented fractures to open, and the strength of anisotropy is consequently enhanced (Figure 10a).Finally, when drainage ceases after the SSEs, the plate boundary is sealed again by cementation in the active hydrothermal environment and the fluid flux toward the upper plate is reduced (e.g., Saishu et al., 2017;Sibson, 2013;Williams & Fagereng, 2022).At this stage, some of the aligned fractures would close and anisotropy becomes weakened in inter-SSE periods (Figure 10b).This process would yield that the fracture density above the plate interface increases immediately after SSE.However, data points of the strong and weak anisotropy periods defined in this study are not completely periodic (Figure 8), and whether the average fracture density always increases after SSEs and decreases in inter-SSE periods is not well quantified.Further investigation is required to quantify the temporal correlation between the increase in the fracture density and the occurrence of SSEs.If we assume that the opening and closure of the aligned fractures occur in response to drainage from the plate interface, it is likely that permeability, which is proportional to the fracture density, changes with time (e.g., Benson et al., 2006;Gueguen & Dienes, 1989;Guéguen & Schubnel, 2003).Nakajima and Uchida (2018) estimated permeability of as high as 10 14 -10 13 m 2 after SSEs, and such high permeability may be caused by the opening process of the fractures due to the inflation of overpressurized fluids into the anisotropic medium.The high permeability can be supported by high fluid transport velocities that are estimated from fluid signatures within fractured rock in field observations (Yoshida et al., 2023).We infer that the reason why the fracture radius does not change during this process is that the radius reflects the size of fractures inherently distributed in rocks and the size of fractures is probably constant over time.

Conclusions
We investigated the long-term temporal variations in S-waves polarization anisotropy beneath Kanto district, Japan, using high-quality waveform data for the period of September 2009-August 2021.We discussed the temporal correlation between the observed delay times between fast and slow S waves and the timing of SSEs inferred from repeating earthquake activity along the plate interface, and compared the observed frequencydependent anisotropy with theoretically predicted ones.The major results obtained in this study are summarized as follows.
1. We estimated delay times at 1-4 Hz, 2-6 Hz, and 4-8 Hz using 306 seismic waveforms and revealed that the delay times exhibit marked temporal variations in all frequency ranges.2. We observed temporal correlations between the delay times and the occurrence of SSEs along the plate interface, where the delay times are enhanced after the SSEs with a lag of 0.0-0.1 years.3. Comparisons of our results with numerical calculations of fracture-induced anisotropy suggested that the average fracture density above the plate interface temporally changes from ∼0.025 to ∼0.035 but the average fracture radius little change with time.The temporal change in the fracture density may be caused by the opening and closure of fractures as a response to drainage from the underlying subducting plate interface.
Our results support the conclusion in the previous studies that seismic structure changes in response to drainage from the plate boundary and provide a more detailed and additional insight into the structural changes in the target area.The fracture properties estimated in this study are important information to understand hydrological process ongoing in the rocks above the plate interface.Future studies could involve a quantitative estimate of fluid migration process into the overlying plate through a combined analysis of observational results, theoretical experiments, and field studies.

Figure 1 .
Figure 1.(a) Tectonic setting in central Japan.Purple lines denote the iso-depth contours (in km) of the Philippine Sea (PHS) plate (Hirose et al., 2008; Nakajima et al., 2009).Black lines with the teeth denote the troughs.The red rectangle shows the study area.(b) Seismicity (colored circles; M ≥ 0.0) during September 2009-August 2021 and the location of three Metropolitan Seismic Observation Network (MeSO-net) stations (triangles).A red triangle denotes station E.YSPM.The seismicity is color-coded according to the focal depth.(c) Cross-sectional view of the earthquake distribution along the A-B profile in (b).Horizontal axis denotes the distance from point B along the A-B profile.Red and blue circles represent supraslab and plate-interface earthquakes, respectively, for which splitting parameters were estimated at 1-4 Hz at station E.YSPM.Dark and light gray stars represent repeating earthquakes, with dark stars denoting the repeaters used to estimate the plate-interface slip rates.The dashed line denotes the Moho depth (Matsubara et al., 2017) and the purple line denotes the upper surface of the PHS plate.(d) Frequency distribution of repeating earthquakes along the A-B profile in (b).Dark gray bins denote the number of the repeaters used for the estimation of the plate-interface slip rates.

Figure 2 .
Figure 2. Examples of shear-wave splitting analysis at E.YSPM station.(a) Solid black lines show the NS and EW components of the observed waveform.Red and blue lines show the fast and slow S waves after rotation, respectively, and dashed line shows the slow S wave corrected for the lag time.Gray shaded regions represent the ranges of start and end times used for setting 250 time windows, and horizontal bar denotes the length of the time window used in the splitting calculation in (a).(b) (left) the particle motion of the original horizontal component, and (right) the particle motion after correcting for the anisotropic effect.The cross-correlation coefficients, the estimated fast direction and delay time are shown at the bottom.(c) Two earthquake examples of delay times and fast directions estimated for the 250 different time windows.The window numbers are assigned in order of the total length of the time window.Vertical bars at each point represent the 95% confidence region obtained from the t-test.(d) Distribution of delay times and fast directions for the 250 time windows (black circles) analyzed for the two earthquakes in (c).Vertical and horizontal bars represent 95% confidence region of the fast direction and delay time at individual points, respectively.Red cross indicates the location of measurement values with minimum estimation errors (95% confidence regions) with the largest number of parameter pairs, which is regarded as the optimal pair of the delay time and fast direction (see Section 2.2).

Figure 3 .
Figure 3. (a) Cross-sectional view of the earthquake distribution used for shear-wave splitting analysis.Black circles denote the plate-interface earthquakes used to estimate temporal variation in delay times (see Section 3.2).Horizontal axis denotes the distance from point B along the A-B profile in Figure 1b.(b) The delay times at E.YSPM station with respect to the horizontal locations of earthquakes at 1-4 Hz.Colored circles denote the delay times for plate-interface earthquakes used for the estimation of the temporal variation, while white circles denote the delay times for supra-slab earthquakes.Vertical error bars denote the 95% confidential regions of the best measurements (see Method).(c) The same as in (b) but at 2-6 Hz.(d) The same as in (b) but at 4-8 Hz.

Figure 4 .
Figure 4. Rose diagrams of the fast direction inferred from (a) plate-interface and (b) supra-slab earthquakes at E.YSPM station for the frequency bands at (left) 1-4 Hz (middle) 2-6 Hz, and (right) 4-8 Hz.The average fast direction is shown at each panel.Purple lines denote the direction of the maximum horizontal compressional axis inferred from active folds in the Kanto region(Sugiyama et al., 1997).

Figure 5 .
Figure5.Temporal variations in individual delay times inferred from plate-interface earthquakes at E.YSPM station for the frequency bands of (top) 1-4 Hz (middle) 2-6 Hz, and (bottom) 4-8 Hz, and in the slip rate along the plate interface (shaded area).The slip rates were averaged in 0.4 year time windows at every year interval.Vertical bars denote 95% confidential intervals of the estimated delay times.Vertical lines denote the timing of the 2011 Tohoku-oki earthquake.

Figure 6 .
Figure 6.Temporal variations in the delay times inferred from plate-interface earthquakes at E.YSPM station that are averaged in 0.4 year time windows at every 0.1 year interval for the frequency bands of (top) 1-4 Hz (middle) 2-6 Hz, and (bottom) 4-8 Hz.Vertical bars denote the uncertainty of the average delay time calculated by the error propagation from individual measurements of delay time, and horizontal bars denote the time windows with 0.4 year length.Other symbols are the same as in Figure 5.

Figure 8 .
Figure 8.(a) Temporal variations in the strength of anisotropy inferred from plate-interface earthquakes for the frequency bands of (top) 1-4 Hz (middle) 2-6 Hz, and (bottom) 4-8 Hz.Colored-filled and color-bordered circles indicate the strength of anisotropy that are larger and smaller than 95% confidence intervals of the mean value (horizontal-colored bands), respectively.(b) Mean values for the anisotropy strength for the three frequency bands for the entire period.Vertical bars indicate the 95% confidence regions.(c) Mean values of the anisotropy strength for a strong-anisotropy period.(d) Mean values of the anisotropy strength for a weak-anisotropy period (see text for details).

Figure 9 .
Figure 9. Examples of comparisons between observed and theoretical anisotropy strength for (a) strong and (b) weak anisotropy periods.Colored lines denote the observed strength of anisotropy at each frequency band.Dashed lines show examples of theoretical strength of anisotropy for various fracture radii with the fixed fracture densities of 0.035 and 0.025 at strong and weak anisotropy periods, respectively.Gray bands are acceptable ranges of the fracture radii and densities that explain the observed strength of anisotropy.(c) Distribution of fracture radii and densities that can explain the observed anisotropy strength in (red) strong and (blue) weak anisotropy periods.Triangles indicate the centroids of the plots at each period, and red and blue bands denote the 95% confidence regions in the mean values of fracture radii and densities for strong and weak anisotropy period, respectively.

Figure 10 .
Figure 10.Schematic interpretation of drainage process in response to SSEs (modified from Nakajima & Uchida, 2018) and associated changes in seismic anisotropy.(a) Opening process of fractures during SSE period and (b) closing process of fractures at an inter-SSE period.