Time-dependent modeling of slow-slip events along the Nankai Trough subduction zone, Japan, within the 2018–2023 period

According to a 2013 report by the Earthquake Research Committee of Japan, it was estimated that the probability of the occurrence of the next Nankai earthquake within the next three decades is 70–80%. Therefore, to realize hazard estimation, it is crucial to estimate the state of the plate interface in the Nankai Trough subduction zone. In this study, we focused on a transient from 2018 in the subduction zone of the Nankai Trough, southwest Japan, detected by the Global Navigation Satellite System (GNSS) network. Thereafter, we estimated the spatiotemporal evolution of the slip at the plate interface by subjecting the GNSS position time series to time-dependent inversion. The results obtained showed that a long-term slow-slip event (l-SSE) possibly occurred on the Kii Peninsula in 2020. The Kii-Chan-nel l-SSE (M w 6.3), with an irregular recurrence interval, was observed from 2019 to 2022. Additionally, the Central Shi-koku l-SSE (M w 6.5) was observed from 2019 to 2023, and there appeared to be a correlation between the Central Shi-koku l-SSE and the northwestern Shikoku short-term slow-slip event (s-SSE). l-SSE occurred from July 2018 to August 2019 in the northern Hyuga-nada and Bungo channel, and in late 2018, another l-SSE occurred in southern Hyuga-nada before the May 2019 Hyuga-nada earthquake. Further, after the 2018–2019 events, the southern Hyuga-nada l-SSE occurred from mid-2020 to early 2021 and in January 2023 and is still ongoing, consistent with the expected recurrence interval. The seemingly transient slip off the Ohsumi Peninsula, detected after the M w 5.7 earthquake of October 2022, continued until July 2023 with M w 6.1. Furthermore, after the 2019 Tanegashima earthquake, an l-SSE occurred for approximately 4 years. The moment magnitude (M w ) of the afterslip of the Tanegashima earthquake was estimated to be 6.7, larger than that of the main shock (M w 6.4). We also noted that the transient slip off Tane-gashima included four possible s-SSEs that occurred in 2020, 2021, 2022, and 2023. Therefore, the sporadic aseismic slips along the Nankai Trough changed the stress state of the areas neighboring the aseismic slip zones in favor of the interplate slip.


Introduction
The Philippine Sea Plate subducts from the Nankai Trough beneath the Amur Plate in an N57°W direction at 5 to 6 cm/year (Fig. 1) (Kreemer et al. 2014).Owing to this subduction of the Philippine Sea Plate, there have been repeated megathrust earthquakes in the Nankai Trough subduction zone (e.g., Kumagai et al. 1996).Further, given that it is already 70 years since the last Nankai earthquake, it is thought that the Nankai Trough subduction zone has accumulated enough strain and is capable of releasing an earthquake with a moment magnitude (M w ) up to 8. Additionally, in 2013, it was reported that the probability of the occurrence of the next Nankai earthquake within the next 30 years is 70-80% (Earthquake Research Committee 2013).Given that slowslip events (SSEs) change the stress state and, in some cases, trigger large-scale earthquakes in nearby areas (e.g., Radiguet et al. 2016), it is crucial to monitor SSEs in the Nankai Trough subduction zone for better hazard estimation.
The Nankai Trough subduction zone is characterized by shallow SSE (Araki et al. 2017;Yokota and Ishikawa 2020), locked, long-term SSE (l-SSE), and short-term SSE (s-SSE) areas from the shallow part to the downdip extension.The locked area is believed to be the source of the next Nankai earthquake (Fig. 1b).Further, as shown in Fig. 1b, the l-SSE area comprises several segments, including the Kii channel, east Shikoku, Central Shikoku, west Shikoku, Bungo Channel, northern Hyuga-nada, and southern Hyuga-nada segments (Kobayashi 2017;Takagi et al. 2019;Yoshioka et al. 2015).The duration of l-SSEs ranges from months to several years.Additionally, in the down-dip extension of the l-SSE area, there is an s-SSE area, with s-SSEs with durations in the order of several days, associated with low-frequency earthquakes (Obara 2002;Hirose and Obara 2005).
In 2018, the Global Navigation Satellite System (GNSS) network in Japan detected transient crustal deformation from Shikoku Island to Kyushu Island (Ozawa et al. 2020).In our previous study, we did not consider postseismic deformation after the 2011 Tohoku and 2016 Kumamoto earthquakes during the treatment of position time series; this resulted in erroneous trend estimates.In the present study, we corrected for the effect of viscoelastic deformation after these two earthquakes based on Suito (2017a, b), as we cannot ignore the effects of these deformations in trend estimates and in the long-term.Furthermore, we removed the displacements associated with s-SSEs, given that they contaminate the accurate estimation of l-SSE, especially in the Shikoku area.Thereafter, we expanded the analytical area reported by Ozawa et al. (2020) and estimated the spatiotemporal evolution of the slip on the plate interface since 2018  (Sagiya and Thatcher 1999), which is locked strongly and expected to slip in the future Nankai earthquake.The green line indicates the area of the plate interface adopted in the time-dependent inversion.The blue circles show the locations of the GNSS sites used in the time-dependent inversion, whereas the white circles show the locations of the GNSS sites whose position time series are shown in Fig. 2. The T-direction and S-direction represent the directions toward which the parameters of the parametric spline surface increase (Ozawa et al. 2020) based on the corrected position time series without viscoelastic relaxation and s-SSE offsets.

