Detection and Characterization of Earthquake Swarms in Nankai and Its Association With Slow Slip Events

Various types of slow earthquakes occur in the Nankai subduction zone, among which the slow slip events (SSEs) seem to connect the shallow and deep part of the subduction plate and further trigger seismic or aseismic transients in surrounding areas. This study uses a version of space‐time epidemic‐type aftershock sequence (ETAS) model with non‐stationary background rate to detect earthquake swarms that may be related to SSEs in the Nankai region. We detect 999 swarm sequences including 13,387 earthquakes, accounting for about 18 percent of all M ≥ 1.0 earthquakes in the selected region along the Nankai trough. The detected swarms are mostly distributed in the Hyuga‐nada, Bungo Channel, Kii Channel, and Tokai regions, forming a complementary pattern to SSEs. Nearly 70 percent of swarm events are located in western Nankai, with the southern part characterized by a higher occurrence rate, which is consistent with the smaller recurrence interval of SSEs in this region. By comparing occurrence times of swarm events to those of SSEs, we find that some SSEs are accompanied by enhanced swarm activities during their occurrence periods. In addition, the swarm rates increase to higher levels in the southern part of western Nankai during the 2003 and 2010 Bungo SSEs, indicating the high sensitivity of earthquake swarms to the migrations of slow earthquakes in this region.

These short-term activations in seismicity related to SSEs, usually referred to as earthquake swarms, are important to understand the aseismic slip process in subduction zones. Therefore, it is necessary to detect earthquake swarms objectively and analyze them in details to clarify the correlation between seismicity and aseismic transients. To extract earthquake swarms possibly related to aseismic slow slip episodes, we should eliminate the effect of earthquake-to-earthquake triggering, and measure the changes in seismicity caused by external mechanisms. Nishikawa and Ide (2017) proposed an objective method to detect earthquake swarms, in which they assumed that the occurrence rate of earthquakes in a region should be fitted by the Epidemic-Type Aftershock Sequence (ETAS) model, and the seismic anomalies could be treated as deviations from the model's expectation. Nishikawa and Ide (2017) used this method to detect and characterize earthquake swarms globally, and related the repeatedly occurrences of swarms to recurring SSEs in subduction zones. By the same method, Nishikawa et al. (2021) detected earthquake swarms along the Hikurangi trench and interpreted the pre-SSE swarm events by SSE-related fluid migration. Apart from the above method, earthquake swarms are also recognized as background seismicity anomalies. Marsan et al. (2013) introduced a time-varying background rate into the ETAS model (Ogata, 1998;Zhuang et al., 2002Zhuang et al., , 2005 and successfully detected seismicity anomalies related to aseismic loading in offshore central Japan. Peng et al. (2021) detected earthquake swarms as events with background rates higher than the 2-year averages in Taiwan and compared the swarm distributions with repeating earthquakes.
Although earthquake swarms associated with aseismic slips have been reported in many places, the spatiotemporal distributions of swarm events along the Nankai subduction zone remain unclear. Since the Nankai region contains various types of slow earthquakes that hold aseismic/seismic slips, it is important to detect and analyze the swarm activity in Nankai to understand the relationship between earthquake swarms and aseismic processes. This study aims at detecting earthquake swarms which may be related to aseismic transients in the Nankai region by following the methods in Nishikawa and Ide (2017) and Peng et al. (2021). Hence, we take the results obtained by both methods as earthquake swarms, and further merge them together to make detailed analyses. The spacetime ETAS model used in this study allows the background rate to fluctuate in time (Jia et al., 2020), and thus we can obtain the temporal migration patterns of background seismicity for the target region.
We organize the rest of this paper as follows: In Section 2, we analyze the catalog data and mention the Nankai SSE data used for comparison with earthquake swarms; In Section 3, we introduce the methodology and swarm detection criterions used in this study; Section 4 analyzes the detected swarms in details and discusses the spatiotemporal distributions of swarm events by comparing them with interplate coupling, SSEs, and migrations of slow earthquakes in the Nankai region; In the last two sections we make further discussions and summarize the conclusions.

Data
In this study, the JMA catalog data from January 2000 to February 2020 in Nankai subduction region is used for swarm detection. A manually defined boundary shown by the dash-black lines in Figure 1a is used in model fitting. To exclude most of the inland shallow earthquakes, we roughly select those events with focal depths close to the plate interface from the longitude-depth and latitude-depth plots of all earthquakes in the target region (Figures S1 and S2 in Supporting Information S1). The magnitude threshold is set as M1.0 to include as many earthquakes as possible, however, it is obvious that some small aftershocks are missing immediately following large earthquakes ( Figure S3 in Supporting Information S1). Here we refer to the method proposed in Helmstetter et al. (2006) to account for transient changes in M c following large mainshocks. From the four mainshocks with significant aftershock sequences in Figure S3 in Supporting Information S1, the M c (t, M) is estimated as M − 5.5 − 0.75 log 10 (t), where M denotes the mainshock magnitude. We use earthquakes with magnitudes The magnitude-frequency plots of earthquakes in different depth ranges. The b value is estimated by the MLE method (Aki, 1965), and the uncertainty is the standard error calculated through b values obtained by setting different magnitude thresholds from M1.0 to M3.5 with a magnitude bin of 0.1.

