Development of a detection method for short-term slow slip events using GNSS data and its application to the Nankai subduction zone

Using global navigation satellite system (GNSS) data to detect millimeter-order signals of short-term slow slip events (S-SSEs) and to estimate their source parameters, especially duration, is challenging because of low signal-to-noise ratio. Although the duration of S-SSEs in the Nankai subduction zone has been estimated using tiltmeters, its regional variation has never been quantitatively studied. We developed an S-SSE detection method to estimate both the fault model and duration with their errors based on the detection methods developed by previous studies and applied it to a 23-year period of GNSS data in the Nankai subduction zone. We extracted S-SSE signals by calculating correlation coefficients between the GNSS time series and a synthetic template representing the time evolution of an S-SSE and by computing the average of correlation coefficients weighted by the predicted S-SSE signals. We enhanced the signals for duration estimation by stacking GNSS time series weighted by displacements calculated from the estimated fault model. By applying the developed method, we detected 284 S-SSEs from 1997 to 2020 in the Nankai subduction zone from Tokai to Kyushu and discussed their regional characteristics. The results include some newly detected S-SSEs, including events accompanying very low-frequency earthquakes and repeating earthquakes in offshore Kyushu. Our study provides the first geodetic evidence for synchronization of S-SSEs and other seismic phenomena in offshore Kyushu. We estimated the cumulative slip and duration, and their error carefully. We also estimated the average slip rate by dividing the cumulative slip by the cumulative duration. This study clarified that the average slip rate in western Shikoku was approximately twice as that in eastern Shikoku and Kyushu. These regional differences were statistically significant at the 95% confidence interval. Multiple factors can influence the regional characteristics of S-SSEs, and we speculate that the subducting plate interface geometry is one of the dominant factors.