Analytical procedure
GNSS data were analyzed daily using Bernese GPS software (version 5.2) with precise ephemeris and earth rotation parameters (Takamatsu et al. 2023).We first removed the offsets of s-SSEs from the longitude, latitude, and height coordinates.To compute the s-SSE offsets, we used the fault parameters published by the Japan Meteorological Agency (JMA) and the National Institute of Advanced Industrial Science and Technology (AIST) of the Nankai Trough Earthquake Assessment Committee (e.g., JMA 2022) together with the catalog reported by Nishimura et al. (2013), which covers the 1996-2012 period.We removed s-SSE coseismic displacements on the onset dates of the events, given that s-SSEs last only several days and have much shorter durations than l-SSEs (several months to years).After that, we removed the viscoelastic deformations of the 2011 Tohoku and 2016 Kumamoto earthquakes based on Suito (2017a, b) (Additional file 1: Fig. S1).Next, we removed the linear trend from the position time series at stations for the periods that did not include the l-SSE transient by regressing a linear function from the position time series.We adopted a trend period for each region (Additional file 1: Fig. S2) based on the event table summarized by Takagi et al. (2019).Coseismic offsets for large earthquakes were estimated by subtracting coordinates averaged over 10-day periods before and after the earthquake and corrected from the detrended time series.Additional file 1: Fig. S3 shows estimated trend vectors, which show smooth spatial distribution, indicating the suitability of the adopted reference trend with an appropriate backslip estimation (Additional file 1: Fig. S4).
We applied time-dependent inversion (McGuire and Segall 2003) to east-west, north-south, and up-down components at 240 selected GPS sites (Fig. 1c) for the January 1, 2018-July 1, 2023 period.For this analysis, we averaged data over 3-day periods and used every 3 days.We also weighted the horizontal and vertical components in a 3:1 ratio for subsequent analysis, considering the extent of data scatter.
For time-dependent inversion, we used a spline surface to represent the fault interface (Ozawa et al. 2020) and adopted the geometric model of the Philippine Sea Plate proposed by Hirose et al. (2008).The range of the adopted plate interface model is shown in Fig. 1c.In our representation of the plate interface, we used 44 and 15 grid knots in the T-and S-directions, respectively, as shown in Fig. 1c.The spacings of the grid knots in the T-and S-directions were approximately 15 and 20 km, respectively.In the following analysis, we assumed Poisson's ratio of 0.25 and rigidity of 30 GPa.We calculated Green's function according to Yabuki and Matsu'ura (1992) and adopted the state x(t) in the time-dependent inversion: where u k and v k represent the slip and slip velocity, respectively, at the kth grid on the plate interface with a total of L grids.e EW , e NS , and e UD represent the common mode errors of the east-west, north-south, and up-down components in the position time series, respectively (Fukuda 2018).We incorporated spatial roughness into the transition equation (Ozawa et al. 2020) and estimated hyperparameters by maximizing the log-likelihood of the system.
We constrained the aseismic slip vector and slip velocity vector on the plate surface in a direction opposite that of the motion of the Philippine Sea Plate relative to the Amur Plate using the positivity method by Simon and Simon (2006).This direction of the Philippine Sea Plate relative to that of the Amur Plate was obtained from Kreemer et al. (2014).
We conducted a checkerboard test by inverting the synthesized data with a Gaussian error, with a 1-mm standard deviation added for the horizontal components and a 5-mm standard deviation added for the vertical components.The well-resolved resolutions near Shikoku and Kyushu using our inversion method are shown in Additional file 1: Fig. S5.Further, the time evolution was well reproduced in this study, as shown in Additional file 1: Fig. S6. (1)