10.1029/2022JB025984
4 of 19 larger than M1.0 in data fitting and seismicity estimation and the magnitude threshold following M5 mainshocks or larger could be higher than M1.0 in the short term. Figure 1b plots the magnitude frequency distributions of earthquakes at different depths in the region for model fitting. Note that the two M7.0+ earthquakes and their aftershocks are not included in the b-value estimation. We can see that within 25-50 km, the number of M < 2.5 earthquakes are slightly lower than that expected from extending the G-R curve for M ≥ 2.5 events. This is because over 70% of M5.0+ earthquakes (33 out of 45) are located in this depth range, and the lower number of smaller aftershocks could be partially due to aftershock incompleteness immediately following the mainshock. For deep earthquakes (≥50 km), the numbers of M ≤ 3.5 events are lower than that expected from larger events, as the magenta and cyan lines indicate in Figure 1b.
Except for the earthquake catalog data fitted to the ETAS model, we also take the catalogs of tremors and repeating earthquakes (Uchida et al., 2020) as references for discussion. In addition, to help understand the earthquake swarms detected by this study, we refer to the information of SSE sources inverted from geodetic data in the Nankai region, including the long-term SSEs (Kobayashi, 2017;Kobayashi & Tsuyuki, 2019;Miyazaki et al., 2006;Ozawa et al., 2016;Takagi et al., 2019) and short-term SSEs (Okada et al., 2022).

