Temporal Variations in QP−1 ${Q}_{P}^{-1}$ and QS−1 ${Q}_{S}^{-1}$ Above a Megathrust Following Episodic Slow‐Slip Events

Recent observations beneath central Japan have shown that periodic fluid drainage occurs during slow‐slip events (SSEs) based on temporal variations in QP−1 ${Q}_{P}^{-1}$ above the megathrust boundary of the subducting Philippine Sea slab. However, no previous studies have estimated associated QS−1 ${Q}_{S}^{-1}$ . A comparison of QP−1 ${Q}_{P}^{-1}$ and QS−1 ${Q}_{S}^{-1}$ can provide clues to the mechanism of seismic attenuation because of different propagation characteristics of the two waves. We estimate temporal variations in QP−1 ${Q}_{P}^{-1}$ and QS−1 ${Q}_{S}^{-1}$ via spectral analyses of waveform data from November 2009 to August 2021 period. The results indicate that both QP−1 ${Q}_{P}^{-1}$ and QS−1 ${Q}_{S}^{-1}$ exhibit temporal variations at about 1‐year periodicity and there are systematic differences between QP−1 ${Q}_{P}^{-1}$ and QS−1 ${Q}_{S}^{-1}$ that QS−1 ${Q}_{S}^{-1}$ have smaller values, less insignificant variation, and weaker correlations with SSEs. Furthermore, QP−1/QS−1 ${Q}_{P}^{-1}/{Q}_{S}^{-1}$ increases concurrent with SSEs. These differences suggest that attenuation is caused by the wave‐induced fluid flow. QP−1/QS−1 ${Q}_{P}^{-1}/{Q}_{S}^{-1}$ could be an important parameter for detecting the presence of fluid.

and thermoelastic relaxation (e.g., Budiansky et al., 1983;Savage, 1966) have been suggested as potential mechanisms that may be influential in generating −1 ∕ −1 > 1.0. Although these attenuation mechanisms have been discussed by theoretical studies and laboratory rock experiments, very few attempts have been made through seismological observations. We investigate long-term temporal variations in −1 and −1 beneath the southwestern part of Ibaraki Prefecture, where periodic drainage from a megathrust boundary is expected to occur. The purpose of this study is to provide simultaneous observations of −1 and −1 and to discuss the temporal characteristics of −1 and −1 with those of slow-slip events (SSEs) along the megathrust and the activity of supra-slab earthquakes. Based on the observed different characteristics in −1 and −1 values and the theoretical studies on the mechanism of seismic attenuation, we propose a plausible physical process that explains our observations.

Data
The study area and seismicity are shown in Figure 1. We selected the supra-slab earthquakes (colored circles at ≤35 km depth in Figure 1c) from the Japan Meteorological Agency (JMA) hypocenter catalog and repeating earthquakes (black and red stars in Figure 1c) from the repeating earthquake catalog by Igarashi (2020). The repeating earthquakes that met the same criteria outlined in Nakajima and Uchida (2018) are shown by the  (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) during November 2009-August 2021 and the distribution of Metropolitan Seismic Observation Network stations (black triangles) in the study area. 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). Black and red stars represent repeating earthquakes, with red stars denoting the repeaters used to estimate the megathrust 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) Frequencymagnitude ( JMA ) distributions of the supra-slab (red) and repeating (blue) earthquakes. red stars in Figure 1c. We analyzed the waveforms of 272 earthquakes (2.0 ≤ M ≤ 5.5) that occurred during November 2009-August 2021 period, which consisted of 262 megathrust earthquakes (blue circles in Figure  S1 of the Supporting Information S1) and 10 supra-slab earthquakes (red circles in Figure S1 of the Supporting Information S1). Eight seismograph stations used in this analysis are MeSO-net (Metropolitan Seismic Observation network, triangles in Figure 1b) stations, which consist of acceleration seismographs that are installed in a ∼20-m-deep borehole and record at a 200-Hz sampling frequency (Aoi et al., 2021;Sakai & Hirata, 2009). We manually picked P-and S-wave arrival times at the eight stations for 4,312 waveforms. For estimation of megathrust slip rate and the number of supra-slab earthquakes, we set magnitude thresholds of M ≥ 2.5 and M ≥ 1.0 for the repeating and supra-slab earthquakes, respectively, because the frequency-magnitude distributions obey the Gutenberg-Richter law (Gutenberg & Richter, 1944) above these magnitudes (Figure 1d).

