Slow Earthquakes Illuminating Interplate Coupling Heterogeneities in Subduction Zones

Slow earthquakes are mainly distributed in the vicinity of seismogenic zones of megathrust earthquakes and relationships between both types of earthquakes are expected. We examined the activity of very low frequency earthquakes (VLFEs), classified as one type of slow earthquake, around Japan because they have the potential to clarify detailed spatiotemporal slip behaviors at the plate boundaries. The distribution of the shallow VLFE activity rate is heterogeneous along trench axes and exhibits an anticorrelation relationship with the spatial distribution of the interplate coupling ratio, whereas deep VLFEs are distributed only in weakly coupled areas, and the spatial variation of the activity rate is small. Furthermore, VLFEs are mainly hosted by low seismic velocity anomalies. Thus, slow earthquakes can be triggered by decreased effective stress due to the high pore fluid pressure within regions with weak interplate coupling, and their activity can be an indicator of interplate slip behavior.


Introduction
Slow earthquakes mainly occur between seismogenic and stable sliding zones along the plate boundaries of subduction zones (Obara & Kato, 2016) and are considered to be transitional phenomena between them. The spatial variation of the slip properties at the plate boundary must be controlled by heterogeneous frictional conditions (Obara & Kato, 2016). Various slow earthquakes, such as low frequency tremors (2-8 Hz), very low frequency earthquakes (VLFEs; 0.02-0.05 Hz), slow slip events (SSEs), and coupled phenomena (episodic tremor and slip, ETS) have been detected in many subduction zones worldwide (e.g., Hutchison & Ghosh, 2019;Ito et al., 2007;Obara, 2002;Obara & Ito, 2005;Rogers & Dragert, 2003;Wallace et al., 2012). Previous studies confirmed that the hypocenters and focal mechanisms of slow earthquakes are consistent with shear slip along the plate boundaries. However, the relationship between slow earthquakes and the neighboring seismogenic zone has not yet been fully understood.
The Philippine Sea Plate and Pacific Plate are subducting beneath the island arc around Japan (Figure 1). The characteristics of the subducting plates completely differ. The Philippine Sea Plate subducting in the Nankai Trough is young and warm, whereas the Pacific Plate subducting in the Japan and Kuril trenches is old and cold (e.g., Syracuse et al., 2010). The plate convergence rates of these plates also differ, 4-5 and 8-9 cm/year in the Nankai Trough and in the Japan and Kuril trenches, respectively (DeMets et al., 1994). Despite these differences, both subduction zones have repeatedly experienced huge regular earthquakes.
Along the Nankai Trough, slow earthquakes occur in both the shallower and deeper extensions of the seismogenic zone. The characteristics of deep slow earthquakes have been extensively investigated using nationwide onshore seismic and geodetic networks (e.g., Ito et al., 2007;Obara, 2002;Obara & Ito, 2005). Shallow slow earthquakes along the Nankai Trough have been investigated using both onshore and offshore seismic records (e.g., Asano et al., 2008;Nakano et al., 2018;Obara & Ito, 2005;Sugioka et al., 2012;. Recent studies revealed that the observed simultaneous occurrence of shallow tremors, VLFEs, and SSEs was similar to that of deep ETS Nakano et al., 2018) and that shallow slow earthquakes are activated by high pore fluid pressure in regions surrounding strongly locked zones .
Along the Japan and Kuril trenches, shallow VLFEs temporarily changed after the 2003 Tokachi-Oki and 2011 Tohoku earthquakes, respectively (Asano et al., 2008;Matsuzawa et al., 2015). In recent studies based on onshore and offshore data, more shallow tremors and VLFEs were detected (Baba et al., 2020;Nishikawa et al., 2019;Tanaka et al., 2019). Results suggested that the slow earthquake activity and large coseismic slip area of a huge earthquake are separated in the along-strike direction. Although the relationships between both types of earthquakes have been extensively investigated in both subduction zones, differences in the spatiotemporal variation of the slow earthquake activity between both subduction zones have not been discussed in detail. This difference may be related to the activity of huge earthquakes or the stress state of the plate boundaries.
Slow earthquakes are inhomogeneously distributed at the plate boundary (Obara & Kato, 2016). Ujiie et al. (2018) identified in a geological study that tremors occur in complex structures, such as subduction  (Koketsu et al., 2012). The black arrow in Figures 1a and 1b indicates the convergence direction of the Philippine Sea Plate, which subducts below the Eurasian Plate in the Nankai Trough, and that of the Pacific Plate, which subducts underneath the North American Plate in the Japan and Kuril trenches, respectively. mélange. Therefore, the spatial variation of their activity can reflect the heterogeneity of the frictional conditions in the fault zone, which can describe a host of different physical properties. Investigations of the activity of slow earthquakes within the subduction zones can provide new insights into the stress accumulation or frictional conditions at the plate boundary. To compare VLFE activities across Japan, we comprehensively detected VLFEs in Southwest Japan using the same method as in our previous studies (Baba et al., 2018(Baba et al., , 2020, which elucidated the distribution of deep VLFEs in Southwest Japan and shallow VLFEs along the Japan and Kuril trenches. The VLFEs were detected using decade-scale onshore seismic records. These records were used because shallow VLFEs can be detected due to the effective propagation of surface waves, the observation period of onshore networks is longer than that of offshore networks, and the comparison of deep and shallow VLFEs is possible using the same data set. Based on the newly constructed catalog, we discussed the characteristics of regions with slow earthquake activity from geodetic and geophysical viewpoints. The shear stress is accumulated at the plate boundary as a result of interplate locking. In addition, the presence of pore fluid, which decreases the seismic velocity, can change the frictional conditions of the plate boundary. Therefore, the comparisons of the VLFE activity with the slip-deficit rate and seismic velocity structure provide insights into the mechanical properties at the plate boundary.

Data
We used continuous seismograms of F-net broadband seismometers (Okada et al., 2004) from January 2003 to June 2019 after removing instrumental responses and resampling at one sample per second. A bandpass filter with a frequency range of 0.02-0.05 Hz was applied to all seismograms to enhance VLFE signals of onshore seismic stations (e.g., Ghosh et al., 2015;Ito et al., 2009).

Detection of VLFEs
Generally, the detection procedure used for VLFEs was the same as that reported in Baba et al. (2020). We placed 196 virtual epicentral grids on the Philippine Sea Plate boundary in Southwest Japan ( Figure S1 in the supporting information) in intervals of 0.3°and computed synthetic waveforms for the 10 stations closest to each virtual source grid using the open-source finite difference method code (OpenSWPC; Maeda et al., 2017) and by using a three-dimensional velocity structure model of the Japan Integrated Velocity Structure Model (JIVSM; Koketsu et al., 2012). We computed waveforms on a 3-D grid with a spacing of 0.2 by 0.2 km. We used the Küpper wavelet with a duration of 10 s and M w of 4.0 as a source time function. The focal mechanisms were assumed to be consistent with the geometry of the plate boundary of the JIVSM and plate motion model, NUVEL-1A (DeMets et al., 1994). We then calculated cross-correlation coefficients between the filtered synthetic template waveforms and F-net seismograms every 1 s. We selected events with the stationand component-averaged coefficients exceeding the threshold defined as nine times the median absolute deviation of the distributions. If the distribution is Gaussian case, 9 × MAD is 6.1 times the standard deviation. The probability of exceeding 9 × MAD for a normally distributed random variable is~6.4 × 10 −10 (Shelly et al., 2007).
False detections by regional regular earthquakes were removed using the catalog of the Japan Meteorological Agency. We removed the events detected between the P wave arrivals of teleseismic events whose M w was more than 4.5 listed in the catalog of the U.S. Geological Survey and 600 s after the S wave arrivals. However, considerable false detections remained, even after removing the teleseismic events based on the catalogs. Although the event amplitudes and cross-correlation coefficients generally are positively correlated (Baba et al., 2020), events with high amplitudes and average cross-correlation coefficients occur, which are considered to be false detections that are mainly caused by teleseismic events. Therefore, we did not count events with average cross-correlation coefficients below 0.4 and relative amplitudes to templates higher than 0.2 or average cross-correlation coefficients below 0.38 and relative amplitudes to templates higher than 0.1, except for Hyuga-nada. In the Hyuga-nada region (south of 32°N in the study area), the events had average cross-correlation coefficients below 0.4 and relative amplitudes to templates higher than 0.8. We established different thresholds for Hyuga-nada because typical VLFE amplitudes are larger than those in other areas.

Estimation of the Moments of Events
We calculated the relative amplitude of an event to synthetic waveforms with source durations of 10 s and M w 4.0 (c): where f i (t) and g i (t) are the observed waveform and synthetic template waveform at the ith station and jth component, respectively. The relative amplitude c was calculated to minimize the variance reduction between the synthetic template waveform and the observed waveform. The moment of each event (M o event ) was estimated from the amplitude of the event relative to the template: where M o syn is the moment of the synthetic waveforms of M w 4.0. Subsequently, we estimated the VLFE magnitude (M event ) using the following relationship between magnitude and moment (Hanks & Kanamori, 1979): The frequency distribution of VLFEs is shown in Figure S2. When we estimated the magnitudes of VLFEs, we excluded the virtual epicentral grids with the number of detected events below 35 or in which most of the events were falsely detected, mainly due to the teleseismic events that remained after discarding false detections using the process described above. The ratio of false detections was examined by visually investigating the detected event waveforms. Although many events were detected near the coast of Kyushu, most of them were false detections. The tendencies of the estimated moment density release rates off Cape Muroto, off the southern Kii Peninsula, and off the southeastern Kii Peninsula are similar to those reported in a previous study .
Regarding the VLFEs along the Japan and Kuril trenches, we evaluated the magnitudes of VLFEs detected in Baba et al. (2020) based on Equations 1-3. The cumulative moment of each grid was calculated using the sum of moments of each VLFE. Error estimation of the cumulative moment of each grid evaluated by using the nonparametric bootstrap method (e.g., Tichelaar & Ruff, 1989) is described in Text S1.

Results
The VLFEs along the Nankai Trough are distributed in the depth ranges of 30-40 km (deep VLFEs) and 5-10 km (shallow VLFEs; Figure S3). Depths of shallow VLFEs along the Nankai Trough are slightly less than those along the Japan and Kuril trenches (10-30 km; Baba et al., 2020). The different source depths of shallow VLFEs in these subduction zones commonly give a similar temperature range (Saffer & Wallace, 2015). We classified deep VLFE activity into four regions (i.e., western Shikoku, eastern Shikoku, Kii Peninsula, and Tokai) and shallow VLFE activity into four regions (i.e., Hyuga-nada, off Cape Muroto, off the southern Kii Peninsula, and off the southeastern Kii Peninsula) according to their spatiotemporal characteristics ( Figure S3).
The number of deep VLFEs detected in western Shikoku, eastern Shikoku, the Kii Peninsula, and Tokai is 895, 243, 594, and 193, respectively ( Figure S4a), whereas the number of shallow VLFEs detected in Hyuga-nada, off Cape Muroto, off the southern Kii Peninsula, and off the southeastern Kii Peninsula regions is 15,249, 1,123, 168, and 1,758, respectively ( Figure S4b). To discuss the relationship between the VLFE activity and interplate coupling, we estimated the cumulative moment of VLFEs for each grid. The temporal change of cumulative moment calculated by the sum of seismic moments of each VLFE, which was estimated using the amplitude magnitudes (details were described in method), yielded results similar to the temporal change of the total number of VLFEs (Figures 2a, 2c, and 2d). The cumulative moment of VLFEs along the Japan and Kuril trenches detected in Baba et al. (2020) was also estimated (Figure 2b). Hirose . The spatial variation of the VLFE activity rate was evaluated using the cumulative moment density release rate, which is obtained by dividing the cumulative moment of the detected VLFEs in each grid by the analysis period and grid area. The rapid increases in the cumulative moment of shallow VLFEs off Cape Muroto can be modulated by shallow SSEs off the Kii channel (Yokota & Ishikawa, 2020) in 2009 (M w 6.2) and 2018 (M w 6.6; Figures 2a and 2d), which were detected by geodetic observation. The temporal correlation between slow earthquake activity and geodetic data was also shown in the Cascadia subduction zone (e.g., Rouet-Leduc et al., 2019). The intervals of VLFE activations are longer for shallow VLFEs than for deep VLFEs, and, unlike deep VLFE activity, shallow VLFE activity has no regular periodicity (Figure 2d).

Geophysical Research Letters
The moment density release rate of deep VLFEs and its spatial variation is smaller than those of shallow VLFEs (Figure 2a). The along-strike spatial pattern of deep VLFEs is generally consistent with the distribution of energy released by deep tremors (Annoura et al., 2016). On the other hand, shallow VLFE activity shows a strong spatial heterogeneity along the Nankai Trough (Figure 2a). The largest moment density release rate was observed in the Hyuga-nada region in which earthquakes with M w > 8 have not been recorded.

Correlation Between the VLFE Activity and Interplate Coupling
The temporal changes in the shallow VLFE activity are synchronous with the interplate coupling change after huge regular earthquakes. To compare the VLFE activity with the interplate coupling in both subduction zones, we determined the coupling ratio by dividing the slip-deficit rate of each grid by the maximum slip-deficit rate in each subduction zone, assuming that the interplate coupling is 100% at the location of the maximum slip deficit (Hashimoto et al., 2012;Noda et al., 2018). The spatial variation of slip-deficit rates was derived from inversions of interseismic GNSS displacement rate data. Along the Japan Trench, the moment density release rate based on VLFEs increased off Ibaraki and off Iwate regions and decreased off

Geophysical Research Letters
Fukushima and off Miyagi regions since the 2011 Tohoku earthquake ( Figure S5a). In addition, a M w 8 earthquake occurred in the off Tokachi region along the Kuril Trench at the beginning of the analysis period, and it has not been confirmed whether the interplate locking has been fully recovered or not (Itoh et al., 2019;Nomura et al., 2017). The moment density release rate off Tokachi continued to decrease until 2013 ( Figure S5b). This tendency may indicate the recovery of the interplate locking around the coseismic slip region.
The strong spatial heterogeneity of shallow VLFE activity correlates well with the spatial distribution of the interseismic slip-deficit rate (Hashimoto et al., 2012;Noda et al., 2018) along the plate boundary (Figures 2a  and 2b). The regions with a high slip-deficit rate and those with VLFE activity are separated, and VLFE activity is typically concentrated in regions surrounding areas with a high slip-deficit rate in both subduction zones. To compare the VLFE activity in preparation for the next huge earthquake, we use VLFEs along the Nankai Trough and VLFEs off Tohoku only before the 2011 Tohoku earthquake. The moment density release rate of shallow VLFEs and the coupling ratio are negatively correlated (Figure 3a). The cross-correlation coefficient between the common logarithm of the moment density release rate and the coupling ratio is −0.44 ± 0.14. Within huge earthquake (strong interplate coupling) areas, such as Nankai (off Muroto, off the southern Kii Peninsula, and off the southeastern Kii Peninsula) and off Tohoku (off Iwate, off Miyagi, off Fukushima, and off Ibaraki), the moment density release rate of shallow VLFEs is low (Figure 3a). In contrast, the coupling ratio in Hyuga-nada is low compared with that of other shallow VLFE regions, and the moment density release rate is the largest. This negative correlation suggested that accumulated stress is released frequently by slow earthquakes in the weakly coupled areas. In contrast, the contribution to stress release by huge earthquakes is large, and slow earthquakes release a small fraction in the strongly coupled areas. Some geodetic studies discuss what percentage of the slip deficits is released by slow earthquakes. Most of the slip deficits beneath the Bungo channel are released by long-term SSEs (Noda et al., 2018), whereas the contribution of SSEs off the Kii peninsula to the release of slip deficits is 30-50% . In some regions off Tohoku, the interplate coupling is strong, but the moment density release rate is relatively high. In 2008, M w 6-7 interplate earthquakes (Nomura et al., 2017) occurred off Fukushima and off Ibaraki regions, which might have activated VLFEs. Because of this triggering process, the negative correlation between the interplate coupling ratio and VLFE activity may be unclear off Tohoku regions. After huge earthquakes, there is no correlation between moment density release rate and interplate coupling ratio ( Figure S6).
In Ecuador, huge earthquakes occur in strongly coupled areas, whereas SSEs release accumulated stress in weakly coupled areas in which no huge earthquakes have been recorded (Vaca et al., 2018). SSEs in the Cascadia subduction zone also occur in the area with weak interplate coupling (e.g., Michel et al., 2019). These tendencies are the same as that in Japan: Accumulated stress can be frequently released by slow earthquakes in weakly coupled areas, whereas the proportion of slow earthquakes is relatively small in terms of stress release in strongly coupled areas. In other words, slow earthquake activity is probably related to the coupling ratio.
On the other hand, deep VLFEs occur only in areas with weak interplate coupling, and the moment density release rate and its variation are small (Figure 2a). Thus, there are no meaningful spatial relationships between the moment density release rate of deep VLFEs and the coupling ratio ( Figure 3b). In areas in which deep VLFEs occur, the proportion of the release of the accumulated stress by deep VLFEs may not be as large as that of shallow VLFEs. The annual slip rate of short-term SSEs associated with ETS in Southwest Japan that was previously estimated to be 2-4 cm/year (Hirose & Obara, 2006) by previous studies is approximately half of the convergence rate of the Philippine Sea Plate. The geodetically estimated weak coupling and small moment density release rate of VLFEs might be affected by such decoupling properties at the plate boundaries in deep slow earthquake source regions next to a stable sliding zone.

VLFE Activity and Seismic Velocity Structure
Based on the comparison between shallow VLFE activity and seismic wave velocity variation along the Nankai Trough (Wang & Zhao, 2006;Yamamoto et al., 2017), shallow VLFEs are mainly distributed within low-velocity anomalies of the bottom of the overriding plate. This tendency is similar to that reported in previous studies (Kitajima & Saffer, 2012;Tonegawa et al., 2017). As for the Japan Trench, there is a high P wave velocity (Vp) area at the bottom of the hanging wall (Zhao et al., 2011).
This high Vp area corresponds to the coseismic slip area of the 2011 Tohoku earthquake; low Vp areas can be observed north and south of the high Vp area (Zhao et al., 2011). These areas correspond to areas with VLFE activity (Baba et al., 2020), such as off Iwate, off Fukushima, and off Ibaraki regions. Within the largest coseismic slip area and at the plate boundary deeper than 35 km, Vp is high, and there are few VLFEs, which was also indicated by the tremor activity (Nishikawa et al., 2019). Although the resolution of S wave tomography is generally lower than that of P wave, the distribution of low S wave velocity areas is almost the same as that of low Vp areas (Huang & Zhao, 2013;Yamamoto et al., 2017).
The existence of low-velocity areas suggests a high pore fluid pressure due to fluids dehydrated from subducting slab (Kamei et al., 2012;Tonegawa et al., 2017). The decrease in the effective normal stress due to the high pore pressure reduces the frictional strength at the plate boundary, which triggers the generation of VLFEs with a low stress drop (Ito & Obara, 2006;Saffer & Wallace, 2015). Undrained conditions, in which much water is contained, could be developed within such regions, similar to the fault planes of deep slow earthquakes (Nakajima & Hasegawa, 2016). The VLFEs may be considered to be indicators of interplate slip delineating the firmly locked portion. Off Aomori and off Tokachi regions, VLFEs actively occur, but the Vp is high. In the region between the Japan and Kuril trenches, regular earthquakes are rare. In addition, the afterslip of the 2003 Tokachi-Oki earthquake has continued in this region, indicating that there might be another factor activating VLFEs.

Mechanical Properties of Regions With VLFE Activity
The VLFEs occur adjacent to large coseismic slip areas of huge earthquakes in both subduction zones of Japan. In the shallow portion, the moment density release rate of VLFEs and geodetically estimated coupling ratio at the plate boundary are negatively correlated (Figure 3a). In strongly coupled areas, which correspond to the largest coseismic slip areas of huge earthquakes, the interplate frictional strength is high, which can explain the occurrence of high-speed ruptures. Because the effective strength of the plate boundary may be high in such rupture areas, the slow earthquake activity rate is low. On the other hand, in weakly coupled areas, the accumulated stress is frequently released by slow earthquakes, and huge earthquake nucleation cannot be favorably initiated.
Although there are a few exceptions, the shallow VLFE activity tends to be high in areas with relatively weak interplate coupling and low seismic velocities. In areas with weak frictional conditions, shallow VLFEs can be activated by the decrease in the effective normal stress due to the high pore fluid pressure. On the other hand, the variation in the moment density release rate of deep VLFEs is smaller than that of shallow VLFEs. This suggests that the horizontal heterogeneity of the frictional properties is stronger in the shallow part of the plate boundary near the seismogenic zone than in the deep part of the plate boundary. The differences in thermal conditions between the deep and shallow parts of the plate boundary (e.g., Saffer & Wallace, 2015) or the complex bathymetry of the plate interface in the shallow part (Okino et al., 1999), such as the existence of seamounts, may be the cause of strong heterogeneity of shallow slow earthquake activity. Hutchison (2020) pointed out that physical complexities can control the distribution slow earthquake activity; therefore, the physical complexity may be stronger in the shallow part than in the deep part.
The results of previous studies (Saffer & Wallace, 2015; strongly suggested that the presence of pore fluid can control the slow earthquake activity. In Baba et al. (2020), we clarified the difference in the VLFE activity inside and outside the largest coseismic slip area of the 2011 Tohoku earthquake. In this study, we suggest that the VLFE activity is strongly related to the distribution of both the interplate coupling and pore fluid, which reflect the frictional properties at the plate boundary. The temporal changes in the VLFE activity are synchronous with the interplate coupling change due to huge earthquakes. Therefore, VLFE activity can reflect the spatiotemporal variation of interplate coupling in subduction zones.

Data Availability Statement
The coseismic slip data of the 1944 Tonankai earthquake, 1946 Nankai earthquake, and 2003 Tokachi-Oki earthquake were derived from the Finite-Source Rupture Model Database (Mai & Thingbaijam, 2014; http://equake-rc.info/SRCMOD/). Catalogs of slow earthquakes along the Nankai Trough and the Japan Trench were downloaded from the Slow Earthquake Database (