Results
Figure 2 shows the detrended position time series for the selected stations (Additional file 1: Fig. S7 for longer periods).We observed an eastward displacement from 2020 at site 021013 on the Kii Peninsula (Fig. 2a).Sites 950,422 and 03P216, adjacent to the Kii channel, showed eastward displacements from 2019 to 2022 (Fig. 2b and  c).This could be attributed to the l-SSE in the Kii channel.Sites 021051 and 940,083 in Central Shikoku showed southeastward displacements of approximately 2 cm from 2019 to 2023, with uplifting of 2-4 cm, indicating an l-SSE in Central Shikoku (Fig. 2d, e).
Site 051142 showed a southeastward transient of up to 4 cm and uplifting up to 4 cm from late 2018 to August 2019 (Fig. 2f ).We observed a transient with an uplift from late 2018 to August 2019 at site 950437 (Fig. 2g).The transients at sites 051142 and 950437 were due to the Bungo-Channel l-SSE, which started in late 2018.In the southern Hyuga-nada, sites 021088 and 02P211 showed southeastward displacements from early 2018 to mid-2019, mid-2020 to 2021, and from January 2023 to July 2023, indicating the appearance of the l-SSE within these three periods (Fig. 2h, i).Site 021088 showed an uplift of up to 3 cm, corresponding to the southeastward transient (Fig. 2h).Site 960716, off the Ohsumi Peninsula, showed a southeastward displacement after the occurrence of the M w 5.7 earthquake on October 2, 2022 (Fig. 2j).Sites 960724 and 950492 in Tanegashima showed southeastward displacements after the 2019 Tanegashima earthquake on January 8, suggesting an afterslip of the main shock (Fig. 2k, l).Further, seemingly short-term eastward displacements were observed at these sites in mid-2020, mid-2021, mid-2022, and early 2023, indicating possible s-SSEs off Tanegashima during these periods (Fig. 2k, l).
Figure 3 shows the spatial pattern of the transient for the January 2018-July 2023 period.As shown in Fig. 3a, the Kii channel area showed an eastward transient of 1 cm, indicating the occurrence of l-SSE in the Kii channel.Further, there was a southeastward transient of up to 2 cm in central Shikou.This transient was distinct from the transient in the Bungo Channel, suggesting the occurrence of an l-SSE in this area.
We observed a southeastward horizontal transient up to 5 cm around the Bungo Channel and northern Hyuganada, with an uplift in southwestern Shikoku (Fig. 3b).The southeastward transient with subsidence in the northern Hyuga-nada and uplifting in the southwest Shikoku was caused by l-SSEs (Ozawa et al. 2020).Furthermore, stations in Miyazaki moved southward up to 5 cm and showed an uplift of up to 2 cm.This indicated the presence of l-SSE in the southern Hyuga-nada Sea.The sites in Tanegashima showed southeastward transients, indicating the afterslip of the January 2019 Tanegashima earthquake with M w 6.4.
Figure 4 shows the time evolution of the aseismic slip in the Nankai Trough from 2018 to 2023, and Fig. 5 shows the total aseismic slip between 2018 and 2023, with the time evolution of the moment for the designated areas.Additional file 1: Fig. S8 shows the transient crustal deformations corresponding to the periods shown in Fig. 4.
In 2018, an aseismic slip occurred in the northern Hyuga-nada Sea and Bungo Channel (Fig. 4a), and by 2019, the aseismic slip magnitude in the Bungo Channel increased as the slip area increased (Fig. 4b).However, in 2020, the slip area in the middle of the Bungo Channel abated.The maximum slip in this area was estimated to be approximately 30 cm (Fig. 5a).The estimated magnitude of the moment between April 2018 and August 2019 was 6.8, (Fig. 5e).As shown in Fig. 5e, the release in the Bungo Channel occurred gradually in 2018, and in 2019, the release rate increased; however, it decreased thereafter, consistent with the evolution of the slip shown in Fig. 4. The maximum slip and magnitude of the 2019 event were similar to those of previous major Bungo-Channel l-SSEs, as discussed by Ozawa et al. (2020).
The northern Hyuga-nada l-SSE started in mid-2018 and leveled off in late 2018, increased in magnitude in early 2019, and thereafter, increased gradually (Fig. 5f ).The gradual moment release from 2020 appeared to be due to the overcorrection of post-seismic deformation, which showed a westward trend.As discussed in the discussion section, we could not conclude that the northern Hyuga-nada l-SSE continued from 2020, as shown in Fig. 4c and d.Further, the northern Hyuganada slip in 2022 was due to the transient before and after the 2023 Hyuga-nada earthquake, which occurred within the subducting Philippine Sea Plate and may not reflect the true slip on the plate interface in this area.
The l-SSE in the southern Hyuga-nada Sea area, shown in Fig. 4, occurred from early 2018, preceding the May 2019 Hyuga-nada earthquake with M w 6.2 (Fig. 4a, b).This event was indicated by a southward transient at sites in the Miyazaki area (Fig. 3 and Additional file 1: Fig. S8).Further, the southern Hyuga-nada l-SSE occurred near the coastal area in 2018, and its slip area moved offshore in 2019 (Fig. 4a, b).We could not rule out the possibility that the removal of offsets from the 2019 Hyuga-nada earthquake was incomplete and, thus, affected the slip estimates for 2019.After the 2018-2019 event, the l-SSE started in mid-2020 (Fig. 4c) and continued into 2021 in the shallower part of the area (Fig. 4d), then from January 2023, the southern Hyuga-nada l-SSE started again (Fig. 4f ).The estimated moment magnitudes of the 2018-2019, 2020-2021, and January-July 2023 l-SSEs were Mw6.4,6.4, and 6.3, respectively.In 2020, The slip area covered an inland area and was larger compared to that of the 2023 event (Additional file 1: Fig. S9 for the shaded periods in Fig. 5g).The time evolution of the aseismic slip was similar to that of the moment (Fig. 5g).There was also an increase before the May 2019 earthquake, and a continued increase followed this.The release leveled off in 2019, increased from mid-2020, and leveled off again between mid-2021 and late 2022.After that, it began to increase into 2023 (Fig. 5g).
An afterslip was observed off Tanegashima after the 2019 earthquake (M w 6.4; Fig. 4b).The estimated moment of this afterslip was M w 6.7 (Fig. 5i).The time evolution of the release occurred logarithmically in 2019 and 2020.Further, the estimated aseismic slip showed short-term events in mid-2020, mid-2021, mid-2022, and late 2022.These events were observed in the position time series at the sites in Tanegashima (Fig. 2k, l) and were reflected in the time evolution of the moment (Fig. 5i).Short-term increases were observed in mid-2020, mid-2021, mid-2022, and early 2023, as shown in Fig. 5i.Although the 1-to 2-month durations of these increases were longer than the several-day durations typical of s-SSEs in the Shikoku area (e.g., Hirose and Obara 2010), in this study, we categorized these short-term events as s-SSEs.The slip areas of the possible s-SSEs corresponding to these short-term moment increases are shown in Additional file 1: Fig. S10, and the slip areas of the possible s-SSEs in 2020, 2021, and 2023 appeared to cover the offshore part of the afterslip area in 2019 as shown in Fig. 4b.
Central Shikoku exhibited an l-SSE in 2019, as shown in Fig. 4b-f.In 2022, we observed an increase in the slip amount and the size of the area (Fig. 4e).The estimated moment of this l-SSE was M w 6.5 (Fig. 5d).The time evolution of the moment indicated several coincidences between the start of the aseismic slip and the start of the s-SSE in northwest Shikoku.The Central Shikoku l-SSE is ongoing as of July 2023, with the moment release rate increasing from 2022.This increase corresponds to the slip increase in 2022, shown in Fig. 4e and f.
The Kii-Channel l-SSE appeared within the 2019-2021 period (Figs.4b-d).The magnitude of the Kii-Channel l-SSE increased, but showed fluctuation from early 2019 to late 2021, as the time evolution trend of the moment shown in Fig. 5c.However, the propagation of the slip area of the Kii-Channel l-SSE remains unclear.
In October 2022, a M w 5.7 interplate earthquake occurred off the Ohsumi Peninsula (Fig. 4e).At the time of and after the earthquake, an interplate aseismic slip was detected, and its release continued until early 2023 (Fig. 5h).Furthermore, an aseismic slip was observed before the 2022 earthquake.However, based on visual inspection of the position time series and goodness-offit to the data, this estimate before the 2022 earthquake was not as reliable as the aseismic slip observed after the 2022 earthquake.The inset of Fig. 5a shows the aseismic slip on the Kii Peninsula.Further, as shown in Fig. 5b, this area showed an increase from 2020, corresponding to M w 6.2 in 2023.