Estimation of Corner Frequency, * and −
The details of our methodology are described in Text S1 of the Supporting Information S1. We here explain the outline of our procedures used in the estimation of corner frequency, Δ * and −1 .
We first estimated P-wave corner frequencies ( ) by applying the spectral ratio method to the vertical components of the acceleration waveforms for 202 megathrust earthquakes and 5 supra-slab earthquakes. We calculated each P-wave spectrum using 2-s time window from the onset of the manually picked P-wave arrival time, then we took spectral ratios at common stations for each earthquake pair with a magnitude difference of ≥0.5 and their epicentral distance of ≤ 5.0 km. We did not include spectra with signal-to-noise ratios (S/N) of <2.0 in the 10-32-Hz frequency range. The noise spectra were calculated using a 2-s time window before P waves. When the spectral ratios were obtained at five or more stations, we calculated average spectral ratio and fitted the theoretical spectral ratio that was obtained from the 2 source model (Brune, 1970) to the averaged spectral ratio via grid-search approach in the 1-32-Hz frequency range. The S-wave corner frequencies ( ) were transformed from by assuming the circular crack model, where = 1.3 (Aki & Richards, 2002;Sato & Hirasawa, 1973). We conducted this analysis separately for the megathrust and supra-slab clusters, and determined and for 175 megathrust earthquakes and 4 supra-slab earthquakes ( Figure S2 in Supporting Information S1).
We used a combined set of values that are estimated in this study and by Nakajima and Uchida (2018) to estimate the temporal variations in −1 . We calculated P-and S-wave spectra from the vertical and transverse components, respectively, using 2-s time window from the onset of the manually picked arrival time for each wave. The noise spectra were calculated using 2-s time window before P waves. We then took pairs of earthquakes, one from the megathrust cluster and the other from the supra-slab cluster, and calculated their spectral ratios at common stations to quantify the differential attenuation (∆ * ) between two earthquakes (see Nakajima & Uchida, 2018 for further details). We applied the source spectrum correction of each spectral ratio using the estimated values and then stacked the spectral ratios at every 0.1-year interval throughout the analysis period (a total of 115 windows for December 2009 (2010.0)-May 2021 (2021.4) period) using a centered 0.4-year time window. Stacking was performed using the earthquake pairs with S/N ≥ 3.0 in 20-45-Hz frequency range to ensure stable estimates of attenuation. We then fitted the theoretical spectral ratio to the stacked ratio, which was discretized at every 5 Hz, via the least-squares method for each 0.4-year time window (∆ * was determined from the slope of the straight lines in Figure S3 of the Supporting Information S1 and their temporal variations are shown in Figures S4a and S4b of the Supporting Information S1). Finally, the temporal variations in −1 were estimated as −1 = ∆ * ∕∆ , where ∆ is arrival time difference of the megathrust and supra-slab earthquakes that were averaged over 0.4-year moving time window ( Figure S4c in Supporting Information S1). These procedures were conducted separately for P and S waves to obtain the temporal variations in −1 and −1 .