Space-Time ETAS Model With Non-Stationary Background Rate
This study uses the space-time ETAS model proposed by Jia et al. (2020) to detect earthquake swarms in the Nankai region. The seismicity intensity function is given by where μ(t, x, y) denotes the background rate at time t and location (x, y). ( ) = ( − ) represents aftershock productivity of event i, in which M i and M c denote the magnitude of event i and the cut-off magnitude, respectively. ( ) = − 1 1 + − is the time probability density function (PDF) of aftershocks, where t denotes the elapsed time after the mainshock. The time PDF is normalized from the Omori-Utsu law (Utsu, 1961). For the spatial component, f(x, y; M i ) takes the following form where (x i , y i ) represents the epicenter location of event i. In the ETAS model, parameters (A, α, p, c, D 2 , q, γ) can be estimated by the maximum likelihood estimation (MLE) procedure (Zhuang et al., 2002).
As shown by Equation 1, the background rate μ is non-stationary, which is different from the conventional versions of the ETAS model (Console & Murru, 2001;Ogata, 1998;Zhuang et al., 2002Zhuang et al., , 2004. Such a change in background rate is useful for detecting seismicity anomalies. Hainzl and Ogata (2005) first applied the ETAS model to a fluid-induced earthquake swarm and quantified the effect of pore pressure changes in triggering earthquakes. For a region that hosts earthquake swarms possibly related to aseismic transients, the seismicity can fluctuate frequently due to the occurrences of aseismic slips. In this case, the conventional ETAS model is no longer suitable for fitting the seismicity because the background rate obtained by the model may be oversmoothed. Therefore, assuming a time-varying background rate enhances the sensitivity of ETAS model in detecting background rate changes.
In this study, we use the 1D and 2D Gaussian kernel functions in time and space, respectively, to estimate the background rate, that is, where h si is the epicentral bandwidth, h ti is the temporal bandwidth, and φ i denotes the probability of event i being a background event. h si is defined as where event j is defined as the event with transformed distance ranking n p th counted from minimum to maximum distance values between event i and other events, n p is an empirical parameter which usually takes values of 3-6 (Zhuang, 2011), ɛ is the minimum bandwidth which is usually set in the range of 0.01°-0.03°, equivalent to the location error of earthquakes, ω is a scaling parameter that transforms time differences between events into distances, and its dimension is distance divided by time. It is worth noting that estimating seismicity rates through the kernel method was also proposed by Choi and Hall (1999), and Zhuang et al. (2002) introduced the above variable bandwidths of the Gaussian kernel in the estimation of background seismicity rates.
Similar to the spatial bandwidth, h ti is defined as The use of the above temporal bandwidth is based on the following two reasons: (a) The bandwidth should be close to the temporal resolution of the data, which is basically at the scale of time difference to the closest event.
(b) Compared to the time lag to the closest event in time, the above selection is more stable. Jia et al. (2020) gave an estimation for the ω value in Equation 4, which equals to √ ∕ , where L, M, and T denote the longitude range, the latitude range, and the period of the fitting data, respectively. However, the fitting region and period are usually subjectively determined and hence an arbitrarily chosen ω value may not be proper for fitting the non-stationary background seismicity. Note that the minimum epicentral bandwidth ɛ in this study is set as 0.02°. Considering a good choice of minimum temporal bandwidth should not be so big that the changes in background seismicity would be oversmoothed, or so small that most aftershocks would be taken as fluctuations in background seismicity. We finally set this value as 20 days to satisfy the above condition and ω equals to 1/3 × 10 −3 .
In Equation 3, the background probability φ i can be calculated by the following equation: The estimation of μ, φ, and other model parameters can be solved by the iterative algorithm in Zhuang et al. (2002) (see also in Guo et al. (2017)).

Criterions for Swarm Detection
In this section, we introduce two methods used for swarm detection in Nishikawa and Ide (2017) and Peng et al. (2021). The first method takes strong earthquake clusterings deviated from the expectations of the ETAS model as earthquake swarms and thus relatively short-term seismic transients can be detected. In contrast, the second method can identify earthquake swarms by extracting temporal variations in background seismicity. Thus, the combination of both methods can yield a broad spectrum of temporal transients reflected by anomalous seismicity.
For the first method, we follow the swarm detection method proposed by Nishikawa and Ide (2017). The number of events in a region S within period (0, T ] expected from the ETAS model can be calculated by Using the above equation we can transform the occurrence time t i of event i in the region S into τ i ≡ Λ(t i ) (Ogata, 1988). The transformed occurrence time τ i should follow a standard Poisson process if the ETAS model fits the seismicity in S, as demonstrated by Nishikawa and Ide (2017), and the plot of cumulative event number against their transformed occurrence times is approximately linear with the unit slope. In fact, the slope values in a region derived from the ETAS model will not be constant and could often deviate from unity. Slope values greater than one suggest that the occurrence frequency is higher than the model expectation, and vice versa.
Here we introduce the "1σ error criterion" used by Nishikawa and Ide (2017). That is, an earthquake swarm should satisfy the condition τ i+1 − τ i + σ < 1 for five successive events or more, where τ i+1 − τ i denotes the time difference between two successive events, and σ is the standard deviation of the event number between τ i and τ i+1 , equaling to √ +1 − . By solving the above inequality we can get its solution: τ i+1 − τ i < 0.382. It is worth pointing out that the space window used for the swarm detection under this criterion is 0.2° × 0.2°, which means that the study region is divided into 0.2° × 0.2° grids and the spatial term in Equation 7 is intergrated over each grid. Therefore, swarms are defined as earthquake sequences which should contain at least five successive events with transformed time differences less than 0.382 in each grid.
Since the background seismicity varies with time in the non-stationary ETAS model, it is possible that earthquake swarms are fitted by temporal variations of background seismicity and thus cannot be successfully detected by the above criterion. Therefore, we need to introduce the second criterion for extracting anomalies with high background rates to complete the whole swarm detection process. Peng et al. (2021) proposed a method to detect earthquake swarms from non-stationary background seismicity rates in Taiwan. Here we use the same method for the Nankai region. First, we divide the whole study region into 0.2° × 0.2° grids, and calculate the background rate μ(t) at each earthquake occurrence time for those grids including over 50 events (over 93% of all earthquakes are selected) by using the 1D Gaussian kernel in Equation 3. Second, the 2-year moving averages of background seismicity are estimated for the selected grids, taken as μ ave (t). As the final step, a background index is defined as and we select earthquakes with background indexes larger than a specific δ value as potential swarm events. Furthermore, to group the selected events, we assume that earthquakes separated by less than 20 km and 10 days belong to the same cluster.

Data Analysis
For data analysis, we fit the nonstationary background space-time ETAS model to the JMA data, select the earthquake swarms by using the two methods outlined in Section 3.2, and discuss the spatiotemporal distribution of those swarms and their correlation to interplate coupling and SSEs.

Fitting Results
The model parameters fitted by the MLE procedure are A = 0.11 events, c = 0.0018 days, α = 1.17, p = 1.05, D 2 = 5.92 × 10 −6 deg 2 , q = 1.69, and γ = 0.61. This set of parameters is used to calculate the transformed times of earthquakes Λ(t) in Equation 7 for a given region. It is worth pointing out that the fitting results of the ETAS model only represent the mean level of aftershock triggering in the region. Hence, mainshock-aftershock sequences may be identified as anomalous swarms by this approach if the productivity of such sequences is larger than the model's expectation ( − ) .
About 84% of all earthquakes are background events according to the background probabilities φ (Equation 6) estimated by the ETAS model. Such a high background ratio is also indicated by the parameter pair A and α. The value of A means that an M1.0 mainshock only triggers 0.11 aftershocks on average, and an M5.0 mainshock directly triggers about 0.11 × e 1.17*(5−1) = 11.9 aftershocks. The low aftershock productivity of large earthquakes could be attributed to the point source assumption in the ETAS model, as Hainzl et al. (2008) demonstrated that the space-time ETAS model yielded a smaller α value and a larger A value compared to the results obtained by the temporal ETAS model. Guo et al. (2017) further verified this and found that the lower estimation of α value could be corrected by incorporating rupture geometries of large mainshocks in the space-time ETAS model. Also, the low estimates of aftershock productivity for large earthquakes might be associated with the interplate coupling, as Hainzl et al. (2019) found that the low aftershock productivity along depth in the subduction zone of Northern Chile could be interpreted by the low coupling ratio in deeper parts of the plate interface. For the Nankai region, stress changes induced by large mainshocks in low coupling areas could be partially released by adjacent aseismic processes, and thus less aftershocks are generated by the mainshocks. Therefore, the low estimate of aftershock productivity could be a common phenomenon in subduction zones. As the result, the detected earthquake swarms in this study may include some aftershocks following large mainshocks where the aftershock rate increases to a much higher level than the model's expectation. More analyses in the categorization of swarm types based on earthquake magnitudes can be referred to Section 4.3.

Detected Swarms
Since this study uses two methods for swarm detection, it is necessary to make a comparison between the results obtained from both methods. In the following we explain more details in swarm detection by the second method. Figure 3 gives an example to illustrate how the background index in Equation 8 is calculated. Note that the minimum temporal bandwidth in Equation 5 for estimating background seismicity is set as 20 days, and hence the deviations of rate values from the 2-year moving averages suggest temporal changes in background seismicity. As presented by Figure 3b, those black circles above the average line mean that the background seismicity rate is at a relatively high level at the occurrence times of these events. To determine an appropriate threshold δ value, we refer to the rule in Peng et al. (2021): the extracted swarm catalog should exclude most of aftershock sequences which generally have larger mainshock magnitudes. The number of swarm clusters and the number of M ≥ 5.0 clusters decrease with increasing δ values (shown in Figure 4). We finally choose δ = 20 as the threshold value to retain more swarm clusters and manually remove 19 clusters including 20 M ≥ 5.0 earthquakes in the results.
To illustrate the differences of results obtained from the two detection methods, we show the results of node (131.6° − 131.8°E, 31.2° − 31.4°N) from the first method (Nishikawa & Ide, 2017) and the second method (Peng et al., 2021) in Figures 2 and 3, respectively. In Figure 2, there were 107 events belonging to 17 swarm sequences detected by the first method. The largest sequence containing 17 earthquakes occurred in December 2007, and it lasted for about 3 hr in a spatial range of 0.06° (longitude), 0.04° (latitude), and 4 km (depth). The maximum magnitude of this sequence was M3.5. In contrast, 288 of all 1147 events in the node were detected as potential swarm events by the second method, as shown by the red circles in Figure 3b. The largest sequence had 47 events including an M4.8 mainshock-aftershock sequence. This sequence had a duration of 30 days within a range of 0.12° × 0.16° × 19 km, and it completely included the largest sequence from the first method. In Figure 3b, some swarm events detected by the first method were characterized by high background index values and recognized as swarm anomalies by the second method as well. Further analyses suggested that there were 53 events detected as swarm anomalies by both methods in the selected node. Additionally, we checked all the swarm sequences and found that substantial differences existed between results obtained the two methods: The swarm sequences from the second method were characterized by larger event numbers (13.6 events), wider spatial extensions (0.22° × 0.21° × 27 km), and longer durations (27.5 days) on average, whereas swarms from the first method merely had 7.7 events, 2.3 days of durations, and spatial extensions of 0.06° × 0.06° × 9 km on average. Therefore, the second method tends to preserve more complete sequences but extends them longer in time and wider in space than the first method.

Spatiotemporal Distribution of Swarms
Following the criterion proposed by Nishikawa and Ide (2017), we initially detect 570 swarm sequences including 4,370 events within the fitting region. In contrast, adopting the method in Peng et al. (2021) yields 865 anomalous sequences including 11,727 earthquakes. To maintain the consistency of event numbers of the sequences detected by both methods, we require that each sequence has at least 5 events. Since both methods adopted in swarm detection are independent from each other, their results are overlapped to some extent. Therefore, to make a consistent   swarm catalog, we merge the results obtained from both methods and regroup all events with the aforementioned rule: earthquakes separated by less than 20 km and 10 days belong to the same sequence. As the final result, the recompiled swarm catalog contains 999 sequences including 13,387 events.
To further categorize the detected swarm sequences, we divide all sequences into two types: the swarm-like sequence and the foreshock/aftershock-like sequence. This categorization is based on the magnitude difference between the largest event and the second largest event ΔM in a sequence. If ΔM is larger than 1.0, then the sequence belongs to the foreshock/aftershock-like type; otherwise, the sequence belongs to the swarm-like type. The value 1.0 is referred to the Båth law (Båth, 1965) which states that the average difference in magnitude between a mainshock and its largest aftershock is 1.2. The total 999 sequences are divided into 834 swarm-like sequences and 165 foreshock/aftershock-like sequences, including 11,363 and 2024 earthquakes, respectively. Note that this categorization is quite rough because some individual foreshock/aftershock sequences could have magnitude-difference values smaller than 1.0, and thus some foreshock/aftershock-like sequences could be identified as swarm-like sequences. However, using a smaller difference value in categorization may lead to an overestimation of the number of foreshock/aftershock-like sequences as well. A difference value of 0.6 results in 683 swarm-like sequences of 9,450 events and 316 foreshock/aftershock-like sequences of 3,937 events. Compared to the results obtained for difference value of 1.0, the number of swarm events decreases to about 83%. Therefore, the difference value 1.0 indicates that we only separate those significant foreshock/aftershock-like sequences from the results. Figure 5 presents the spatial distribution of all swarm sequences. We can see that the swarm events are mostly distributed in off-Kyushu island, Bungo channel, western Shikoku island, Kii channel, southwest Kii peninsula, and Tokai regions. To investigate the spatio-temporal distribution of earthquake swarms, we select four regions whose boundaries are shown in magenta dash lines in Figure 5. Regions H, B, K, and T represent the Hyuga-nada region, Bungo channel, Kii channel, and Tokai region, respectively. The temporal distributions of swarms are presented in Figure 6, from which we can see that these four regions are characterized by different occurrence rates of earthquake swarms. In western Nankai, 569 swarm sequences including 481 swarm-like and 88 foreshock/aftershock-like sequences are detected in Regions H and B, making up more than 56% of all sequences (Table 1). Specifically, there are 320 and 249 sequences in Regions H and B, respectively. While for the Kii region, 190 swarm-like sequences consisting of 1944 events and 39 foreshock/aftershock-like sequences consisting of 407 events are detected as earthquake swarms. The occurrence rate in Tokai becomes even lower as the numbers in Table 1 indicate. More than 90% of all swarm events occur at depths shallower than 50 km, and only 571 swarm-like events and 168 foreshock/aftershock-like events occur at 50+ km depths, as we can see from the distributions of red filled circles in Figure 5.
From Figure 6 we can see that earthquake swarms accumulate at an almost constant rate in each region for the long term, which can be inferred from the slopes of red curves in the subfigures. However, perturbations in short times exist when a long swarm sequence occurs, reflected by the steep increases of cumulative numbers of swarm events. Among all the swarm sequences, we present four swarm-like sequences with each containing over 100 events in Figures S4 and S5 in Supporting Information S1. The magnitude-time plots of the selected four sequences, ordered by occurrence times, are shown by Figure S4 in Supporting Information S1. The four sequences started on 2001/03/28, 2003/07/26, 2009/11/11, and 2019/9/13, and lasted for 28, 7, 110, and 171 days, respectively. The first and second sequences were more concentrated in space, whereas the third and fourth ones Figure 5. The spatial distribution of detected anomalous earthquake sequences: (a) swarms; (b) aftershock sequences. The colored circles denote the focal depths of swarms and aftershock sequences, respectively, with sizes representing their magnitudes. The light blue filled circles mark the locations of repeating earthquakes from Uchida et al. (2020). The pink ellipses mark the regions that host long-term slow slip events (SSEs) in the Nankai region (Obara, 2020). The gray dots represent the epicenters of LFEs. The polygons in black and red dash lines indicate the low and high interplate coupling areas in Hyuga-nada region inverted by Kimura et al. (2019). The depth contours of subducting plate interface shown by gay curves are referred to F. Hirose et al. (2008), Nakajima and Hasegawa (2007), and Baba et al. (2002).
were very scattered and extended to more than 60 km in the horizontal components ( Figure S5 in Supporting Information S1).
Although the model used in this study only incorporates the focal epicenters and thus has no resolution in depth, we plot the occurrence depths of swarm events with respect to the depths of subducting plate along seven profiles in Figure 7 to investigate the depth distribution of swarms. In Region H, swarm events are mostly located in areas between the 10 and 30 km contour lines of the subducting plate Figure 5, however, from Figures 7a and 7b we can see that earthquake swarms extend from 10 to 50 km in depth within this region, forming a very thick seismogenic layer. In Region B, the focal depths of swarm events are mainly presented by Figure 7d, showing one significant feature of earthquake swarms in Bungo channel that swarm events become more scattered as depth goes deeper. Relatively consistent seismogenic layers exist along Profiles EE' and FF' in Regions K and T, with swarm depths extending to about 60 and 50 km, respectively. Along Profile GG' some swarm events occur in depths shallower than the plate interface (Figure 7g), which could include inland earthquakes that are not completely removed from the catalog.

Earthquake Swarms and Interplate Coupling
In Figure 5 (a) and (b), the swarm-like and foreshock/aftershock-like events are overlapped to a large extent. However, a significant difference exists in the along-strike direction off the Kyushu island. In the southern part (Region H), 283 swarm-like and 37 foreshock/aftershock-like sequences occur, while in the northern part (Region B), the number of swarm-like sequences decreases to 198 whereas the number of foreshock/aftershock-like sequences increases to 51. Such difference could result from the coupling heterogeneity of the subducting plate in the off-shore Kyushu region. As the coupling  ratio distribution along the Nankai trough inverted by Kimura et al. (2019) shows that the southern off-Kyushu region is characterized by a relatively low coupling ratio whereas most areas in the northern off-Kyushu region are highly coupled.
We select two regions respectively characterized by low (R1) and high (R2) coupling ratios whose boundaries are indicated by the black and red dash lines in Figure 5. The boundaries of two polygons are taken from the results of coupling ratio distribution in Kimura et al. (2019). 2910 swarm-like events and 243 foreshock/aftershock-like events are located in R1. In contrast, the event numbers of two swarm types in R2 are 1666 and 420, respectively. Given consideration to the area of R1 is much larger than that of R2, the difference in distributions of swarm types becomes larger in these two selected regions. A possible interpretation for this phenomenon could be as follows: The weak interplate coupling allows for a high rate of stress loading in the downdip direction and thus leads to frequent occurrences of aseismic transients such as SSEs, who cause stress perturbations in surrounding areas and further alter local seismicity, partially reflected by anomalous swarm sequences.
Following Hainzl et al. (2019), which investigated the seismicity in northern Chile subduction zone and found a linear relationship between aftershock productivity and interplate coupling in depth, we reconstructed aftershock productivity at different depths and in the two regions with different coupling ratios to investigate the aftershock productivity in Nankai. Results from Figure 8a show that aftershock productivity significantly decreases as the depth becomes larger, which could be partially attributed to the low interplate coupling in deep parts. In Figure 8b, the aftershock productivity of M ≤ 3.0 earthquakes is relatively high in low coupling areas, while the aftershock productivity of M ≥ 4.0 earthquakes in high coupling areas is relatively high. This indicates that aftershock productivity positively correlates with the coupling ratio along subducting depths also holds true for the Nankai region. However, the along-strike aftershock productivity of small earthquakes seem to contradict with the effect of interplate coupling. Note that quite a proportion of earthquakes belong to intraplate events as we can see from Figure 7, and it is difficult to rule out these events by selecting events closer to the plate interface due to the location errors. Moreover, cutting off events farther from the plate interface leads to the underestimate of aftershock productivity since the ETAS model in this study has no resolution in depth. Therefore, the interplate  Figure 5. The x axis denotes the distance measured from the northern end of the profile. Events with horizontal distances less than 30 km from each projecting profile are plotted. Black dots mark the depths of all earthquakes in the selected catalog. Red and blue dots denote events from the swarm-like and foreshock/aftershock-like sequences, respectively. The black solid line in each plot represents the upper interface of subducting plate, which is referred to F. Hirose et al. (2008), Nakajima and Hasegawa (2007), and Baba et al. (2002).
coupling is just one of the factors which seems to control aftershock productivity along the subducting plate.

Relationship Between Earthquake Swarms and SSEs
To inspect the relationship between SSEs and earthquake swarms in the Nankai region, we plot the spatial distributions of swarm numbers and SSEs in Figure 9. The source parameters of the long-term and short-term SSEs are referred to Takagi et al. (2019) and Okada et al. (2022), respectively. From Figure 9 we can see that earthquake swarms generally complement with SSEs in Hyuga-nada, southwest Kii peninsula, and Tokai regions, whereas areas of high swarm rates overlap with SSE ruptures in Bungo channel and west Shikoku island regions. Swarm events in Hyuga-nada mainly distribute in the updip direction of SSEs, although their depths extend to 50 km (Figures 7a-7c). While in western Shikoku and southwestern Kii, swarm events are detected at more than 50 km depth, downdip the SSEs. To quantify the spatial correlations between swarm events and SSEs, we count the event numbers within the rupture sources of SSEs and find that about 52% of swarm events (6,907 out of 13,387) occur within the long-term SSE ruptures. In contrast, about 42% of all events (31,860 out of 76,621) are located in the source areas of long-term SSEs. For shortterm SSEs, there are respectively 84% of swarm events (11,215) and 80% of all events (61,343) within the short-term SSE ruptures. This suggests that swarm events are overlapped with SSEs to a larger extent than all the events.
In the western Nankai region, earthquake swarms, SSEs, and tremors are widely overlapped as shown by Figures 5 and 9. Hence, to better understand the temporal distributions of these three types of events in western Nankai, we plot the time-latitude distributions of swarm events, short-term and long-term SSEs, repeating earthquakes, and tectonic tremors in the background seismicity patterns in Figure 10. Note that the swarm events in Figure 10 are mostly located in areas of high background seismicity rates because the second criterion used for detecting swarms requires that the background rate at the occurrence time of one event is higher than its 2-year moving averaged background rate. Also, we find that repeating earthquakes in the off-Kyushu region (between 31°N and 33°N) from Uchida et al. (2020) are mainly distributed in areas of high background rate. Specifically, 94 of all 525 repeating earthquakes are included in the detected swarm sequences, mostly distributed off the Kyushu island.
To compare the spatio-temporal distribution of earthquake swarms to that of SSEs, we define a swarm sequence that is correlated with SSEs should satisfy the following two conditions: (a) it should have one or more events with epicentral distances less than 10 km from the fault boundaries of a SSE; (b) the events in (a) should occur within the time interval of (Ts − 30, Ts + dT + 30), where Ts and dT in days respectively represent the starting time and duration of the SSE. By this standard, there are 794 and 908 earthquakes, belonging to 135 and 107 swarm sequences, correlated with short-term and long-term SSEs, respectively. Note that the comparison between earthquake swarms and long-term SSEs in Takagi et al. (2019) is restricted to the western Nankai region.
To further quantify the spatiotemporal correlations between earthquake swarms and SSEs, we divide the whole study region into small space-time volumes (0.2° × 0.2° × 10 days), and we count the numbers of swarm events inside (N in ) and outside (N out ) of the volumes covered by the SSEs. For comparison, N in and N out are normalized by the numbers of volumes covered (A in ) or not covered (A out ) by the SSEs. For long-term SSEs, the ratio N in /A in is about 1.8, whereas N out /A out equals to 2.1. A similar result is obtained for the short-term SSEs: N in /A in (2.2) is smaller than N out /A out (2.7). These results suggest that earthquake swarms detected in this study are not spatiotemporally correlated with SSEs in general. In the case we consider the aforementioned wider definition: we count swarm events with occurrence times ±30 days with respect to SSEs which are distributed within 10 km from the source boundaries of SSEs as inside events, N in /A in is still smaller than N out /A out (1.9 vs. 2.2 for long-term SSEs, 2.5 vs. 2.7 for short-term SSEs). However, for some individual SSEs, N in /A in is larger than N out /A out (Figures S6 and S8 in Supporting Information S1). 4 out of 19 long-term SSEs and 32 out of 257 short-term SSEs are characterized by high N in /A in values, indicating that only a small proportion of SSEs are accompanied by enhanced swarm activities.
Relatively rare long-term SSEs are observed near the Kii channel and Tokai regions. Kobayashi (2017) reported two long-term SSEs in the Kii Channel during 2000-2002 and 2014-2016, with moment magnitudes of M6.6 and M6.7, respectively. During these two SSEs, 303 and 366 swarm events belonging to 28 and 39 sequences occurred within Region K. For the Tokai region, two long-term SSEs between mid-2000 and mid-2005, and between 2013 and 2015 were reported (Miyazaki et al., 2006;Ozawa et al., 2016), with moment magnitudes of over M7.0 and M6.6. The spatial range of the two Tokai SSEs were almost limited in 34.5°-35.5°N, 137°−138°E, in which we respectively detected 141 and 180 swarm events belonging to 17 and 10 sequences during these two SSEs. In addition, Kobayashi and Tsuyuki (2019) detected a Mw6.4 long-term SSE in the eastern-most part of the Kii Peninsula from April 2017 to October 2018, with slips mainly distributed in the range of 34°-35°N, 136°−137°E. 27 swarm events from 4 sequences occurred during this SSE event. The N in /A in and N out /A out ratios for the above five long-term SSEs are 1.7 and 1.8, respectively. The N in /A in values of first Tokai SSE (1.85) and the east-Kii SSE (2.4) are relatively higher, suggesting stronger correlations between earthquake swarms and these two SSEs ( Figure S7 in Supporting Information S1). Uchida et al. (2020) carried out a comprehensive analysis of the seismic and geodetic data and revealed the existences of deep slow and shallow fast migration patterns in western Nankai, manifested by deep long-term SSEs, and shallow VLFEs, repeating earthquakes and tremors, respectively. The two slow migration processes were associated with the 2003 and 2010 long-term SSEs in Bungo channel. In Figure 10, the slow migrations of long-term  Hirose et al., 2010;Uchida et al., 2020), are indicated by the thick green lines in Figure 10.