Discussion
The time-dependent inversion showed an aseismic interplate slip at several sites between 2018 and 2023 (Fig. 5a).In summary, we reported the l-SSEs detected in this study as well as those observed in previous studies (Fig. 6).
With regard to the annual components in position time series at several sites, the estimated moment evolution did not show clear annual changes.Thus, we reasoned that the effect of annual components at several sites had little effect on our results.
The main l-SSEs in the Bungo Channel detected via the GNSS network had an approximate recurrence interval of 5 to 6 years.With a slip deficit of 4 to 5 cm/year in the Bungo Channel area, the slip deficit was released near the area with the maximum slip.Correction of viscoelastic relaxation and s-SSE did not change the results of the slip estimates reported by Ozawa et al. (2020).However, owing to changes in grid spacing, our model showed that the northern Hyuga-nada and Bungo-Channel l-SSEs were separated from each other compared with those reported by Ozawa et al. (2020).
Northern Hyuga-nada l-SSEs occurred in 2002, 2005-2006, 2008-2009, 2014, 2016, and 2018(Takagi et al. 2019)).Given that the recurrence interval of these northern Hyuga-nada l-SSEs ranges from 2 to 5 years, the associated slip budget remains unclear.The obvious changes in the increasing trend in the moment observed in 2018-2019 indicated the occurrence of SSEs, while the changes observed within other periods (Additional file 1: Fig. S11) could not be defined (Figs. 5, Additional file 1: Fig. S11).Further, the position time series at site 021081 (Additional file 1: Fig. S12) showed a linear-like eastward movement when annual components were corrected, indicating that the northern Hyuga-nada l-SSEs from 2020 were probably due to the correction of viscoelastic relaxation.A realistic structural model in Finite Element Method may be necessary to solve this problem.Regardless, this remains to be further clarified in future studies.
The southern Hyuga-nada l-SSE has a recurrence interval of 2 years (Yarai and Ozawa 2013).Specifically, it occurred shortly before the 2016 Kumamoto earthquake, in 2018, 2020, and 2023.Therefore, it seemed its recurrence interval was not affected by the 2016 Kumamoto earthquake, which apparently affected the recurrence interval of the northern Hyuga-nada and Bungo-Channel l-SSEs (Ozawa et al. 2020).
We also observed a clear afterslip of the 2019 Tanegashima earthquake that continued until 2023, and its estimated moment was M w 6.7, greater than that of the main shock.Previous studies have shown that the interquartile range of the ratio of the moments of afterslips to coseismic slips varies in the range of 9-32% (Churchill et al. 2022).Thus, the afterslip of the Tanegashima earthquake released large amounts of energy, even though several other factors, such as detrending and s-SSEs, need to be considered before arriving at such a conclusion.It is also possible that the long-term trend was contaminated by transients, which resulted in an increase in the magnitude of the afterslip.Further, in the Tanegashima area, we observed a seemingly short-term transient in this region.This event occurred four times after the 2019 earthquake.The total moment of this afterslip was greater than that of the main seismic event, even after excluding s-SSEs.The time evolution of the moment release of the afterslip followed a logarithmic pattern in 2019 and 2020 and could be modeled by the logarithmic decay expression shown in Eq. ( 2) (Marone et al. 1991).
(See figure on next page.)Fig. 5 a Estimated total slip for the 2018-2023 period.The colors indicate the magnitude of the aseismic slip, whereas the arrows show the aseismic slip vector.The contour interval is 10 cm; that in the inset is 5 cm.The descriptions of the other symbols are the same as shown in Fig. 4. b-i Time evolution of moment release.Each area of moment calculation is designated by solid lines in (a).b Kii Peninsula of (1).c Kii channel of (2).d central Shikoku of (3).The vertical broken lines indicate the dates of the s-SSEs in northwest Shikoku.e Bungo of (4).f Northern Hyuga-nada of (5) g Southern Hyuga-nada of (6).The vertical line indicates the date of the 2019 Hyuga-nada earthquake.The slip distributions for the shaded periods are shown in Additional file 1: Fig. S9.h Off the Ohsumi Peninsula of (7).The vertical line indicates the date of the 2022 off-Ohsumi earthquake.i Off Tanegashima of (8).The vertical line indicates the date of the 2019 Tanegashima earthquake.The green translucent box shows the period of short-term slow-slip events off Tanegashima where M(t) represents the moment of the afterslip, μ represents its rigidity, S represents the area of the fault slip, t represents time, t 0 represents the time of the mainshock, and T r represents the relaxation time.Further, A is a constant related to stiffness, frictional parameters, and normal stress.We estimated μSA and T r (57 × 10 17 Nm and 0.3 years, respectively) by regressing Eq. ( 2) on the moment evolution data (Additional file 1: Fig. S13a).When we did not include possible s-SSE periods in the regression, the values of μSA and T r were 38 × 10 17 Nm and 0.14 years, respectively (Additional file 1: Fig. S13b).In both cases, the theoretical model showed good agreement between the observation and computation for the analytical period.
The Central Shikoku l-SSE occurred independently or was associated with the western Shikoku l-SSE in 1978, 2010, 2012, and 2015(Kobayashi 2012;Takagi et al. 2019).Further, its source area in 2019-2023 was found to be similar to that of past events, and its estimated moment was M w 6.5, consistent with that of the 2015 event that occurred independently.As shown in Fig. 5d, the aseismic slip often coincided with or began with the occurrence of s-SSEs in northwest Shikoku.Without correction for s-SSEs in northwest Shikoku, the Central Shikoku SSE showed an increase in the moment at the time of occurrence of the northwestern Shikoku s-SSEs probably owing to an overlap of s-SSE and l-SSE areas due to spatial smoothing.Thus, we cannot rule out the possibility of the incomplete removal of s-SSE offsets.
Although position time series at sites in Central Shikoku indicated time-varying evolution of SSE (Fig. 2d and e and Additional file 1: Fig. S14), we regarded the Central Shikoku SSE in this study period as one continuous event.
Kii-channel l-SSEs occurred in 1996, 2000-2004, and 2014-2016 with irregular recurrence intervals.The Kii-Channel l-SSE started in 2019 and continued in 2022, with fluctuations.Further, the magnitude of the Kii-Channel l-SSE increased in 2021, and its estimated moment was M w 6.3.The 2019-2022 events lasted longer than the 2014-2016 and 1996 events (Kobayashi 2014(Kobayashi , 2017)).Long-term fluctuations in the release moment of the 2019-2022 l-SSEs were possibly affected by the complicated geometry of the plate interface in this area.
After the M w 5.7 interplate earthquake off the Ohsumi Peninsula on October 2, 2022, we detected an interplate aseismic slip, which continued in 2023 with M w 6.1 near the epicenter area (Figs. 4e,f,5h).During the mainshock (M w 5.7), we observed small offsets, and GNSS sites gradually moved eastward after the small offsets of the (2) M(t) = µSAlog( t − t 0 T r + 1), event (Fig. 2j).The relationship between the main shock and the gradual transient remains to be examined in future studies.Aseismic slip has occurred at least once in this region (1998)(1999), as pointed out by Takagi et al. (2019).
We estimated the aseismic slip in the Kii Peninsula based on the displacements observed at several stations (Additional file 1: Figs.S15 and S16).The three stations moved eastward from 2020 to 2023.We also observed good agreement between the observations and calculations, except for the vertical displacements at site 970826.Therefore, owing to this inconsistency, we could not clearly state whether an aseismic slip occurred on the Kii Peninsula.
Most of the l-SSEs detected in this study were distinct from s-SSEs and locked areas, consistent with previously reported findings (e.g., Obara and Kato 2016).The southern Hyuga-nada l-SSE in this study seemed to coexist with the afterslip of the 2019 earthquake in regions at depths deeper than that of the hypocenter (Additional file 1: Fig. S9).Via simulation based on the rate and statedependent friction law (Dieterich 1979) in the velocityweakening region (A-B < 0) with cutoff velocity, Nakata et al. (2012) reproduced the coexistence of an afterslip and SSEs in Hyuga-nada.Their results were inconsistent with the observations made in several aspects, e.g., the observed recurrence intervals of the SSEs were inconsistent.Thus, modifications are required to interpret the cases observed in this study.
The Tanegashima afterslip was followed by s-SSEs that occurred in almost the same areas at depths shallower than or equal to that of the hypocenter (Additional file 1: Fig. S10).The moments of the afterslip and s-SSEs in the Tanegashima area follow a logarithmic decay trend, indicative of the stress relaxation process in the velocity-strengthening region (A-B > 0) (Marone et al. 1991).However, a velocity-strengthening region with a onestate variable in the rate-and state-dependent friction law cannot cause an SSE (Nakata et al. 2012).Therefore, further studies are required to examine whether A-B is positive or negative in the Tanegaeshima afterslip area.The cases of southern Hyuga-nada and Tanegashima possibly suggest the existence of heterogeneity in frictional properties in short ranges in Kyushu (Additional file 2: Tables S1, S2).
The total slip in the Nankai Trough showed a spatially sporadic pattern between 2018 and 2023 (Fig. 5a).Thus, we calculated the change in Coulomb failure stress (ΔCFS) on the plate interface with slip direction opposite to that of the motion of the Philippine Sea Plate, with friction coefficient of 0.2.Consequently, we observed positive ΔCFS values for the areas neighboring the spatially sporadic slip zones, which increased the possibility of SSE occurrence in these areas (Additional file 1: Figs.S17 and S18).