Estimated Corner Frequencies
The values estimated in this study are shown in Figure S2 of the Supporting Information S1. The values that were determined for 175 megathrust earthquakes and 4 supra-slab earthquakes are inversely proportional to the cube root of seismic moment. We calculated the stress drop (Δ ) from using the Eshelby (1957) formula (see Nakajima & Uchida, 2018, for details). The obtained Δ values are in the 0.1-10 MPa range (average of 3.4 MPa) for supra-slab cluster and the 0.1-100 MPa range (average of 22 MPa) for megathrust cluster. The average stress drops for these two clusters are almost equivalent to the stress drops estimated by Tsuda et al. (2010), Nakajima and Uchida (2018), and Kinoshita and Ohike (2002) for the same region. We will discuss the effect of the uncertainty in on the results in Text S2 of the Supporting Information S1. Figure 2a and Figure S5 in Supporting Information S1 show temporal variations in −1 (circles) and −1 (squares) for December 2009 (2010.0)-May 2021 (2021.4) period. The obtained attenuation values range from almost zero to 4.0-6.3 × 10 −3 for −1 and to 2.0-3.3 × 10 −3 for −1 , and the average values for the entire period are estimated to be 1.9 × 10 −3 for −1 and 1.3 × 10 −3 for −1 (dashed lines in Figure S5 of the Supporting Information S1), which are comparable to estimated values from tomographic inversions beneath Kanto (Nakajima, 2014;Nakamura et al., 2006). The results show −1 have generally smaller values than −1 ( −1 ∕ −1 > 1.0). The power spectral density plots of the −1 time series, which are calculated for the period after 2013 to avoid the effects of the 2011 Tohoku-oki earthquake, indicate that both −1 and −1 possess a 0.9-year periodicity (see red and blue lines in Figure 2b). We assessed the reliability of this periodicity with statistical analysis of 1,000 bootstrap samples (see Text S3 in Supporting Information S1 for details) and confirmed that the 0.9-year periodicity is robust (see black lines in Figure 2b).

Temporal Variations in − and − and Relationships With SSEs
We evaluated the temporal variations in the slip rates along the megathrust, which can serve as a proxy for the occurrence of SSEs, from the individual sequences of repeating earthquakes of Igarashi (2020) Figure 2a and Figure S6a in Supporting Information S1) (e.g., Nakajima & Uchida, 2018;Uchida et al., 2016). We performed cross-correlation analyses between −1 and the megathrust slip rates to investigate temporal relationships between −1 and SSEs. We analyzed the data after 2013 to remove the effect of the 2011 Tohoku-oki earthquake (dashed green line in Figure 2a) on the local seismic activity.
The cross-correlation analyses indicate that −1 and the occurrence of SSEs are positively correlated, with −1 lagging the SSEs response by 0-0.1 years (red and blue circles in Figure 2c) with a statistical significance of P value < 0.05. We evaluated the effect of uncertainty in −1 on the temporal correlations by applying the bootstrap method (see Text S3 in Supporting Information S1 for details) and confirmed that more than 96% of −1 and 74% of −1 have a statistically significant correlation (P value ≤ 0.05) (see dashed black lines in Figures 2c and 2d and Text S3 in Supporting Information S1). Therefore, we consider the correlation between −1 and SSEs to be slightly weaker than −1 , thereby −1 is less affected by the occurrence of SSEs. This finding highlights essential differences between −1 and −1 associated with the periodic drainage from the megathrust. We confirmed that −1 and the supra-slab seismicity (gray shaded regions in Figures S6b and S7a of the Supporting Informa tion S1) are also positively correlated, with the seismicity lagging −1 by 0.0-0.1 years ( Figure S7b in Supporting Information S1), suggesting that the three independent phenomena, SSEs along the megathrust, −1 above the megathrust, and supra-slab seismicity, are highly correlated in time and space. It is noted that each of these three phenomena has the periodicity of ∼1 year ( Figure S8 in Supporting Information S1).  Figure 3a). To characterize the relationship between −1 ∕ −1 and megathrust slip rates, we classified the entire period into low- (Figure 3b) and high-slip-rate periods (Figure 3c), where we defined the low-or high-slip-rate periods as the periods with slip rates that are slower or faster than the average PHS plate subduction rate (5 cm/year), respectively (Heki & Miyazaki, 2001;Seno et al., 1993). The −1 ∕ −1 tends to have higher values during high-slip-rate period than low-slip-rate period (see rose histograms in Figures 3b and 3c). Therefore, we infer that the systematic variations in −1 ∕ −1 occur as a response to the occurrence of SSEs.

