Seismicity distribution in the Tonankai and Nankai seismogenic zones and its spatiotemporal relationship with interplate coupling and slow earthquakes

We conducted seismic tomography to estimate the seismic velocity structure and to evaluate the spatiotemporal distribution of interplate earthquakes of the Kii Peninsula, central Honshu, Japan, where the Tonankai and Nankai megathrusts are located. Microearthquakes were quantitatively detected by using the data from a cable-type seafloor seismic observation network, completed in 2015. Our velocity model was consistent with the previous 2-D active-source surveys, which reported the areal extent of key structures: a high-velocity zone beneath Cape Shionomisaki, a subducted seamount off Cape Muroto, and the subducted Paleo-Zenith Ridge. The absence of any other subducted seamount with the same or larger spatial scale, than the identified key structures, was confirmed. Our velocity model also revealed that there was not a simple relationship between areas of large coseismic slip or strong interplate coupling and areas of high velocity in the overriding plate. Relocated hypocenters widely ranged from the upper plate to within the slab, while the most active region was attributed to the oceanic crust in the aftershock region of 2004 off-Kii earthquake. Compared with the results from the land-based observation network, the accuracy of the focal depth estimation was substantially improved. Furthermore, we identified the seismic activity in the vicinity of the plate boundary and determined 14 locations for interplate seismicity areas. They were primarily distributed in the range of seismogenic zone temperature (150–350 °C) along the plate boundary and were located outside of the strong interplate coupling zone. Several active areas of interplate earthquakes exhibited clustered activity during the periods of slow-slip events, observed and accompanied with shallow very-low-frequency earthquakes. Thus, regular interplate microearthquakes became active at the plate boundary in the conjunction with slow slip. In summary, as regular earthquakes provide a more accurate source location than slow earthquakes and can detect events of smaller magnitude, monitoring such interplate earthquakes may reveal spatiotemporal variations in the stick–slip conditions on the plate boundary.