Conclusions
We estimated aseismic interplate slippage on the plate interface between the Philippine Sea Plate and Amur Plate by applying time-dependent inversion to the GNSS position time series.The results obtained showed the occurrence of several independent l-SSEs along the Nankai Trough subduction zone between 2018 and 2023.However, it remains unclear whether these spatially isolated l-SSEs are coincidentally or physically related.Spatially, the Nankai Trough area slipped sparsely between 2018 and 2023, resulting in an increase in ΔCFS in the neighboring areas.Further, we noted that in southern Hyuga-nada, an aseismic slip preceded the M w 6.2 Hyuga-nada earthquake in 2019, indicating the potential of monitoring interplate slip with respect to the estimation of seismic activity.Our results contribute to enhancing understanding regarding the subduction processes in the Nankai Trough subduction zone.

Fig. 1 a
Fig. 1 a Tectonic setting of four converging plates around Japan.b Enlarged map of the rectangular area in (a).The broken contours represent the iso-depth contours of the plate interface between the subducting Philippine Sea Plate and the overriding Amur Plate.The contour interval is 20 km.The red dots indicate the epicenters of low-frequency earthquakes, and the numbered ellipses schematically indicate long-term slow-slip events (l-SSEs) along the Shikoku and Kyushu coasts.(1) Kii channel, (2) Eastern Shikoku, (3) Central Shikoku, (4) Western Shikoku SSE, (5) Bungo SSE, (6) Northern Hyuga-nada SSE, and (7) Southern Hyuga-nada SSE.The stars show the epicenters of the 2019 Hyuga-nada Sea, 2019 Tanegashima, and 2022 Off Ohsumi earthquakes with M w 6.2, 6.4, and 5.7, respectively.c The 1946 Nankai area showing the source region of the 1946 Nankai earthquake(Sagiya and Thatcher 1999), which is locked strongly and expected to slip in the future Nankai earthquake.The green line indicates the area of the plate interface adopted in the time-dependent inversion.The blue circles show the locations of the GNSS sites used in the time-dependent inversion, whereas the white circles show the locations of the GNSS sites whose position time series are shown in Fig.2.The T-direction and S-direction represent the directions toward which the parameters of the parametric spline surface increase(Ozawa et al. 2020)