Migrations of Slow Earthquakes and Accompanied Swarms in Western Nankai
From Figure 10 we can see that earthquake swarms distribute extensively in western Nankai, and three belts with high background rates exist in the western Nankai region: Near 31.2°N, near 32.1°N, and around 32.6°N. However, no significant migrations of earthquake swarms were found, neither all long-term SSEs were accompanied by enhanced swarm activities during the migrations. For discussion, we roughly divided the region into three parts: The southern Part 31° ∼ 32°N, the middle Part 32° ∼ 33°N, and the northern Part 33° ∼ 34°N. The swarm rate in the southern part was 0.47 event/day during the whole study period, and it had increased to 0.62 events/day from 2002 to 2005 and 0.53 events/day from 2007 to 2012. In contrast, the swarm rate in the middle part was 0.50 events/day, and decreased to 0.39 events/day from 2002 to 2005, and then increased to 0.50 events/ day from 2007 to 2012. For the northern part, no significant variations in the swarm rate were observed during the above three periods, and the mean rate was 0.24 events/day with a standard error of 0.005 events/day. Moreover, we investigated the swarm rates during two short periods (from June 2003 to November 2003 and from November 2009 to May 2010) when the two fast migrations occurred according to Uchida et al. (2020). In the southern part, the swarm rate increased to 1.0 events/day during the first period and 0.83 events/day during the second period, much higher than the averaged 0.47 events/day; In the middle part, the swarm rate decreased to 0.42 events/day during the first period and then increased to 0.70 events/day during the second period; In the northern part, the swarm rate first decreased to 0.21 events/day and second increased to 0.29 events/day.
In Figures 11a and 11b, the spatiotemporal distributions of earthquake swarms during two long periods (from August 2001to May 2004and from November 2006to February 2011 are presented, and the two slow migrations indicated by long-term SSEs in Takagi et al. (2019) are plotted. Although the SSEs showed clear migrations from offshore Kyushu in the south to west Shikoku in the north, earthquake swarms extended widely from the south to the north in the whole region and showed no significant temporal migrations. Note that the three SSEs in Figure 11b characterized by high N in /A in ratios (rectangles in solid red lines in Figure 10) were mainly located in the south part around 32°N,namely SSE No. 1,No. 3,and No. 6, indicating a stronger correlation between earthquake swarms and SSEs in this region. Two fast migrations indicated by the occurrences of shallow VLFEs, repeating earthquakes, and tremors occurred at the beginnings of 2003 and 2010 Bungo channel SSEs (Uchida et al., 2020). To analyze the swarm sequences that The rectangles in solid and dashed lines represent long-term and short-term SSEs inverted by Takagi et al. (2019) and Okada et al. (2022), respectively. The rectangles in blue and red respectively denote SSEs with larger N out /A out and N in /A in values. The inverted cyan triangles mark repeating earthquakes located by Uchida et al. (2020). The thick green lines mark the occurrence times of short-term migrations in 2003 and 2010. are possibly related to these two migrations in details, we present the epicenter-time distributions of swarm events in the above two short periods in Figures 11c and 11d, together with repeating earthquakes, long-term SSEs, and tremors. The depth and magnitude distributions of these swarm events versus time in the southern, middle, and northern parts are shown by Figures S9-S14 in Supporting Information S1. In Figure 11c, 155 events existed in the south of 32°N, among which 113 events occurred within 5 days after 2003/7/26, spatially concentrated in a 0.08° × 0.06° range centered on (131.95°E, 31.26°N). The depth range of these 155 events was between 12 and 44 km, and the maximum magnitude was M4.4 ( Figure S9 in Supporting Information S1). In the middle part between 32°N and 33°N, 63 events occurred from July 11th to November 6th in 2003, among which 35 events occurred between August 18th and September 19th. The swarm events in this part had a largest magnitude of M2.5 and extended from 6 to 37 km in depth ( Figure S10 in Supporting Information S1). In the northern part, only 31 swarm events occurred in the north of 33°N, which were located deeper than 40 km. These events were scattered in space and time (Figure 11c and Figure S11 in Supporting Information S1), and might not be associated with the fast migration process. However, a large number of tectonic tremors started in the end of August 2003, downdip the 2003 Bungo SSE, suggesting that tremors were more sensitive to SSE-related migrations because of their smaller magnitudes.
In Figure 11d, There were 161 events in the south of 32°N, and 66 of them were located in a 0.16° × 0.06° range centered on (131.95°E, 31.29°N). This sequence began on 2010/2/4 and lasted for about 4.0 days, during which period 5 repeating earthquakes occurred. The largest event in this sequence was an M4.3 earthquake and the highly concentrated events extended from 17 to 31 km in depth, about half of the cluster extension ( Figure  S12 in Supporting Information S1). 137 swarm events were located in the middle part between 32°N and 33°N in Figure 11d, which extended widely in space: more than 60 km horizontally and more than 35 km in depth. No highly concentrated sequences were found in this region and the largest magnitude of these events was M3.5 Figure 11. Spatio-temporal distributions of earthquake swarms in the corresponding periods, with colors denoting the occurrence times. In (a and b), the rectangles in black dot lines represent the long-term slow slip events (SSEs) detected in Takagi et al. (2019), and the thick black arrows indicate the slow migration of SSEs. In (c and d), the rectangles in black dot lines represent the 2003 and 2010 Bungo channel SSEs, respectively. The thin black arrows indicate the fast migration of shallow very-low-frequency earthquakes, repeating earthquakes and tremors demonstrated by Uchida et al. (2020). The inverted triangles denote repeating earthquakes in the corresponding periods.
( Figure S13 in Supporting Information S1). In the north of 33°N, 56 events deeper than 37 km extended widely in space and time with no densely distributed clusters observed, whereas some sequences formed by tectonic tremors were activated on 2010/2/20 downdip the 2010 Bungo SSE ( Figure S14 in Supporting Information S1).

Possible Triggering Mechanisms for Earthquake Swarms
In the Nankai region, SSEs are not accompanied by activated earthquake swarms for most of the cases, however, some SSEs indeed have enhanced the swarm rate during their occurrences. Stronger spatiotemporal correlations between earthquake swarms and SSEs have been reported in other subduction zones. For example, the seismic transients detected by Marsan et al. (2013) in central Japan mostly occurred in decoupled zone with a recurrence interval of 5.9 years, which could be related to the Boso SSEs (H. Hirose et al., 2012). Fukuda (2018) demonstrated that the occurrence rates of earthquake swarms during the Boso SSEs were positively correlated with SSE slip rates, suggesting that SSE-induced stress loading could be the triggering mechanism for swarms.
For the Nankai region, we detected 355 swarm events belonging to 78 sequences and 74 swarm events belonging to 12 sequences within 30 days before the short-term and long-term SSEs, respectively. Although there is a lack of focal mechanism data, swarm depth distributions in Figure 7 suggest that earthquake swarms detected in this study mostly belong to intraplate events, which are consistent with earthquake swarms along Hikurangi detected by Nishikawa et al. (2021). To interpret the pre-SSE swarms in Hikurangi, Nishikawa et al. (2021) proposed that SSE-related fluid migrations could trigger these pre-SSE events. As Warren-Smith et al. (2019) found that there could be an accumulation and release process of fluid pressure within the subducting slab before and during SSEs according to the changes of stress ratio retrieved from earthquake focal mechanisms in the Hikurangi subduction zone. Therefore, these pre-SSE swarms in Nankai could also indicate that fluid migration occurs before SSEs and may enhance nearby seismicity rates in the short term. Furthermore, fluid drainage from megathrust to overlying plate may be induced by SSEs, possibly leading to the occurrences of co-and post-SSE earthquake swarms (Nakajima & Uchida, 2018;Warren-Smith et al., 2019). Hence, those intraplate swarm events occurred during and after the SSEs in Nankai could be partially explained by the fluid drainage process during SSEs.
However, the SSE-related triggering mechanisms can explain the occurrences of only a small proportion of earthquake swarms. According to the N in /A in and N out /A out ratios, many swarm events are located outside of the space-time windows covered by the SSEs. This indicates that the earthquake swarms in this study reflect shortterm seismic transients by selecting deviations to the ETAS model and relatively long-term seismic transients by selecting background rates higher than the 2-year moving averages. These changes in seismicity include some mainshock-aftershock sequences with high aftershock rates because the model parameters only represent the mean level of aftershock triggering within the study region. Also, some swarms caused by undetected smaller-scale SSEs or other types of slow earthquakes could lead to larger N out /A out values. Consequently, the insignificant correlation between swarms and SSEs might be caused by the intrinsic limitation in the methods used for swarm detection and the incompleteness of the SSE catalogs.

More on Earthquake Swarms and Interplate Coupling
The occurrence rates of earthquake swarms seem to correlate with interplate coupling to some extent, as Figure 5 shows that there are more swarm events observed in the low-coupling region inverted by Kimura et al. (2019). This could be explained by the correlation between SSE recurrence intervals and updip interplate locking. Takagi et al. (2019) detected 24 long-term SSEs in western Nankai and found that more frequent SSEs occurred downdip the low interplate coupling areas. Since earthquake swarms can be activated by SSE-related fluid migrations and stress changes, we attribute the high occurrence rate of earthquake swarms in the low coupling region to the more frequent downdip SSE occurrences. However, more foreshock/aftershock-like swarm sequences exist in high coupling region (Figure 5b), which contradicts with the above interpretation. Given consideration to the fact that foreshock/aftershock-like sequences detected in this study are characterized by larger maximum magnitudes (Table 2), the frequent occurrences of this swarm type in high coupling regions indicate larger aftershock productivity of earthquakes. Hainzl et al. (2019) proved that aftershock productivity was linearly correlated with seismic coupling in the Northern Chile subduction zone, which could support our findings in Nankai. Nevertheless, stochastic reconstruction of aftershock productivity in Figure 8 shows a systematic decrease of aftershock productivity along depth but the larger aftershock productivity in high coupling region only exists in M ≥ 4.0 earthquakes. Consequently, the difference in aftershock productivity along the subduction strike due to interplate coupling is less significant than that along depth.

Correlation Between Earthquake Swarms and Repeating Earthquakes
326 of 525 repeating earthquakes in Uchida et al. (2020) are located in offshore Kyushu and west Shikoku, mostly covered by earthquake swarms detected by this study (Figure 5). A detailed comparison shows that 94 repeating earthquakes are included in swarm sequences, and 324 repeating earthquakes are located within 5 km from their closest swarm events. Such a strong spatial correlation between earthquake swarms and repeating earthquakes was also presented by Peng et al. (2021), which detected earthquake swarms in Taiwan and found that 86% of repeating earthquakes were located within 5 km from the closest swarm events in southeast Taiwan. The analysis of spatiotemporal distributions of SSEs, repeating earthquakes, tremors, and earthquake swarms in western Nankai (Figures 10 and 11) shows that the deep slow and shallow fast migration modes proposed by Uchida et al. (2020) are accompanied by seismic activities including not only repeating earthquakes but also earthquake swarms in some localized areas.

Conclusion
We use an non-stationary space-time ETAS model to detect anomalous seismic transients in Nankai subduction zone. By following the methods proposed by Nishikawa and Ide (2017) and Peng et al. (2021), 13,387 events belonging to 999 swarm sequences including 834 swarm-like and 165 foreshock/aftershock-like sequences are detected, accounting for 17.5% of all earthquakes in the target region. Almost 70% of swarm events (9,068 of 13,387) are located in western Nankai with longitudes less than 133°, suggesting a higher seismicity rate in this region. The earthquake swarms spatially complement with SSEs and tectonic tremors in a large scale, and partially overlap with SSEs in the Bungo Channel, southwest Kii Peninsula, and Tokai region. Detailed analyses in the temporal distributions of selected swarm sequences in western Nankai suggest that some individual SSEs are accompanied by enhanced swarm activities. Although the correlation between earthquake swarms and SSEs is not a common phenomenon in the Nankai region, detecting and characterizing anomalous seismic transients is important to understand the seismic and aseismic processes along the Nankai subduction zone.