Introduction
The plate boundary between the subducted Philippine Sea Plate and overriding Amur Plate is located in the Nankai Trough subduction zone, central Japan (Bird 2003;Argus et al. 2011). The Nankai trough subduction zone is well known for repeated M 8 class megathrust earthquakes, which have occurred at intervals of
The Japanese government has stated that the possibility for the occurrence of M > 8-9 class earthquake within 30 years is as high as 70-80% (The Headquarters for Earthquake Research Promotion 2022). As a result, various seismological and geodetic studies have been conducted on the Nankai Trough subduction zone: For example, a number of these studies demonstrated the subsurface structural heterogeneities that can control rupture propagation (Kodaira et al. 2002(Kodaira et al. , 2006

K i i P e n i n s u l a Cape Muroto
T o n a n k a i T o n a n k a i N a n k a i N a n k a i Tokai Tokai Fig. 1 Maps of the study area. The inset regional map shows the location of the study area with the megathrust seismogenic segments (Tokai, Tonankai, and Nankai) of the Nankai Trough in southwestern Japan with plate boundaries. Two letter abbreviations reflect the name of the Plate; AM: Amur Plate, PS: Philippine Sea Plate, OK: Okhotsk Plate, and ON: Okinawa Plate (Bird 2003). Yellow allows are the relative motion of the PS with AM (Argus et al. 2011). In the main topographic map, the megathrust seismogenic zones along the Nankai subduction zone are shown by the red line. The Shionomisaki igneous complex (SIC; Kodaira et al. 2006), subducted seamount identified off Cape Muroto (SSM; Yamamoto et al. 2017), and subducted Paleo-Zenith Ridge (PZR; Park et al. 2004) are outlined in black. Blue diamonds show the location of the seismic stations of the Dense Oceanfloor Network system for Earthquakes and Tsunamis (DONET). The epicenters of the 1944 Tonankai and 1946 Nankai earthquakes were taken from Kanamori (1972) (yellow stars), and two M > 7 earthquakes of the Kii Peninsula occurring on September 5, 2004 (red stars) and the M JMA 6.5 off-Mie earthquake occurring on April 1, 2016 (orange star) are from the JMA unified catalog and indicated the location of the strong coupling regions (e.g., Yokota et al. 2016;Nishimura et al. 2018;Noda et al. 2018).
Moreover, various types of slow earthquakes, longterm and short-term slow-slip events (SSE), very-low-frequency earthquakes (VLFE), tremors, and low-frequency earthquakes have been identified in both the downdip and updip sides of the Nankai Trough seismogenic zone (Obara 2002;Ito et al. 2007;Obara and Kato 2016;Nakano et al. 2018a). The spatiotemporal variation of slow earthquake activity can be expressed as an index of the stress perturbation around the seismogenic zone. For instance, simulation studies have shown that the intervals between the occurrence of slow-slip events along the downdip side of the seismogenic zone are shortened once the onset of a large earthquake approaches . Thus, slow earthquakes are being thoroughly monitored to estimate the status and future temporal evolution of the stick-slip behaviors of the Nankai seismogenic zone by the Japan Meteorological Agency (JMA).
It is essential to monitor the spatiotemporal distribution of both slow earthquakes and regular interplate earthquakes for determining the occurrence processes of large earthquakes. However, studies on regular earthquake activity are rare because of their low activity level. Moreover, the accurate determination of the hypocenter location is hampered by the offshore conditions. Unlike the Japan and Ryukyu Trench forearcs, small repeating earthquakes along the plate interface are also inactive in the Nankai Trough forearc (Igarashi 2020). On the one hand, previous offshore seismic observation studies indicated that most of the earthquakes in the Nankai forearc are intraplate events, located within the Philippine Sea slab or overriding Amur Plate (Mochizuki et al. 2010;Yamamoto et al. 2017). On the other hand, based on the offshore seismic observation data, the Mw 5.9 (M JMA 6.5) off-Mie earthquake, which occurred on April 1, 2016, was classified as an interplate earthquake (Nakano et al. 2018b) (Fig. 1). The occurrence of this event indicated the existence of other interplate microearthquakes that could not be detected or have been plagued by in-precise location detection by onshore seismic observation.
From 2010, the Dense Oceanfloor Network system for Earthquakes and Tsunamis (DONET1 and 2; Kawaguchi et al. 2015;Kaneda et al. 2015) was started to establish ( Fig. 1). DONET1 and DONET2 have been fully operational since 2012 and 2015, respectively. Broadband seismometers were installed at 51 stations in total, based on which continuous seismic records were obtained. These instruments exhibit superior detection capabilities compared with the land-based observation network. Moreover, offshore observation data can improve the accuracy in epicenter location. However, it remains crucial to identify the effect of velocity structure heterogeneity for accurate epicenter determination. Although previous amphibious seismic tomography based on the offshore temporal seismic observation and active source survey established a three-dimensional (3-D) velocity structure model in this region (Yamamoto et al. 2017), they used only a limited portion of DONET1 stations (Fig. 2a). A recent amphibious tomographic study based on large datasets of both active and passive sources also established a 3-D velocity structure for the entire Nankai Trough seismogenic zone (Arnulf et al. 2022), but they analyzed only P-wave structures and used no offshore seismic stations for passive source data.
Thus, the main aim of this study was to investigate the seismic activity of interplate earthquakes by using seismic data from the DONET. To achieve high accuracy of the estimated hypocenter location, we performed seismic tomography by adding DONET data (Fig. 2b) to a previous 3-D tomographic study (Yamamoto et al. 2017) and updated the velocity structure and hypocenter location simultaneously. Then, we discussed the spatiotemporal patterns of microseismicity around the plate boundary and their relationship with the areal extent of structural heterogeneities, interplate coupling, and slow earthquakes.

Methods
We established the dataset for seismic tomography in this study by combining data from Yamamoto et al. (2017) with the additional first-arrival data of microseismicity obtained by DONET. We extracted 1999 relocated events and 465 seismic stations as the initial passive source data from Yamamoto et al. (2017) (Fig. 2a). P-wave arrivaltime data from previous 11,168 active-source shots along 13 lines at 545 stations were also used. Following the methods of Yamamoto et al. (2017), we used the time difference between S-and P-wave arrivals as input data for two of our ocean bottom seismographs in the off Kumano region that could not be calibrated with Global Navigation Satellite System clock information.
For additional data, we first detected the event by using the STA/LTA method (the ratio between short-time average and long-time average of the waveform amplitude). Both the WIN system (Urabe and Tsukada 1992) and the method of Horiuchi et al. (2009) were applied for all DONET stations (Fig. 2b). Second, we combined these results and deleted duplicates. We manually picked first-arrival times of the detected events on the continuous seismic records at both selected land and all DONET stations. Owing to the processing, the dataset for events occurring from September 2015 to March 2016 and after April 2019 was incomplete and partially analyzed. As a result, we detected 21,368 events in total from May 2012 to February 2020. Note that the magnitude for each event was quantified by the using maximum amplitude of a vertical component seismogram (Watanabe 1971).
Then, we conducted hypocenter relocation using the 3-D velocity structure of Yamamoto et al. (2017) to identify the initial location of an additional dataset for tomographic analysis. This procedure was conducted by using the tomoFDD code (Zhang and Thurber 2006). For S-wave first arrivals of the additional dataset, we added a station correction (ΔT S ), calculated from the traveltime differences between the P-to-S converted phases at the basement of the sediment layer and direct P-waves (T PS-P ). This procedure was based on Yamamoto et al. (2017).
At the next step, we selected the dataset analyzed in this study using the same criteria from Yamamoto et al. (2017): (1) At least six first-arrival times of P-or S-waves existed, and (2) either a gap in azimuthal coverage was < 180 degrees or the minimum epicentral distance to the nearest station with phase picking was < 30 km. Overall, we obtained 16,042 additional earthquakes for use in the tomographic study (Fig. 2b). As a result, we used a total of 11,168 active source and 18,041 earthquakes at 715 stations (Fig. 3). By utilizing these data, we conducted double-difference seismic tomography. We further applied the tomoFDD package (Zhang and Thurber 2006) and accepted the estimated 3-D velocity model from the previous study (Yamamoto et al. 2017) as an initial 3-D velocity structure. As in a previous study, we set velocity grid nodes at 10 km intervals in the horizontal direction and 2-10 km intervals in the vertical direction. We calculated double-difference data for event pairs whose epicentral distances were < 10 km. As a result, we acquired 423,050 P-wave arrivals, 414,811 S-wave arrivals, 141 S-P intervals, and double-difference data for 2,109,460 P-waves and 1,673,567 S-waves.
We determined the damping parameter for tomographic inversion by considering the recovered amplitudes of a checkerboard resolution test (CRT) (Fig. 4), while the weights and smoothing parameters were established with the same values as in Yamamoto et al. (2017). The sizes of the checkerboard pattern were 20 km, 40 km, and 10-20 km in trough-normal (X-axis), trough-parallel (Y-axis), and vertical directions (Fig. 4f ). We assumed that the amplitude of perturbation was ± 3%. We added 0.1 and 0.2 s of Gaussian noise for the P-and S-waves, respectively, to synthetic travel times. We further compared the CRT results with the spatial distribution of the derivative weight sum (DWS; Thurber and Eberhart-Phillips 1999), and defined the resolved area based on the following criteria: (1) a recovery rate of CRT, determined by dividing the result of the CRT amplitude by the assumed CRT amplitude, of > 0, and (2) the DWS > 1000 for P-waves and > 500 for S-waves.

Results
After 10 iterations, the root mean square of the absolute travel-time residuals decreased from 0.28 to 0.22 s for P-wave arrivals and from 0.48 to 0.36 s for S-wave arrivals. As a result, all 18,041 earthquakes were relocated. The average errors of hypocenter relocation were found to be 0.3 km in both the horizontal and vertical directions. Figure 5 shows the relocated hypocenters. As seen, the offshore earthquakes were distributed to ~ 40 km depth. To investigate the depth distribution relative to the plate interface, we applied the plate geometry mode of Nakanishi et al (2018). Although there are several other plate geometry models (e.g., Hashimoto et al. 2004;Hirose et al. 2008;Hayes et al. 2018), the model of Nakanishi et al. (2018) included the largest number of offshore seismic surveys; thus, we consider their model to be the most suitable. The result for all the relocated events indicates that most of these events were intraslab events, while the major active areas corresponded to oceanic crust (Fig. 5b). However, most events within the oceanic crust were concentrated within the aftershock region of the 2004 off-Kii earthquake. In addition to the intraslab events, we identified events that occurred close to the plate interface.
Next, we compared our relocation results with the JMA unified catalog (Japan Meteorological Agency 2022). Figure 6a shows both horizontal and vertical differences for M JMA > 2 events. We did not reveal any large differences between our findings and the JMA catalog, despite several events exhibiting ~ 10 km difference. However, the offshore relocated hypocenters were generally 10-15 km shallower than the focal depths of the JMA catalog. These findings are similar to the trends reported for the aftershock of 2004 off-Kii earthquakes based on offshore seismic observations (Sakai et al. 2005;Nakano et al. 2015). Thus, retrieval of precise focal depth by using only onshore seismic network is hampered. We also compared Distance along X-axis (km) Distance along X-axis (km) Distance along X-axis (km) the estimated magnitude in this study and that of the JMA (Fig. 6b). Our magnitude was slightly larger than that of the JMA catalog in general. Finally, we examined the distribution of the number of events by magnitude (Fig. 6c). The number of events was evaluated in increments of 0.1, based on the magnitudes from this study. We found that all events above M2.8 were those listed in the JMA catalog, whereas the most of earthquakes whose magnitude was < 2 were not. Thus, the detection capability was strongly improved by the offshore observation. The estimated velocity model clearly illustrates the subducted Philippine Sea plate down to 40 km depth (Fig. 7). The comparison of these results with the plate interface model of Nakanishi et al (2018) revealed the high velocity slab mantle (Vp > 8.0 km/s, Vs > 4.5 km/s) at ~ 10 km deeper than the plate interface. Moreover, several lowvelocity anomaly zones within the slab were identified.
The most significant zone was related to the Nankai segment (gray arrow in Fig. 7a; Y = −80 km, X = 20-50 km) for both the P-and S-wave models. We also found that the high-velocity zone just above the plate interface around the toe of the Kii Peninsula (blue arrow in Fig. 7b, X = 0 km, Y = −20 to 50 km; X = 30 km, Y = −10 to 30 km).
The analysis of the velocity structure along the plate boundary ( Fig. 8) revealed that the Vp = 5 km/s and Vs = 3 km/s iso-velocity contours were nearly parallel to the trough axis. However, larger iso-velocity contours (e.g., Vp = 6 km/s; Vs = 3.5 km/s) tended to be convex to the south, centered at the tip of the Kii Peninsula. The earthquakes within 2.5 km of the plate boundary were clustered in several locations and distributed in the area from near the trough axis to ~ 20 km depth at the plate boundary. A high velocity region with Vp > 8 km/s and Vs > 4.25 km/s was imaged beneath the Kii Peninsula (34°N,135.5°E), and deep low-frequency tremors (Maeda and Obara 2009;Obara et al. 2010) were distributed on the updip side of this region.

Location and spatial extent of structural anomalies
First, we evaluated how the structural heterogeneities, estimated by previous studies, were reconstructed by our velocity model. In the overriding Amur Plate, high-density and high-seismic velocity body called the Shionomisaki igneous complex (SIC) exists just beneath the toe of the Kii Peninsula (Honda and Kono 2005;Kodaira et al. 2006;Qin et al. 2021). In both the average P-and S-velocity perturbations within a 5-km-thick layer just above the plate interface ( Fig. 9a), the high-velocity zone is clearly imaged. In particular, the center of the high-velocity zone (A in Fig. 9a), located just north of the epicenter of the 1946 Nankai earthquake, was attributed to the highest gravity anomaly (Honda and Kono 2005). This finding leads to the conclusion that this high-velocity region represents the spatial extent of the SIC.
Based on the interpretation of Kimura et al (2014), this high-velocity region reaches the vicinity of the wedge mantle. The high-velocity zone between 30 and 40 km of the depth contour of the plate interface (B in Fig. 9a; 34° N, 135.5° E) corresponds to the landward mantle because the depth of Moho discontinuity in this area is shallower than 30 km (Matsubara et al. 2017). Although areas A and B were seemingly separated in the S-wave velocity model, their boundary was unclear in the P-wave velocity model (Fig. 9a). Arnulf et al. (2022) also investigated the spatial extent of SIC (called the "Kumano Pluton" in their paper) based on the high (> 6.5 km/s) P-wave velocity zone. The P-wave velocity contour of 6.5 km/s from our model at 15 km depth (green line in Fig. 10a) and the spatial extent  . Black, red, and gray correspond to all data, outside of the aftershock area denoted in a, and from Yamamoto et al. (2017), respectively of the Vp > 6.5 km/s areas within the overriding plate along Y-axis direction (− 30 ≤ X ≤ 40 km; black lines in Fig. 10a) are consistent with their result. Note that the Vp > 6.5 km/s areas of our model located south of the overlapped area with Arnulf et al. (2022) correspond to the oceanic crust of the subducting plate. Based on this consistency, we determined that the southern limit of the SIC was imaged by our model. However, our result is inconsistent with theirs at north of 34° N (Fig. 10a).
In addition, it is difficult to estimate the spatial extent of Vp > 6.5 km/s areas along the northern side (X ≤ −40) owing to the extent of the wedge mantle. Therefore, the north extent of the SIC was not determined in this study.
In certain subduction zones, the relationship between the structural heterogeneity of the upper plate and the seismic slip distribution has been previously investigated, thereby indicating that the high-velocity zone corresponds to the large coseismic slip zone (e.g., Zhao  For the area of interest in this study, we compared the structure with the distribution of interplate coupling strength (Yokota et al. 2016) and did not find a significant relationship, thereby agreeing with Yamamoto et al. (2017). Note there are various models for the slip distribution of the 1944 Tonankai and 1946 Nankai earthquakes whose characteristics were not similar. For example, the largest slip areas of the 1946 Nankai earthquake in Baba and Cummins (2005) do not overlap with those reported by Murotani et al. (2015). Thus, the comparison between our velocity structure and the slip distribution of past megathrusts was not conducted in this study.
Within the slab, the existence of subducted topographic anomalies, such as seamounts (Kodaira et al. 2000) and the Paleo-Zenith Ridge (Park et al. 2004), has been previously reported. We attempted to identify their location from our velocity model. As these features are characterized by thick crust, they are illustrated as lowvelocity regions in the deep structure of the slab relative to the surrounding mantle. Figure 9b shows the average velocity perturbation within a 5-km-thick layer between 10 and 15 km beneath the plate interface. We identified two low-velocity zones (C and D in Fig. 9b). Here, lowvelocity zone C corresponded to the subducted seamount (Kodaira et al. 2000), and its spatial extent was nearly the same as that from Yamamoto et al. (2017). Low-velocity zone D was located close to the area where the Paleo-Zenith Ridge subducted (Park et al. 2004), thereby indicating no overlap with the resolved area from Yamamoto et al. (2017). Although the location of low-velocity zone D from P-wave perturbations was slightly south of that of the S-wave perturbations, this low-velocity zone was attributed to the existence of the Paleo-Zenith Ridge. Although smaller-size subducted seamounts than our grid spacing can theoretically exist, we could not find other similar low-velocity zones in the study area. Thus, it is reasonable to suggest that there was seemingly no other subducted seamount or ridge as large as them. Arnulf et al. (2022) observe a significant (Vp = 6.5-7.5 km/s) low-velocity anomaly within the subducting mantle around the aftershock region of the 2004 off-Kii earthquake and interpreted it as a combination of serpentinization and enhanced porosity from bending stress. Although their model has the lowest point of P-wave velocity at 30 km depth, there is no such low-velocity zone in our model (Fig. 10b). Note that the low-velocity zones imaged between the 20 and 30 km contour of the plate interface depth correspond to the oceanic crust. Instead, the location of the significant low P-wave area (Vp < 6.5 km/s) from their model was spatially close to the subducted Paleo-Zenith Ridge, where we found the lowvelocity anomaly 10-15 km beneath the plate interface (Fig. 9b). That is, the depth of the low-velocity anomaly obtained in our study was ~ 10 km shallower than that of Arnulf et al. (2022). We consider this difference to be due to the poor hypocenter location accuracy in the dataset of Arnulf et al. (2022). As shown in Fig. 6, there is a large difference in hypocenter depth of offshore seismicity between this study and the JMA unified catalog, which Arnulf et al. (2022) applied. Because their analysis was performed by earthquakes whose calculated depths were   (Oleskevich et al. 1999). Other symbols are the same as Fig. 1a more than 10 km deeper than their actual depth, the lowvelocity region shown in Arnulf et al. (2022) was imaged at a deeper position than the actual location.

Spatial relationship among interplate earthquakes, interplate coupling, and slow earthquakes
We confirmed the occurrence of interplate earthquakes in the study area (Fig. 8). The comparison with the thermal structure (Oleskevich et al. 1999) revealed that most interplate earthquakes occurred between 150 and 350 °C, thereby resonating with the range of seismogenic zone depth (Hyndman et al. 1997). In addition, they were seemingly located outside of high Vp/Vs (> 1.9) zone except for two areas, where a subducted topographic high was identified (Fig. 8c). The downdip limit of the Vp/Vs > 1.9 zone was located ~ 100 km from the trough axis in line with the downdip limit of the dehydration occurred within the sediment layer (Hyndman and Peacock 2003). The high Vp/Vs ratio around the plate interface normally indicates the existence of pore fluids and is attributed to the weak coupling zone as pore fluids decrease the effective normal stress (Moreno et al. 2014). Thus, interplate seismicity can be interpreted as occurring in areas with relatively less pore fluid and not fully creeping zones. However, the high Vp/Vs zone in this study area contained at least one strong interplate

R e c u r r e n c e R e c u r r e n c e s h a ll o w S S E s s h a ll o w S S E s
Strong Strong coupling coupling zone zone Fig. 9 Velocity perturbation within the 5-km-thick layer. a At just above the plate interface and b between 10 and 15 km below the plate interface. White contours are depth to the plate boundary . Dashed contours reflect velocity differences from − 10 to + 10% at 2.5% increments. The areas of strong (> 5 cm/year) interplate coupling estimated from onshore and offshore geodetic data are outlined by purple lines (Yokota et al. 2016). Areas outlined by gray, green, and orange lines reflect the fault location of deep SSE beneath the Kii channel (Kobayashi 2014), shallow SSE between 2017 and 2018 (Yokota and Ishikawa 2020), and recurrence shallow SSE (Araki et al. 2017;Ariyoshi et al. 2021), respectively. A, B, C, and D are the locations of substantial velocity anomalies discussed in the main text. Other symbols are the same as in Fig. 1 coupling area (Yokota et al. 2016), indicating that the relationship between Vp/Vs and coupling strength along the plate boundary proposed by Moreno et al. (2014) cannot be applied to the Nankai Trough. Thus, the existence of interplate microearthquakes cannot be simply attributed to the Vp/Vs distribution. Rather, there are other potential factors, such as topographic anomalies in the subducted plate (Nakamura et al. 2022) or thickness variation of the subducted sediment layer (Akuhara et al. 2017), that may control the occurrence of interplate seismicity in the Nankai trough.
To identify the areas where interplate earthquakes prevailed, we selected the earthquakes within 1 km of the plate interface . On this basis, we defined 14 areas as active areas of interplate earthquakes (Fig. 11). These active areas were generally located outside the especially strong coupling zones (> 5 cm/yr; Yokota et al. 2016), except Area 11, where the Mw 5.9 off-Mie earthquake occurred on April 1, 2016 (Nakano et al. 2018b) (Fig. 11a). Two areas in the vicinity of the trough axis (Areas 13 and 14) corresponded to the aftershock area of the 2004 off southeastern Kii earthquakes.
At the next step, we examined the relationship between these 14 active areas and slow earthquakes. We did not find overlap with the SSE region in the Kii Channel (Kobayashi 2014). The segregation of slow earthquakes and regular earthquakes at the plate boundary has been previously demonstrated in the off Boso (Ito et al. 2019), Ryukyu Trench 2020), Hikurangi (Bartlow et al. 2014), and other areas. However, in the offshore, there was an area of overlap with the fault location of the SSE between 2017 and 2018 (Yokota and Ishikawa 2020) (Area 9; Fig. 11a). Moreover, three areas (Areas 8, 9 and 12) corresponded to the area where shallow VLFE was identified based on the offshore observation (Sugioka et al. 2012;Nakano et al. 2018a;Toh et al. 2020;Yamamoto et al. 2022). The overlapping areas were all located at the area, where the plate interface was located at the depth of ~ 7 km and was attributed to the subduction position of a seamount or its vicinity (Fig. 11a). The source location accuracy of VLFE was found to be lower (e.g., ~ 0.03° in horizontal and 1.3 km in depth for VLFE in Yamamoto et al. 2022) than that of the relocated hypocenters in this study (~ 0.3 km in both the horizontal and vertical direction in this study). Thus, a fine-scale analysis is challenging with the current methodology, but it is reasonable to suggest that several of the active areas of interplate earthquakes in the Nankai Trough subduction zone overlapped with the source areas of the slow earthquakes. Compared with the estimation of the VLFE location based on onshore data (Takemura et al. 2019), the active areas in Areas 6, 7, and 13 seemingly overlapped with VLFE (Fig. 11b). However, as the VLFE location, estimated by Takemura et al. (2019), was plagued by large uncertainty in the across-trough direction (see Fig. S1 of Takemura et al. 2019), their result may not reveal spatial differences between Areas 6 and 8, Areas 7 and 9, and Areas 12 and 13.
Both slow earthquakes and microearthquakes along the fault have been previously thought to occur in the transition zone between the strong coupling seismogenic zone and fully creeping zone (e.g., Sholtz 2002;Obara and Kato 2016). The magnitude of completeness for microearthquakes in this study is less than 2 (Fig. 6), which is much lower than that for VLFE (~ Mw = 3; Nakano et al. 2018a, b;  Yamamoto et al. 2022). In addition, the improvement of the accuracy in the hypocenter relocation using the offshore seismic network with a 3-D velocity model enabled the distinction between interplate and intraplate earthquakes. As both the detection capability and location accuracy of regular earthquakes were superior to those of shallow slow earthquakes, a comprehensive detection of interplate microearthquakes is valuable for monitoring the current stick-slip status of the strong coupling zone in the Nankai Trough.

Temporal variation of interplate earthquake activities
For each active zone, we examined the temporal changes in seismic activity and the number of integrated earthquakes at the plate boundary ± 25 km by using the entire relocation result from this study (Fig. 12). Because the network-based observations exhibited temporal variability, we focused on the temporal distribution after March 2015 at the area close to the DONET2 network (Areas 1-10) and after 2012 at the area within the DONET1 network (Areas 11-14), respectively. No remarkable changes were identified in Areas 1-5, which were close to land, or in Areas 10, 13, and 14, which were close to the trough axis. However, in other areas, characteristic activity was identified, even while considering the periods when the dataset was incomplete (gray masked periods in Fig. 12). In Areas 6, 8, and 9, the cumulative number of earthquakes increased sharply in 2018. Seafloor crustal displacement observations in this area suggest that SSEs occurred between 2017 and 2018 (Yokota and Ishikawa 2020). As the activations in these areas during the SSE period seemingly occurred between − 5 and 5 km in relative depth (Fig. 12), we examined the temporal variation of five areas close to this SSE region. Therefore, we narrowed the time window to two years and the epicenter to ± 5 km from the plate interface (Fig. 13a). As a result, swarm-like activation of interplate seismicity in Areas 8 and 9 was identified in March and May 2018, respectively.

Fig. 11
The location of the active areas for interplate seismicity and very-low-frequency earthquakes (VLFEs). Relocated hypocenters within the 1 km from the plate interface are shown by pink circles. a Numbers from 1 to 14 reflect area code. VLFEs estimated by using the offshore seismic network are shown by green circles (Sugioka et al. 2012;Nakano et al. 2018a;Toh et al. 2020;Yamamoto et al. 2022). Location of the subducted seamount off Muroto and the Paleo-Zenith Ridge is shown by light-blue areas. b VLFEs estimated by using onshore seismic network are shown by yellow circles (Takemura et al. 2019). Other symbols are the same as in Fig. 9 (See figure on next page.) Fig. 12 Temporal variations of seismicity in each area shown in Fig. 11. Blue dots reflect relocated hypocenters, and their vertical axis is relative depth from the plate interface. The red line shows the cumulative number of earthquakes within the relative depth range between − 25 and 25 km. Operation periods for each seismic observation are shown on the upper side of the panel by colored arrows; OBS1,2: ocean bottom observation data from Yamamoto et al. (2017), DONET1, and DONET2. Masked areas of gray, green, and orange reflect the period of incomplete dataset, shallow slow-slip events (SSEs) estimated by Yokota and Ishikawa (2020), and recurrence SSEs  Page 15 of 20 Yamamoto et al. Progress in Earth and Planetary Science 2022, 9( Although detailed discussion is not possible owing to the lack of an available VLFE catalog from offshore observation corresponding to this SSE, VLFE activity has been previously reported by Takemura et al. (2019), with most activity near the western side (Area 8) in March and the eastern side (Area 9) in May (black line in Fig. 13a).
However, the activation of interplate seismicity seen in Area 6 occurred approximately six months later than the timing of the VLFE activity observed by Takemura et al. (2019). Moreover, VLFE activity was not accompanied by the interplate seismicity in Area 7, and no swarm-like seismicity or VLFE activity was identified in Area 10. The onshore network-based VLFE activity periods were the same between Areas 6 and 8 and between Areas 7 and 9. This finding might indicate the previously indicated large location uncertainty in across-trough direction. Thu, it is necessary to determine the source locations of the VLFE by using the seafloor observation network during this activity period to determine the difference in activity between regular earthquakes and VLFE and to understand the state of stick-slip behavior in the shallow part of the subduction zone.
In Area 11, a sharp increase in the cumulative number of earthquakes related to the 2016 off-Mie earthquake was identified (Fig. 12). Unlike the other areas, this area was located within the strongly coupled zone (Yokota et al. 2016), where we also observed several seismic  activities occurring near the plate boundary before the 2016 earthquake (e.g., in 2013 and 2015). This indicates that the 2016 off-Mie earthquake did not occur at the seismic gap but in an area of background seismicity. The investigation of the location of background interplate seismicity within the strongly coupled area can unravel a potential nucleation point for moderate-size or, even large interplate earthquakes. Area 12 seemingly had no significant fluctuations during this period, except a slight increase in the cumulative number of earthquakes in the first half of 2018 (Fig. 12). However, a deeper investigation of Area 12 suggests that while intraslab earthquakes (> 5 km deeper than plate interface) generally prevail, earthquakes occur at relatively shallow depths close to the plate boundary (within 5 km from the plate interface) in four episodes during the observation period: the first half of 2014 (episode-1), the second half of 2015 (episode-2), the first half of 2017 (episode-3), and the first half of 2018 (episode 4) (Fig. 13b). The first and third episodes exhibited fewer earthquakes, while the second and fourth episodes experienced more earthquakes. Based on borehole pore pressure changes, recurrent SSEs were shown to occur in the area adjacent to the western region (Araki et al. 2017;Ariyoshi et al. 2021) (Fig. 11a). Of these events, SSEs with large amounts of slip were reported during episodes 2 and 4. VLFE activity was also confirmed by the offshore observation data during episode 2 (Nakano et al. 2018a), although no VLFE was observed during episode 4.
This study proved the existence of interplate earthquake swarms, synchronized with SSEs at Areas 8, 9, and 12 in the Nankai Trough. Such earthquake swarms, associated with SSE, have been previously reported in subduction zones such as off Boso (Hirose et al. 2014;Fukuda 2018), Ecuador (Vaca et al. 2018), and Hikurangi (Bartlow et al. 2014), where slow-slip stress loading and stress triggering outside the SSE region were the main drivers of earthquake swarms. In this study, the seismic swarm activity in Areas 8, 9, and 12 may have been activated by the SSE, because the timing of the SSE and VLFE activity was the same.
Near the subducted seamount in the Hikurangi forearc, similar regular earthquake swarms collocated with slow earthquakes were reported, but they were located within the overriding plate (Shaddox and Schwartz 2019). The observed waveforms near the swarm areas generally indicate an increasing difference between S-and P-wave arrivals with respect to hypocenter depth (Fig. 14). Considering this result and the accuracy of our hypocenter depth in this study (~ 0.3 km), we suggest that seismic swarms observed in this study are located within the slab and across the upper plate and include seismic activity at the plate boundary, rather than being located solely within the overriding plate as in Hikurangi. The focal mechanisms of low angle thrust solutions should be a strong piece of evidence for interplate seismicity. Although we could estimate the focal mechanism for the 2016 off-Mie earthquake based on P-wave polarity data by applying the FOCMEC package (Snoke 2003) as a low angle thrust (Fig. 14b), focal mechanism solutions could not be obtained for swarm activities collocated with slow earthquakes owing to their small magnitudes. Further study by adding the information of the S-phase waveforms may allow us to estimate a number of their focal mechanisms.
In Area 8, seismic activity swarms occurred twice with an interval of ~ 10 days (Fig. 13a). This finding suggests that monitoring of seismic swarm activity around the plate boundary may lay the foundation for estimating the state of interplate stick-slip changes with higher temporal resolution. Furthermore, episode 4 in Area 12 (with the largest number of interplate seismic activity but composed of smaller magnitude than that in episode 2) was not accompanied by any VLFE (Fig. 13b). Based on the relationship between the seismic energy release and magnitude (Gutenberg and Richter 1956), the summation of the seismic energy release for the swarm activity during episode 2 was estimated as 2.53 × 10 8 N m, much larger than that during episode 4 (6.77 × 10 7 N m), respectively. This, in turn, indicates that the total energy release during swarm activity was one of the indicators of the magnitude of SSEs.
Note that in most previous studies, the SSE areas and seismic swarm activity areas did not spatially overlap but were rather adjacent to each other. In this study, the occurrence of swarm earthquakes within the area of SSE and/or VLFE activity was also confirmed, but large uncertainties were identified in the location of shallow slow earthquakes. Hence, the detection of seismic swarm activity along the plate interface can strengthen the accuracy, thereby improving the estimation of slow-slip fault models. Furthermore, if the number of detections can be increased, as in the case of Nicaragua (Thorwart et al. 2013), one can retrieve physical parameters on the fault, such as diffusion coefficients and porosity based on the swarm seismicity catalog. As the minimum magnitude of the earthquake swarms obtained in this study was very small, compared to previous studies, it is essential to take advantage of the offshore seismic networks. The use of seafloor observation networks and hypocenter determination by using 3-D velocity model should extend the prospects for estimating spatiotemporal changes of b-values focused on the plate interface in the Nankai megathrust seismogenic zone, which has only previously been discussed with regard to spatial changes based on the JMA catalog without distinction between interplate and intraplate earthquakes (Nanjo and Yoshida 2018).

Conclusions
Simultaneous estimation of 3-D seismic wave velocity structure and hypocenter relocation in the Nankai Trough of the Kii Peninsula was performed using both onshore and offshore observation network data to investigate the spatiotemporal distribution of interplate earthquakes and their relationship with interplate coupling and slow earthquakes. First, the estimated velocity model revealed three large structural heterogeneities, indicated by previous studies: high-velocity plutonic rocks beneath Cape Shionomisaki, a subducted seamount off Muroto, and the subducted Paleo-Zenith Ridge. Moreover, our velocity model indicated the absence of other subducted seamounts, whose spatial scale was similar or larger than the subducted seamount off Muroto. Second, the hypocenter relocation unraveled the existence of patch-like interplate seismicity zones. Most of such zones were located within the seismogenic zone depth, outside both the strong interplate coupling areas and high (> 1.9) Vp/