Fig. 2
Fig. 2Detrended position time series at selected Global Navigation Satellite System (GNSS) sites relative to the 950388 site (the locations of these sites are shown in Fig.1a and c) averaged over a 3-day period.EW, NS, and UD indicate east-west, north-south, and up-down components with eastward, northward, and upward positive directions, respectively.The red lines show the transient computed from our optimal model.The vertical broken line in (j) indicates the date of the 2022 Ohsumi earthquake.The blue translucent boxes show the periods of long-term slow-slip events (l-SSEs), while the green translucent box in (k) and (l) shows the periods of short-term slow-slip events (s-SSEs) off Tanegashima.a 021013, b 950,422, c 03P216, d 021051, e 940,083, f 051142, g 950,437, h 021088, i 02P211, j 960,716, k 960,724, l 950,492 (See figure on next page.)

Fig. 3
Fig. 3 Observation (black arrows) and optimal model-based computed (white arrows) detrended crustal deformation.The error ellipses indicate 1σ error.a Horizontal displacements within the January 2018-July 2023 period.b Vertical displacements within the January 2018-July 2023 period

Fig. 4
Fig. 4 Estimated interplate slip on the plate interface.The colors are indicative of the magnitude of the aseismic slip.The contour interval is 4 cm.The solid arrows represent aseismic slips larger than 3σ, whereas the gray arrows represent slips less than 3σ.The broken contours indicate the iso-depth contours of the plate interface between the subducting Philippine Sea plate and the overriding Amur plate with a contour interval of 20 km.The blue dots show the epicenters of low-frequency earthquakes for the corresponding period.The stars show the epicenters of the M w 6.2 Hyuga-nada earthquake on May 10, 2019, and the M w 6.4 Tanegashima earthquake on January 8, 2019, in (b). a January 1, 2018-January 1, 2019.b January 1, 2019-January 1, 2020.c January 1, 2020-January 1, 2021.d January 1, 2021-January 1, 2022.e January 1, 2022-January 1, 2023.The star indicates the epicenter of the M w 5.7 earthquake of October 10, 2022.f January 1, 2023-July 1, 2023 (See figure on next page.)