Introduction
Slow slip events (SSEs) are transient aseismic fault slip events occurring in a transition zone between seismogenic and stable creep zones (Obara and Kato 2016). SSEs have been observed in several plate boundaries along the Pacific Rim owing to the dense global navigation satellite system (GNSS) network since the 1990s (Hirose et al. 1999;Dragert et al. 2001;Lowry et al. 2001;Ozawa et al. 2003;Douglas et al. 2005;Ohta et al. 2006;Heki and Kataoka 2008;Ruiz et al. 2014;Rousset et al. 2019). SSEs are occasionally classified into the following categories based on their duration: long-term SSE (L-SSE) lasting several months to years, and short-term SSE (S-SSE) lasting several days to weeks (Obara and Kato 2016). A remarkable feature of SSEs is their synchronization with other seismic phenomena, including low-frequency tectonic tremors (Rogers and Dragert 2003), very low-frequency earthquakes (VLFEs) (Ito et al. 2007), and repeating earthquakes (REs) (Gardonio et al. 2018). Synchronized events between SSEs and tremors are occasionally referred to as episodic tremor and slip (ETS) (Rogers and Dragert 2003).
The Nankai subduction zone, where the Philippine Sea plate subducts beneath the Amurian plate with a velocity of 63-68 mm/year (Miyazaki and Heki 2001), is one of the most active regions for SSEs worldwide (Fig. 1). Both L-and S-SSEs occur at the down dip extension of the rupture zones of large megathrust earthquakes (Obara and Kato 2016). L-SSEs with M w 6.8 and 6.6 have been observed by GNSS twice at Tokai, the eastern end of the Nankai subduction zone, over the past 28 years (e.g., Miyazaki et al. 2006;Ozawa et al. 2016) (Fig. 1). Other major L-SSEs occur approximately every 6 years beneath the Bungo Channel, located between the Shikoku and Kyushu Islands, with M w 6.8-7.0 (e.g., Hirose et al. 1999;Ozawa et al. 2020) (Fig. 1). In addition to these major L-SSEs, minor L-SSEs have been detected in the entire Nankai subduction zone between Kyushu and Tokai (Kobayashi 2017;Takagi et al. 2019).
S-SSEs have also been widely observed throughout the Nankai subduction zone (Obara et al. 2004;Hirose and Obara 2005;Nishimura 2014). S-SSEs occurring from Shikoku to Tokai coincide with tremors (Obara 2002) and have been analyzed using tiltmeters (e.g., Obara et al. 2004) and strainmeters (e.g., Kobayashi et al. 2006), since the early 2000s. Previous studies reported that these S-SSEs last for several days with approximately M w 6.0 (e.g., Obara et al. 2004;Hirose and Obara 2005). These studies manually detected S-SSEs and estimated their duration based on coincident tremors. In the last decade, S-SSEs in this region have been detected and analyzed using GNSS data (e.g., Nishimura et al. 2013;Kano et al. 2019;Sakaue et al. 2019). GNSS data are advantageous for analyzing S-SSEs because Japan has a dense GNSS network since the mid-1990s, with a signal-to-noise ratio that is rather stable over seasons and years as compared with tiltmeters and strainmeters. However, the duration of S-SSEs from Shikoku to Tokai has not been estimated systematically using GNSS data because of its low signalto-noise ratio. S-SSEs with approximately M w 6.0 have been detected around Kyushu (Nishimura 2014), but their duration has never been estimated. Episodic activations of VLFEs and REs have been reported in offshore Kyushu (e.g., Asano et al. 2015;Uchida et al. 2020), and their migration occasionally synchronizes with that of L-SSEs (Uchida et al. 2020). However, the relationship between S-SSEs and these seismic phenomena in Kyushu remains unclear.
SSEs have recently been detected in the up dip extension of the megathrust source region (Fig. 1) owing to improvements in seafloor geodetic measurements. Araki et al. (2017) observed changes in the pore pressure at seafloor borehole observatories off the Kii Peninsula and indicated that these changes were related to SSEs because they coincided with swarms of tremors. Yokota and Ishikawa (2020) detected year-long transient displacements using GNSS-Acoustic (GNSS-A) measurements from offshore Cape Muroto and suggested that these SSEs were accompanied by low-frequency events, including VLFEs (Asano et al. 2015).
Several methods for detecting S-SSEs using GNSS data have been developed. Nishimura et al. (2013) extracted transient signals from a single component of a GNSS time series, approximately parallel to the plate convergence, and estimated the fault models of S-SSEs. This method successfully detected S-SSEs that occur in and around the GNSS network (Nishimura et al. 2013;Nishimura 2014); however, it did not estimate the duration of S-SSEs. Distribution of slow slip events and the related phenomena in the Nankai subduction zone. The solid and dashed lines represent the Nankai Trough and iso-depth contours of the upper surface of the subducting Philippine Sea plate, respectively (Iwasaki et al. 2015). The pink arrow represents relative plate motion between the Philippine Sea and the Amurian plates (Miyazaki and Heki 2001). The gray shaded area represents the ≥ 2 m slip area of the 1944 Tonankai and the 1946 Nankai earthquakes (Sagiya and Thatcher 1999). Blue and green shaded areas represent the slip areas of L-SSEs (Kobayashi 2017;Takagi et al. 2019) and S-SSEs (Nishimura et al. 2013;Nishimura 2014), respectively. Red dots and orange shaded areas represent epicenters of tectonic tremors (Maeda and Obara 2009;Obara et al. 2010) and source areas of shallow VLFEs (Baba et al. 2020) Rousset et al. (2017) developed the geodetic matched filter method by using two horizontal components of the GNSS velocity time series (time derivative of the coordinate time series) with the spatial pattern of S-SSE displacements considered. This method successfully estimated the duration and magnitude of S-SSEs, but did not distinguish between fault size and slip amount.
As previously mentioned, the duration of S-SSEs from Shikoku to Tokai was estimated based on the accompanying tremors. This is not applicable to all S-SSEs in the Nankai subduction zone because S-SSEs do not always accompany tremors, as observed around Kyushu. The duration of S-SSEs around Kyushu has not been estimated yet and its regional characteristic along the entire Nankai subduction zone remains unclear. However, the duration as well as fault models of SSEs are important parameters for characterizing SSEs. Therefore, we developed an S-SSE detection method that builds on previously developed methods (Nishimura et al. 2013;Rousset et al. 2017). Using this method, we detected S-SSEs and estimated their causative fault models and duration in the Nankai subduction zone, using all available GNSS data. We also attempted to compare the detected S-SSEs to other phenomena, including tremors, VLFEs, and REs, and verified the regional characteristics of their slip amount and duration.