Possible Model for − > − and the Temporal Variations in −
Although the major cause of seismic attenuation is considered to be associated with intergranular deformation with −1 ∕ −1 < 1.0 (e.g., Jackson, 2015;Jackson & Anderson, 1970;Jackson et al., 2002), we found that −1 ∕ −1 > 1.0 for most of the analysis period. Therefore, our results suggest that another dominant process is responsible for observed −1 ∕ −1 . Here, we consider two plausible mechanisms, WIFF (e.g., Ba et al., 2008; Power spectral density plots of −1 (red line) and −1 (blue line), respectively. Gray lines denote 1,000 bootstrap samples of the power spectral density obtained by fluctuating −1 within the 2 range of the Gaussian distribution (see Text S3 in Supporting Information S1), and the black line represents their average value. (c) Cross-correlation coefficients between the slip rates and −1 (red circles) and −1 (blue circles). The maximum correlation coefficients are 0.34 at lag = 0.0 and 0.1 year for −1 (red solid circles) and 0.31 at lag = 0.1 year for −1 (blue solid circle). The dashed horizontal lines indicate the cross-correlation coefficients where the null hypothesis that attenuation is randomly enhanced, irrespective of megathrust slip, is rejected at a statistical significance P = 0.05. The black lines are the mean values of correlation coefficients of 1,000 bootstrap samples (gray lines). (d) Histograms of correlation coefficients at the lag with the highest correlation. When the lag time between slip rate and −1 is 0.1 year, 96.1% and 74.7% have statistical significance at P ≤ 0.05 for −1 (red) and −1 (blue), respectively. the observed −1 ∕ −1 > 1.0 trend. Hauksson and Shearer (2006) obtained −1 ∕ −1 ∼ 1.3 in the lower crust of southern California and inferred that WIFF in fluid-saturated rocks is a likely mechanism generating high −1 ∕ −1 . WIFF occurs along the pore-pressure gradient at a mesoscopic scale, which is larger than the average grain size but smaller than the wavelength of propagating seismic waves (see Figure 4a). The fluids that are trapped in the rock mass would flow from low-aspect-ratio pores at higher pore-fluid pressures to high-aspect-ratio pores at lower pore-fluid pressures by compressional stress (0 < t < T/2 in Figure 4a) and versa vice by tensional stress (T/2 < t < T in Figure 4a; flow indicated by blue arrows in Figure 4a) (Müller et al., 2010). The volumetric deformation is cyclically induced by seismic-wave propagation through the rock mass, with the degree of volumetric deformation due to P waves being more dominant compared with that due to S waves. WIFF can therefore explain observed −1 ∕ −1 > 1.0 trends (Quintal et al., 2012;Zheng et al., 2017). Stachnik et al. (2004) identified areas with −1 ∕ −1 > 1.0 along the Moho in the Alaska subduction zone and inferred that thermoelastic relaxation is a dominant process generating these high −1 ∕ −1 (e.g., Budiansky et al., 1983;Savage, 1966). Thermoelastic relaxation occurs as a rock mass approaches thermal equilibrium along the temperature gradient owing to either compression or tension during seismic-wave propagation, with this process being dominant in rocks that possess mesoscopic heterogeneities in their bulk modulus and thermal expansion, and can produce substantially enhanced −1 (e.g., Jackson & Anderson, 1970;Savage, 1966).
The observed −1 ∕ −1 > 1.0 trends in this study can be explained by either WIFF or thermoelastic relaxation. However, given the interpretation of Nakajima and Uchida (2018) that SSEs induce drainage from the megathrust, we infer that WIFF is a more likely mechanism to explain the observed −1 ∕ −1 > 1.0 trends and temporal variations in −1 . Once a large amount of fluids has drained from the megathrust boundary owing to SSEs (processes [1] and [2] in Figure 4b), cracks are forced to open owing to the high pore-fluid pressures (process [2] in Figure 4b). A coupled process of fluid inflow with WIFF enhances −1 relative to −1 during seismic-wave propagation through the rock mass (Figure 4a). The drained fluids continue to migrate upward over time and trigger supra-slab seismicity (process [3] in Figure 4b); some cracks then close, causing subsequent reduction in −1 . This process is then repeated at a 0.9-1.1-year periodicity (Figure 4c).

Crack Size Estimation
The crack size that accommodates WIFF can be estimated with the model of Hudson et al. (1996Hudson et al. ( , 2001 by assuming the rock and fluid properties. If fluid-filled cracks at a mesoscopic scale are distributed in a medium, as shown in Figures 4a and 4b, then their characteristic radius ( ) can be expressed as a function of the char acteristic frequency ( ), permeability ( ), shear modulus ( ), porosity of the rock matrix ( ), fluid viscosity ( ),  (Dziewonski & Anderson, 1981), = 3.5 × 10 −4 Pa·s, = 8.5 × 10 9 Pa, = 0.01. and were taken from the fluid thermo-physical database provided on the National Institute of Standards and Technology website (Lemmon et al., 2018), assuming the temperature of 450-500 K (Wada & He, 2017) and confining pressure of 1 GPa. We adopted based on the velocity-porosity relation (Marquis & Hyndman, 1992), assuming crack aspect ratio of 0.01 and P-wave velocity of 6.5 km/s (Hasegawa et al., 2007). We estimated characteristic crack size (2 ) to be 10-20 cm in the 20-45-Hz frequency range. The estimated crack size is comparable to the size of fractures that are generated due to fluid penetration from the megathrust in pressure conditions comparable to the 30-40 km depths (Padrón-Navarta et al., 2010).

Conclusions
We investigated the long-term temporal variations in −1 , −1 , and −1 ∕ −1 beneath the Kanto district, Japan, during November 2009-August 2021. We discussed the temporal correlation of −1 with both the timing of SSEs Figure 4. (a) Schematic mechanism of wave-induced fluid flow due to P-wave propagation. When a volumetric stress acts on a medium containing strong mesoscopic-scale heterogeneity, WIFF-induced relaxation from the pores with low aspect ratios to those with high aspect ratios occurs, and the energy loss ( −1 ) occurs by internal friction due to the fluid relaxation (modified from Müller et al. (2010)). T is the characteristic period of P wave. (b) Schematic interpretation of the fluid drainage process (modified from Nakajima and Uchida (2018)). (c) Periodic variations (0.9-1.1 years) in the slip rate (green) along the megathrust, attenuation (red and blue) above the megathrust, and supra-slab seismicity (black). Temporal evolutions of each phenomena, which are interpreted to occur due to periodic drainage from the megathrust, are indicated by blue arrows.
along the megathrust and the seismic activity in the overlying plate and proposed a plausible attenuation mechanism. The major results in this study are summarized as follows.
1. We estimated the temporal variations in −1 and −1 at 20-45 Hz using 4,312 seismic waveforms that were recorded during the November 2009-August 2021 period. 2. The results revealed that both −1 and −1 exhibit marked temporal variations, with a periodicity of approximately 1 year but −1 variation is less significant compared to −1 . 3. We observed spatiotemporal correlations between the occurrence of SSEs along the megathrust, −1 above the megathrust, and the supra-slab seismicity. −1 ∕ −1 is enhanced accompanied by the occurrence of SSEs. 4. We attributed the observed −1 ∕ −1 > 1.0 trends to WIFF and estimated the characteristic crack size to be 10-20 cm.
Our results will lead to a better understanding of seismic-wave attenuation mechanisms and fluid-water interactions around megathrust boundaries. Further studies should include seismic-velocity and shear-wave-splitting analyses to provide additional constraints on the fluid-rock interactions above the megathrust during SSEs.