Data and preprocess
This study used the daily coordinates of 734 GNSS stations located in the Nankai subduction zone from July 1, 1996 to June 19, 2020 (Fig. 2). These stations include not only the GNSS Earth Observation NETwork system (GEONET) stations operated by the Geospatial Information Authority of Japan (GSI), but also stations operated by the Japan Coast Guard, Kyoto University, Kobe University, Kochi University, Kyushu University, Japan Crustal Activity Science Consortium, and International GNSS Service. The daily coordinates were calculated using the precise point processing with ambiguity resolution strategy of Gipsy-X version 1.3 (Bertiger et al. 2020). These coordinate data are based on the International Terrestrial Reference Frame 2014 (ITRF2014) (Altamimi et al. 2016).
We removed noise and signals unrelated to S-SSEs to enhance detectability prior to S-SSE detection. Firstly, we corrected the artificial offsets caused by station maintenance, by removing the difference in the 10-day averages before and after the offsets using the maintenance log provided by the GSI. Secondly, we removed For the three large earthquakes; the 2004 M w 7.5 off the Kii Peninsula earthquake, the 2011 M w 9.0 Tohoku-Oki earthquake, and the 2016 M w 7.0 Kumamoto earthquake, we corrected both the coseismic and postseismic deformations. We removed coseismic offsets using the same method as that of the artificial offsets for the other earthquakes where coseismic displacements were recorded at GEONET stations (GSI 2019). We estimated the secular velocity in a stable period listed in Additional file 1: Table S1, before removing the coseismic and postseismic deformations for the three large earthquakes. For the 2011 M w 9.0 Tohoku-Oki earthquake, the coseismic and postseismic deformations were modeled using Eq. (1), based on Tobita (2016): where t is the number of days since the earthquake, x 1 (t) is a displacement component, a 1 , b 1 , and c 1 represent the amplitudes of the coseismic and postseismic deformation, H(t) represents the Heaviside step function, v 1 is the secular velocity, and d 1 is the intercept at t = 0 . We defined t = 0 as the day of earthquake occurrence. In the case of the other two large earthquakes, we fit Eq. (2) expressed as: where x 2 (t) is a displacement component, a 2 and b 2 represent the amplitude of the postseismic deformation and its time constant, respectively, c 2 is the amplitude of coseismic deformation, v 2 is the secular velocity, and d 2 is an intercept at t = 0 . We fit the above functions for each period in Additional file 1: Table S1 (Additional file 1: Figs. S1ac). Thirdly, we corrected the offsets caused by the 2011 Kirishima Shinmoe-dake eruption (Nakao et al. 2013) and the 2015 dike intrusion in the Sakurajima volcano (Hotta et al. 2016) using the same method as that used for the artificial offset. Fourthly, we removed long-term trends including seasonal oscillations by calculating the 91-day moving averages (Additional file 1: Figs. S1a-c). Finally, we removed common-mode errors, presumably caused by an estimation error of the satellite orbits and the reference frame, by applying a spatial filter (Wdowinski et al. 1997). The spatial filter effectively reduced the time series scatter (Additional file 1: Fig. S1d).

SSE detection method
Our method for detecting S-SSEs can be divided into three steps. Firstly, we detected transient displacement events from GNSS time series by calculating weighted averages of correlation coefficients. Secondly, we categorized the transient events into S-SSEs using the estimated results of their fault models and duration. Finally, we applied a bootstrap method to estimate the duration error margins for each detected S-SSE. The details of each step are described in the following sections.

Detection of transient displacement events
Firstly, we calculated correlation coefficients between preprocessed GNSS data and a synthetic template while shifting the time window to detect the time of occurrence of the transient events in each GNSS time series. It was assumed that the temporal evolution function of S-SSEs can be expressed as a ramp function. Because the 91-day moving average of the GNSS time series was removed, static offsets are expected to vanish in the long-term average. Therefore, we assumed a synthetic template as the sum of a ramp function and a linear function that is zero at both ends of a time window (Fig. 3a). In this step, the duration of the ramp and the total length of the time window were set to 4 and 121 days, respectively. The template length and the duration of the ramp was chosen based on the recurrence interval and the mode value of the S-SSE duration reported by Sekine et al. (2010). The template length is taken to be as long as possible but not to include multiple S-SSEs. This calculation was executed for each horizontal component at all stations. Although the correlation coefficients of the single component significantly fluctuated, we found that the negative peak of correlation coefficients around early February 2007 approximately corresponded with the tremors near the station (Fig. 3b). The amplitude and sign of the correlation coefficient approximate those of the transient displacement, therefore, the negative peak in Fig. 3b indicates a large southward motion. This suggests that a correlation coefficient is an effective tool for detecting transient signals in GNSS time series. Secondly, we calculated the weighted average of correlation coefficients in order to consider the spatial pattern of displacements predicted from SSEs, originally proposed by Rousset et al. (2017). Daily weighted averages of correlation coefficients were calculated for each rectangular sub-fault located at 5-50 km depth range on the Philippine Sea plate interface ( Fig. 2) (Iwasaki et al. 2015). The weighted average of correlation coefficients at sub-fault f is defined by where C s,i (t) is the correlation coefficient of i th component, that is, east-west (EW) and north-south (NS) components, at station s . The weight G f s,i in Eq. (3) is set for each sub-fault f , station s , and component i and is defined by where g f s,i is the predicted surface displacement of i th component at station s caused by a unit slip at sub-fault f , and m is a constant for the minimum weight in normalization. We assumed that the slip of the sub-faults was approximately opposite to the plate convergence direction (N55°W; Miyazaki and Heki 2001). We computed surface displacements using Green's function for an elastic half-space (Okada 1992). This weighting process enabled us to emphasize the correlation coefficient of a component at a station with a larger displacement predicted from the slip on the plate interface (Eq. 3). In Fig. 4, sub-faults with high weighted average of correlation coefficients are spatially coincident with coherent south-eastward displacements and tremors. This indicates that the weighted average of correlation coefficients is effective for detecting the timing and location of a transient displacement event. However, sub-faults with high weighted averages of correlation coefficients were occasionally distributed extending toward a trench normal direction without a clear peak (Additional file 1: Fig. S2a). This is attributed to neglecting the contribution of some stations that potentially constrain the fault location, even if no displacements are predicted from a sub-fault at the stations. For example, if the actual fault locates inland, the weighted averages in offshore sub-faults sometimes can be high because a sub-fault located offshore sometimes predicts a displacement pattern partially similar to those located inland. In addition, very small displacements can be predicted even if the stations are close to the fault, depending on the fault dip and depth. Thus, we introduced m in Eq. (4) as the minimum weight after normalization. This term can maintain the minimum weight of the correlation coefficient for all components. However, too large m reduces the amplitude of the weighted average of correlation coefficients and hides the peak of the weighted average of correlation coefficients relating to an S-SSE. We empirically decided the minimum weight to be 0.02 to enhance the spatial resolution in the trench normal direction and to detect many events.
Thirdly, we extracted peaks that were larger than the threshold from a weighted average of correlation coefficients. The threshold of the weighted averages of correlation coefficients was set to 0.0215, which was the sum of the average and standard deviation for all of them. This gentle threshold extracts many peaks potentially related to S-SSEs including those related to non-SSE events. In the next step, a further examination is conducted to judge whether the extracted displacement event is an actual S-SSE. The extracted peaks can be regarded as the middle date of the duration and location of the transient displacement event. We then selected the maximum peak within a range of 100 km and 41 days from the extracted peaks of the weighted average of correlation coefficients satisfying the threshold, to avoid overlapping events. In this study, we regarded the date and location of the maximum peaks as those of the transient displacement events.

Categorization of transient displacement events into S-SSEs
After detecting the transient displacement events, we estimated their rectangular fault models and also estimated the duration of the transient displacement events from the stacked GNSS time series. We then applied criteria to categorize them as S-SSEs.
Firstly, we calculated displacements assuming a duration from 1 to 40 days. Displacements for each duration were calculated by fitting a synthetic template to 121-day GNSS coordinates, while varying the duration of the ramp function. We then estimated the fault models from the displacements for each duration. For the fault model estimation, we used the nonlinear inversion method of Matsu'ura and Hasegawa (1987), assuming a rectangular fault model in an elastic half-space (Okada 1992) following Nishimura et al. (2013). We constrained the fault model by the interface between the subducting Philippine Sea plate and the overriding continental plates (Iwasaki et al. 2015), computing the depth, strike, and dip based on the longitude and latitude of the center of the fault. Therefore, the estimated parameters are longitude and latitude of the fault center, length, width, rake, and slip amount. In addition to these parameters, we also estimated the translation parameters representing arbitrary common offsets at all stations, which are often recognized in the detected displacements in the EW, NS, and up-down (UD) directions. Thus, nine parameters were estimated for the fault model. In the inversion, we used the longitude, latitude, and rake of the sub-fault as initial values for the transient displacement event. The other initial parameters, length, width, and slip amount, were set as 50 km, 30 km, and 0.01 m, respectively. We also set the initial standard deviation of longitude, latitude, length, width, rake, and slip amount as 3.0°, 3.0°, 20 km, 10 km, 15°, and 99 m, respectively. Once the fault model is estimated, we calculate a reduction of chi-square to assess the fault estimation result. Previous studies (Nishimura et al. 2013;Nishimura 2014) defined the reduction of chi-square as: where O j and σ j represent the observed displacement of j th data and its error, C all j is the calculated displacement from both the estimated fault model and the translation of the coordinate, and C trans j is the calculated displacement from the translation. The reduction of chi-square is larger when the fault model accurately reproduces the observed displacement by considering the contribution from the translation. We used the reduction of chi-square in the final process to categorize the event into S-SSEs.
(5) Secondly, we estimated the duration of the event based on the method proposed by Rousset et al. (2017). Before the duration estimation, we calculated the weighted stacking of the horizontal GNSS coordinates in a 121day time window, with the center date being the date of the detected event (Fig. 5a). Previous studies stacked weighted GNSS time series by either displacement signals (Rousset et al. 2017(Rousset et al. , 2019 or the noise level of the time series (Miyaoka and Yokota 2012;Nishimura 2021), whereas this study used both signal and noise levels for weighting (Eq. 6): In Eq. (6), x s,i (t) and x d (t) represent the GNSS coordinate data of component i at station s and the weighted stack of coordinate data, respectively. The weight c d s,i is defined as where u d s,i is the displacement computed from the estimated fault assuming an event duration of d days, σ s,i is the noise level of the GNSS time series of component i at station s , σ is the average of σ s,i , and A d is a constant depending on the duration. We selected the GNSS coordinates of the first and last 30 days in the 121-day window and regarded their standard deviation as the noise level of the GNSS data. The coordinates were stacked individually in descending order with a weight factor c d s,i for all used data and then the correlation coefficient between the stacked time series and a synthetic template was calculated. This procedure enabled us to determine the number of stacked coordinates using the maximum correlation coefficient in a set duration (cf. Miyaoka and Yokota 2012). For instance, the correlation coefficient maximized when 17 components of the coordinates were stacked in the case of Fig. 5a. We then conducted this process for varying duration from 1 to 40 days. We compiled and maximized the correlation coefficient for all duration and determined the event duration (Fig. 5b). Finally, we applied the following four criteria to the estimated fault models and duration to categorize the transient displacement events into S-SSEs. We assumed that the SSEs were reverse faulting and that their slip direction was opposite to plate convergence. Thus, we extracted the events whose fault model had rake and slip azimuth (that is, strike minus rake) in the range of 20° to 160° and N100° E to N170° E, respectively. In addition, we extracted an event which had a correlation coefficient, to the stacked time series, larger than 0.4 in order to retain a significant signal of transient displacement caused by S-SSEs. We regarded the events satisfying the above three criteria as S-SSEs and further categorized them into classes 1 and 2 based on a reduction of chi-square calculated in the fault model estimation process. The reduction of chi-square is ≥ 150 and 50 -150 for class 1 and 2 SSEs, respectively. We subjectively decided these thresholds by comparing the detected events to the previously detected S-SSEs and tremors. Events that do not satisfy these criteria or whose reduction of chi-square is smaller than 50 are categorized as class 3 events. We also manually categorized several class 3 events. These include tectonic transients, including the 2000 Kozushima-Miyakejima seismo-volcanic event (Nishimura et al. 2001) and large regular earthquakes. They also include events whose fault model is estimated further offshore than the trough axis. Additionally, we chose the event with the largest reduction of chi-square for overlapped duration events.

Estimation of duration error
The duration of S-SSEs detected in GNSS data has a large uncertainty because of the small displacements caused by S-SSEs. Therefore, assessing the uncertainty in the duration for further analysis is essential. We estimated the estimation error of the duration of the detected S-SSEs by applying a bootstrap method. In the duration estimation process, the GNSS time series with the first to N th largest weights were stacked until the correlation coefficient was maximized. The value of N maximizing the correlation coefficients ranges from 1 to 819 for all the detected events and its median is 195. In this step, we randomly resampled GNSS time series N times from those with the first to N th largest weights, allowing resampled data overlap. We then estimated the duration from the resampled stacked time series using the same method as that in the duration estimation process. We iterated this process 2,000 times and obtained the frequency distribution of the duration (Fig. 5c). We used these frequency distributions for a detailed analysis of the duration in the following chapters. We also computed the 70% confidence intervals (that is, an interval between 15 and 85th percentiles) for each frequency distribution (Fig. 5c) and used it as the duration error.

Result
We detected 143 class 1 SSEs and 141 class 2 SSEs from February 1, 1997 to January 31, 2020 using our new method. The detectability of our method is approximately equal to that of Nishimura et al. (2013) because we detected 153 SSEs in a region east of 132°E from February 1997 to January 2012 comparable to the 155 detected by Nishimura et al. (2013) in the same region and period.
Most of the detected S-SSEs shown in Fig. 6a are distributed at a depth of 30-40 km, and almost coincide with tremors, especially from Shikoku to Tokai (Fig. 1). This distribution is similar to that in previous studies (Sekine et al. 2010;Nishimura et al. 2013;Nishimura 2014), but some of the detected SSEs are located shallower or deeper than the previously recognized ETS. It is notable that 10 or more S-SSEs were detected in the up dip of the rupture regions of the megathrust earthquakes ( Figs. 1 and 6a) and their clusters were found offshore Kyushu, offshore Cape Muroto, and offshore the Shima Peninsula. Several S-SSEs offshore Kyushu were contiguous with REs (Uchida et al. 2020) and VLFEs (Baba et al. 2020). We discuss the details of these S-SSEs and related phenomena in the "Examples of individual events" section. One of three S-SSEs offshore Cape Muroto was detected alongside the shallow SSE detected by offshore geodetic measurements (Yokota and Ishikawa 2020) and was also accompanied by VLFEs (Baba et al. 2020). However, the other two S-SSEs offshore Cape Muroto are possibly false-detections because we were unable to find any coincidence with other phenomena, including VLFEs and SSEs. Moreover, S-SSEs detected offshore the Shima Peninsula coincided with tremors that occurred in the down dip extension of the megathrust earthquakes. This suggests an inaccurate estimation of the source location caused by the poor spatial resolution in the trench normal direction in our method, although the method successfully detected transient displacement events. Furthermore, we detected 10 or more S-SSEs scattered at a depth of > 40 km, a depth at which no SSEs have been previously identified (Sekine et al. 2010;Nishimura et al. 2013). They may also be falsely detected events because their scattered distribution is inconsistent with SSEs characteristics, which repeat at similar locations with a relatively short recurrence interval (Sekine et al. 2010).
Although the frequency distribution of the moment magnitude of the detected S-SSEs (Fig. 6b) is approximately similar to that of previous studies, the duration is rather different from that of previous studies (e.g., Sekine et al. 2010). Most of the S-SSEs detected by Sekine et al. (2010) had a duration of ≤ 10 days, whereas 58 events detected in this study lasted for more than 20 days (Fig. 6c). However, the ranges of the 70% confidence interval for the estimated duration were rather wide (that is, 10 to 20 days) (Additional file 1: Fig. S3). A comparison of the duration of 29 events detected in both this study and Sekine et al. (2010) demonstrates that the duration of 23 events is concordant with each other within the 70% confidence interval. We concluded that the difference in duration between our study and the previous studies using tiltmeters (Sekine et al. 2010) is insignificant within the uncertainty of the duration estimated in our study. This indicates the importance of considering the error for the detailed duration analysis of SSEs, especially S-SSEs in the Nankai subduction zone.
We studied the relationship between the moment and duration of the detected S-SSEs (Additional file 1: Fig.  S4). Ide et al. (2007) proposed that the seismic moment of slow earthquakes, which include SSEs, tremors, and VLFEs is proportional to their duration. The detected S-SSEs as well as L-SSEs of Takagi et al. (2019) are approximately consistent with the scaling law for slow earthquakes (Ide et al. 2007). However, discussing the moment-duration relationship of detected S-SSEs in detail is complicated because of the large uncertainty and limited range of the estimated duration and moment.

Examples of individual events
The fault model of the S-SSE in Tokai corresponding to the event shown in Figs. 3, 4 and 5, reasonably reproduced the observed GNSS displacements (Fig. 7a). The magnitude and duration of the event are similar to those in previous studies (Sekine et al. 2010;Nishimura et al. 2013). This event was coincident with tremors and was recognized as an ETS in previous studies (e.g., Sekine et al. 2010). An S-SSE in western Shikoku (Fig. 7b) was estimated to last longer than that in previous studies (Hirose and Obara 2005;Sekine et al. 2010). This event can be interpreted as an accelerated part of an L-SSE because it occurred during the most active period of the 2003 Bungo Channel L-SSE (Yoshioka et al. 2015). The SSEs in eastern Shikoku have a right-lateral strikeslip component as compared to those in other regions (Fig. 7c, Additional file 2). This is probably caused by the difference in the angle between the iso-depth contours of the subducting plate and its convergence direction. A spatio-temporal comparison between detected S-SSEs and tremors suggests that most of the S-SSEs between western Shikoku and Tokai accompanied tremors, though some large tremor bursts occurred without detected SSEs around the Shima Peninsula (Additional file 1: Fig.  S6). We speculate that these tremors are accompanied by small SSEs and their signals are too small to be detected by GNSS data.
The S-SSE in offshore Kyushu in January 2010 is approximately coincident with the spatial and temporal burst of VLFEs (Baba et al. 2020) and REs (Uchida et al. 2020) (Fig. 7d). We also detected an S-SSE in December 2017 with a major VLFE burst though we could not study an RE activity because the analysis period of Uchida et al. (2020) ended in mid-2015 (Additional file 3). Focusing on the location of each phenomenon, VLFEs and REs appear to occur in the up dip extension of these S-SSEs.

136˚E
137˚E  However, we could not discuss a detailed spatial relationship between these S-SSEs and other phenomena because of the poor resolving power of GNSS data for offshore regions. We also compared the temporal relationships of each phenomenon during the January 2010 and December 2017 events ( Fig. 8 and Additional file 1: Fig. S7). A major activity of VLFEs started approximately with the January 2010 S-SSE and was sustained after the S-SSE, while a burst of REs occurred later during the S-SSE (Fig. 8). Therefore, we concluded that the January 2010 S-SSE has some relationship with VLFEs and REs. Uchida et al. (2020) indicated the occurrence of an S-SSE due to the sudden increase in the number of VLFEs (Asano et al. 2015) and the cumulative slip of the REs. The January 2010 S-SSE manifests the first geodetic evidence on S-SSEs related to VLFEs and REs in offshore Kyushu. The bursts of VLFEs synchronizing with the January 2010 and December 2017 S-SSEs in offshore Kyushu are more active than those around a depth of 30-40 km (Baba et al. 2020). This indicates that S-SSEs in offshore Kyushu activate more VLFEs than the deep ETS in the Nankai subduction zone. However, it is difficult to discuss details of their relationship because of the limited spatial and temporal resolution of the estimated S-SSEs. The relationship is expected to be clarified by improving the resolution using seafloor geodetic measurements in the near future.

Regional characteristics of overall S-SSEs
We evaluated the regional characteristics of the detected S-SSEs by calculating the cumulative slip and duration, and the average duration and slip rate (that is, the quotient of the cumulative slip divided by the cumulative duration) over 23 years. This is the first attempt to examine the duration and slip rate of S-SSEs estimated from geodetic data over the Nankai subduction zone. We used a Monte Carlo method to quantitatively compare them with their uncertainties. The calculation was performed at each grid point with an interval of 0.02°, in both the longitudinal and latitudinal directions. Firstly, we randomly selected the duration d i of the i th SSE from the frequency distribution obtained for each SSE in the duration error estimation process (Fig. 5c). Secondly, we drew the slip amount s i of the i th SSE from a normal probability distribution. The mean of the normal distribution is the estimated slip of the i th SSE corresponding to the selected duration d i and the standard deviation of the distribution is the error estimated by the nonlinear inversion in the second step of the developed method. Our detection method estimated the fault model of an S-SSE for each duration. The selection of the duration and slip amount was performed for all 284 S-SSEs. Thirdly, we calculated the cumulative duration and slip by adding d i and s i of the i th SSE to the grid point within the surface projection of the SSE fault. We also computed the average duration by dividing the cumulative duration by the cumulative number, while the average slip rate was also calculated by dividing the cumulative slip by the cumulative duration. We iterated this calculation process 10,000 times and obtained the frequency distribution of the cumulative number, slip, duration, average duration, and average slip rate (Additional file 1: Figs. S8-11). Finally, we adopted averages as the representative values of each parameter and twice the standard deviations as their uncertainties ( Fig. 9 and Additional file 1: Fig. S12) because the frequency distribution can be approximately regarded as a normal distribution (Additional file 1: Figs. S8-11).
The distribution of the cumulative number shows four major clusters in western, central, eastern Shikoku, and around Ise Bay (Fig. 9a). These clusters have been suggested to be the most active S-SSEs in previous studies (e.g., Sekine et al. 2010;Nishimura et al. 2013). Conversely, there were no notable clusters around Kyushu, but S-SSEs were zonally distributed along the coastline (Fig. 9a). This is slightly different from the results of Nishimura (2014) who detected clusters in offshore northern Kyushu and offshore Tanegashima Island (approximately 131°E and 30.5°N). We did not assign any sub-faults south of 31°N (Fig. 2) because the offshore Tanegashima Island region is beyond the scope of this  (Fig. 7d). a Stacked time series for the January 2010 S-SSE. Red dots and a green shaded area represent a stacked GNSS time series and the estimated duration of the S-SSE, respectively. b Activity of VLFEs (Baba et al. 2020) and REs (Uchida et al. 2020) in the region of 130-133°E and 30-32.5°N. Brown stars show a magnitude for REs, and yellow bars represent the daily count number of VLFEs study. However, it is unclear why our study did not show the cluster in offshore northern Kyushu. We speculated that the different distribution of the cumulative number around Kyushu between our study and Nishimura (2014) was caused by the different sensitivities depending on the duration of the detected SSEs. We also removed longterm trends from GNSS coordinates during preprocessing, unlike Nishimura (2014), therefore, the catalog of   Nishimura (2014) probably records SSEs with a longer duration, and our method did not detect them as S-SSEs. The distribution of the cumulative slip also shows a similar pattern to that in a previous study (Nishimura et al. 2013) from western Shikoku to Tokai (Fig. 9b). Three clusters, central and eastern Shikoku and around Ise Bay, slipped by 0.2 m over 23 years and their location corresponds to the region where the coupling ratio is small, during -2009(Nishimura et al. 2018. The cluster in western Shikoku recorded the largest slip among the four clusters. In western Shikoku, the coupling ratio is also larger than that in other S-SSE source regions (Nishimura et al. 2018). The relationship between the cumulative slip and coupling ratio in western Shikoku appears inconsistent because active SSEs generally reduce the coupling ratio over a period of SSE recurrence. Another factor reducing the coupling ratio is a stable creep rate. If the creep rate in western Shikoku is much lower than that in the other S-SSE regions, it explains both the high cumulative slip of S-SSEs and high coupling ratio. We speculate that the stress shadow caused by the high coupled area which locates up dip of the SSE cluster in western Shikoku (Nishimura et al. 2018) reduces the stable creep and keeps the high coupling ratio despite the large cumulative slip of S-SSEs.
The distribution of the cumulative duration shows a spatial pattern different from that of the cumulative slip (Fig. 9c). The cumulative duration in eastern Shikoku is longer than 400 days, the longest in the Nankai subduction zone (Fig. 9c). In addition, the cumulative duration in Kyushu (> 200 days) was comparable to that in the major clusters in Shikoku and around Ise Bay (Fig. 9c). Although the estimated duration has a large uncertainty owing to a low signal-to-noise ratio, the total duration of S-SSEs has a different regional characteristic from the total slip in the Nankai subduction zone. The average duration also shows a regional characteristic. The average duration around Kyushu is the longest among the Nankai subduction zone (Fig. 9d). The average duration in western Shikoku and around Ise Bay is shorter than that in eastern Shikoku.
The spatial distribution of the average slip rate shows that the slip rate is the highest in western Shikoku and gradually decreases toward eastern Shikoku (Fig. 9e). We compared these estimated values at the points aligned on the 30 km iso-depth contour (Iwasaki et al. 2015) (Additional file 1: Fig. S5) to assess the regional characteristics at the same depth along the Nankai Trough (Fig. 10). There appears to be a clear regional difference in the average slip rate between western Shikoku (approximately 132.5°E) and eastern Shikoku (approximately 134°E) (Fig. 10d). This difference is statistically significant because the error bars of the 95% confidence interval do not overlap (Fig. 10d). The average slip rate around Kyushu (between 131°E and 132°E) is roughly similar to that in eastern Shikoku (Figs. 9e and 10d). The slip rate around Ise Bay (around 136.7°E) is larger than that in eastern Shikoku and smaller than that in western Shikoku, although these differences are not statistically significant (Figs. 9e and 10d).
We compared the spatial pattern of the average slip rate of S-SSEs to tremors. Daiku et al. (2018) estimated the conversion factors between the apparent moment of tremors and the seismic moment of S-SSEs, the inverse of which represents the excitation efficiency of tremors by S-SSEs, in the Nankai subduction zone. They showed that the conversion factor in western Shikoku was smaller than that in eastern Shikoku (Fig. 10d). S-SSEs in western Shikoku enhance tremor activities more than those in eastern Shikoku (Daiku et al. 2018). We infer that the conversion factor roughly anticorrelates with the average slip rate estimated in this study and a higher slip rate enhances tremors, as suggested by Wech and Bartlow (2014) (Fig. 10d). However, this correlation is not retained throughout the Nankai subduction zone because little or no tremors are observed in areas with rather high slip rates, including Ise Bay and Kyushu (Fig. 10d). The conversion factor is low around the Shima Peninsula while the average slip rate in this region is not significantly higher than that in the other regions (Fig. 10d). We speculate that the excitation efficiency of tremors relates to both the slip rate of S-SSEs and fault properties, including a high pore fluid pressure (e.g., Shelly et al. 2006).
Finally, we discuss the control of the average slip rate of the S-SSEs. Previous studies have suggested that the occurrence of SSEs is controlled by the lateral change of the subducting plate interface geometry (Mitsui and Hirahara 2006), the subduction of the seamount (Sun et al. 2020), the fluid pressure, and frictional properties (Saffer and Wallace 2015). The plate interface geometry seems to be the dominant factor because the slip rate in a region with a high dip (e.g., around Kyushu) is significantly lower than that in a region with a low dip (e.g., western Shikoku) (Figs. 9e and 10d) and a correlation coefficient between the slip rate and dip angles of the plate interface (Iwasaki et al. 2015) is −0.85 (Additional file 1: Fig. S13). However, plate geometry cannot completely explain the slip rate characteristics. For example, the slip rate in eastern Shikoku is similar to that around Kyushu, while the dip of the plate interface in eastern Shikoku is significantly smaller than that around Kyushu. Therefore, we conclude that the activity of S-SSEs in the Nankai subduction zone is governed by multiple factors, including the oceanic plate interface geometry.

Conclusions
We successfully detected 284 S-SSEs by applying our newly developed method in the Nankai subduction zone between February 1997 and January 2020. Careful preprocessing, correlation coefficients between GNSS time series and a synthetic template, and weighting by displacements predicted from assumed faults are key processes in extracting small transient signals from noisy observed data. We detected S-SSEs synchronizing with VLFEs and REs in offshore Kyushu. Although their detailed spatio-temporal relationship remains unclear, using both inland and offshore geodetic measurements for analyzing and checking either geodetic or seismic phenomena, may assist in revealing their relationship in offshore Kyushu. The quantitative calculation considering the errors of both the fault models and duration enables us to discuss the averaged regional characteristics of SSEs. In addition to the cumulative number and slip discussed in previous studies, we clarified the regional characteristics of the cumulative duration and the average slip rate in the Nankai subduction zone. The average slip rate in western Shikoku is approximately twice as large as that in eastern Shikoku and Kyushu. We speculate that the subducting plate interface geometry is a mechanism for the regional differences in S-SSEs. Statistical examinations with SSE observations increased by our enhanced detection method will help clarify the regional characteristics of SSEs worldwide. the approximate location of Kyushu, western Shikoku, eastern Shikoku, the Kii Peninsula, the Shima Peninsula, and Tokai, respectively. Note that the calculation points experiencing less than five events are shaded. a Cumulative slip and its error. b Cumulative duration and its error. c Average duration and its error. d Average slip rates and their error. Green triangles and blue squares represent a dip angle of the subducting plate (Iwasaki et al. 2015) and the conversion factor between the tremor moment and S-SSE moment (Daiku et al. 2018)