Long-term slow slip events along the Nankai trough delayed by the 2016 Kumamoto earthquake, Japan

The Global Navigation Satellite System network in Japan detected transient crustal deformation along the Nankai trough, Japan, from June 2018. Time-dependent inversion analysis shows that a long-term slow slip event in northern Hyuga-nada Sea along the Nankai trough, Japan, started in June 2018 and decayed in October 2018. From October 2018, a slip area appeared in the Bungo channel and expanded to northern Hyuga-nada Sea and southwest Shikoku at the time of the maximum slip rate. The slip area in the middle of the Bungo channel started to abate around May 2019, with the slip in neighboring areas nearly stopping in August 2019. The estimated rupture propagation is different from those in the past Bungo channel SSEs, in which rupture propagated southwestward from the Shikoku side to the Kyushu Island side at the time of the maximum slip rate. Furthermore, the rupture in northern Hyuga-nada Sea preceded the Bungo channel SSE and reappeared together with the Bungo channel SSE at the time of the 2018–2019 event, though the northern Hyuga-nada Sea SSE followed the 2009–2011 Bungo channel SSE. There is a possibility that the differences in the rupture propagation and recurrence interval from the past events are due to the 2016 Kumamoto earthquake. The adjacent locked area along the Nankai trough subduction zone is a well-known seismic gap and the 2018–2019 SSE changed the stress state in favor of the occurrence of nearby subduction earthquakes.


Introduction
Geodetic networks such as the Global Navigation Satellite System (GNSS) have detected many slow slip events (SSEs) worldwide (e.g., Schwartz and Rokosky 2007). SSEs usually occur in the transient area between a locked area and an aseismically creeping area (e.g., Schwartz and Rokosky 2007). Many SSEs have approximate recurrence intervals such as ~ 6 years for the Bungo channel longterm SSE (Yoshioka et al. 2015), ~ 14 months for the Cascadia SSE (Rogers and Dragert 2003), 3-4.5 years for the Guerrero SSE, Mexico (Graham et al. 2016), ~ 20 months for the Costa Rica SSEs (Voss et al. 2017), 2-5 years for the Hikurangi SSEs, New Zealand (Wallace and Beaven 2010), and 8 years for the Alaska SSE (Fu and Freymueller 2013). SSEs have a variety of durations. We call SSEs with durations over several months long-term SSEs and those with durations of several days short-term SSEs in this study. The mechanism behind the observed variation of SSEs still remains to be clarified. Furthermore, the relationship between large earthquakes and SSEs is unclear. Since some SSEs have triggered nearby large earthquakes (e.g., Ruiz et al. 2017), monitoring of SSEs is indispensable to assess the hazard potential of nearby seismogenic areas.
The area along Shikoku and Kyushu Islands, Japan, in the Nankai trough subduction zone has experienced long-term SSEs many times and is a good place to investigate the characteristics of long-term SSEs (e.g., Takagi et al. 2019). Furthermore, the Nankai trough subduction zone is a well-known seismic gap (e.g., Kumagai 1996), and it is important to investigate SSEs along the Nankai Open Access *Correspondence: ozawa-s96sa@mlit.go.jp Geospatial Information Authority of Japan, Ibaraki 305-0811, Japan trough to estimate the slip budget and properties of the plate interface. Figure 1 shows the tectonic setting in Kyushu and Shikoku along the Nankai trough, southwest Japan. The Philippine Sea plate subducts beneath the Amur plate from the Nankai trough at a rate of 6 to 7 cm/year (Kreemer et al. 2014). Because of the subduction of the Philippine Sea plate, large interplate earthquakes have repeatedly occurred off Shikoku with a recurrence interval of around 140 years (e.g., Kumagai 1996) and in Hyuganada Sea with a recurrence interval of 20-30 years (e.g., Shiono et al. 1980). In down-dip areas of a seismogenic zone, long-term SSEs repeatedly occur, with short-term slow slip events located in the down-dip extension of the long-term slow slip area (e.g., Obara 2002;Hirose and Obara 2005). Long-term slow slip areas are approximately segmented into western Shikoku, Bungo channel, northern Hyuga-nada Sea, and southern Hyuga-nada Sea in this region (see Fig. 1b). Western Shikoku SSEs have a recurrence interval of around 6 years and a duration of several months to 1 year (Takagi et al. 2016(Takagi et al. , 2019. Major Bungo channel SSEs occurred in 1996SSEs occurred in -1998SSEs occurred in , 2003SSEs occurred in -2004SSEs occurred in , and 2009SSEs occurred in -2011 with a duration of 6 months to 2 years (Hirose et al. 1999;Ozawa et al. 2001;Yoshioka et al. 2015;Takagi et al. 2019). Southern Hyuga-nada Sea SSEs and northern Hyuga-nada Sea SSEs have occurred independently or simultaneously (e.g., Ozawa 2017;Takagi et al. 2019). Northern Hyuga-nada Sea SSEs occurred in 1996-1997, 2002, 2008, 2014 with a duration of around 6 months (Ozawa 2017;Takagi et al. 2019 Southern Hyuga-nada Sea SSEs occur with a recurrence interval of 2 years and a duration of up to 1 year (e.g., Yarai and Ozawa 2013). Before the 2016 Kumamoto earthquake sequence, long-term slow slip events were occurring in the Hyuganada Sea and southwestern Shikoku (Ozawa 2017) (Additional file 1: Figure S1). Since the Kumamoto earthquake, it has been unclear whether these SSEs have ceased or continued because of the large coseismic and postseismic crustal deformation due to the Kumamoto earthquake.
Two years after the Kumamoto earthquake, the GNSS network in Japan (GEONET) detected a transient from June 2018 along the Nankai trough subduction zone. The spatial pattern of this transient suggests the occurrence of a long-term SSE in the Nankai trough subduction zone. In this study, we investigate the spatio-temporal evolution of aseismic slip along the Nankai trough subduction zone from 2018 to 2019 on the basis of GNSS data. We compare the 2018-2019 event with the past major longterm SSEs in this area. We also discuss the possibility of stress perturbation caused by the 2016 Kumamoto earthquake affecting the long-term SSEs in 2018-2019.

Analytical method
GNSS data were analyzed using Bernese GPS software (version 5.0) on a daily basis with precise ephemeris and earth rotation parameters (Hatanaka et al. 2003). We transformed coordinates of longitudes, latitudes, and heights to east-west, north-south, and up-down components in a local framework, respectively. After this transformation, we removed annual components estimated for the period between 2000 and 2018 by regressing trigonometric and polynomial functions. We removed the linear trend for the period between January 1, 2017 and January 1, 2018 by regressing a linear function from the position time series without annual components.
The position time series show that the effect of the 2016 Kumamoto earthquake seems to have become steady after January 2017. Thus, we think that the above detrending is not significantly affected by the postseismic deformation following the Kumamoto earthquake. The adopted steady crustal deformation rate from January 2017 to January 2018 is shown in Additional file 1: Figure  S2. In addition to the 2018-2019 SSE, we also analyzed the past major Bungo channel SSEs. In the past events, we adopted the linear trend from October 2007 to May 2009 in the detrending of position time series, when there was no transient.
We applied time-dependent inversion (McGuire and Segall 2003;Ozawa 2017) to the east-west, northsouth, and up-down components at 240 selected GPS sites (see Fig. 1b) for the period between January 1, 2017 and August, 1, 2019. We averaged the data over 3-day windows and used the data every 3 days in timedependent inversion. We weighted horizontal and vertical components at a ratio of 5:1 in the following analysis, considering repeatability.
In the time-dependent inversion, we used a parametric spline surface to represent a fault interface as described in Ozawa et al. (2001). We 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. 1b. The parametric spline surface is represented by two parameters, ξ and η. The direction of increasing ξ is named T-direction, and the direction of increasing η is named S-direction as shown in Fig. 1b. In our representation of the spline surface, we took 35 grid knots in the T-direction and 15 grid knots in the S-direction in Fig. 1b. The spacing of the grid knot in the T-direction and S-direction is approximately 10 km and 20 km, respectively. Green's function was calculated following Ozawa et al. (2001).
We adopted the following state x(t) in the timedependent inversion: where u k indicates the slip and v k indicates the slip velocity at the kth grid knot on the plate interface, which is represented by the spline surface of a total of L grid knots. e EW , e NS , and e UD indicate common mode errors in the position time series (Fukuda 2018). We incorporated spatial roughness in the transition equation as follows (Ozawa 2017): where M is the roughness matrix described in Ozawa et al. (2001) and λ is the hyperparameter of spatial smoothing. u n|n−1 and u n−1|n−1 are the predicted slip at time n based on the data until time n − 1 and the estimated slip at time n − 1 based on the data until time n − 1, respectively. Δt is the time lapse. v n|n−1 is the predicted slip velocity at time n based on the data until time n − 1. e 1 and e 2 are Gaussian noise with 0 mean and a standard deviation of 1.
We constrained the aseismic slip vector and slip velocity vector on the plate surface to within 10°from the direction opposite to the motion of the Philippine Sea plate relative to the Amur plate using the method of Simon and Simon (2006). The direction of the Philippine sea plate relative to the Amur plate was taken from Kreemer et al. (2014). (1) We conducted a checkerboard test by inverting the synthesized data with Gaussian errors with 1 mm standard deviation added for horizontal components and 5 mm standard deviation added for vertical components. The resolution near Shikoku and Kyushu in our inversion method is shown in Additional file 1: Figure S3, and the resolution is well resolved near Shikoku and Kyushu. The time-dependent inversion shows the aseismic interplate slip between June and August 2018 in northern Hyuga-nada Sea (Fig. 4b). This area was experiencing long-term slow slip shortly before the 2016 Kumamoto earthquake (Additional file 1: Figure S1). No low-frequency earthquakes were observed in the area neighboring the estimated slow slip area. The aseismic slip area in northern Hyuga-nada Sea ended in October 2018 and the slip area moved to the Bungo channel from October 2018 (Fig. 4c, d). During this period, the number of low-frequency earthquakes increased in the down-dip extension of the long-term slow slip area in the Bungo channel. The magnitude of aseismic slip increased between December 2018 and February 2019. Low-frequency earthquakes continued to occur in the Bungo channel. The aseismic slip area propagated to Shikoku and to northern Hyuganada Sea with the slip magnitude increasing between February and April 2019. The slip rate is maximum for this period. The slip area in the middle of the Bungo channel abated between April and June 2019. The neighboring area decayed between June and August 2019. For this period, the number of low-frequency earthquakes decreased. The estimated moment magnitude between April 2018 and August 2019 was 7.0 with a rigidity of 30 GPa. The maximum slip was estimated to be approximately 30 cm in the Bungo channel area and 36 cm in northern Hyuga-nada Sea (see Table 1 and Additional file 1: Figure S4). The estimated magnitude and maximum slip in 2018-2019 are near the values for the past major Bungo channel SSEs (see Table 1).

Results and discussion
In the last major Bungo channel long-term SSEs, the slip area moved southwestward at the time of the maximum slip rate (Yoshioka et al. 2015). Our inversion also shows southwestward rupture propagation at the time of the maximum slip rate for the past three events (Additional file 1: Figures S5-S7). Thus, the 2018 event shows different slip propagation in that the slip area expanded to southwest Shikoku and to northern Hyuga-nada Sea. Long-term SSEs were occurring in northern Hyuganada Sea and southwestern Shikoku before the 2016 Kumamoto earthquake (Ozawa 2017). Additional file 1: Figure S8 shows ΔCFS on the plate interface from the 2016 Kumamoto earthquake with a friction coefficient of 0.2 and slip direction opposite to the motion of the Philippine Sea plate. The fault parameters of the foreshock and main shocks are from Kobayashi et al. (2018). We can see a negative ΔCFS in northern Hyuga-nada Sea and the Bungo channel up to around 10 kPa. Takagi et al. (2019) estimated the average recurrence interval of the northern Hyuga-nada Sea SSEs to be 6 years. The latest events occurred in 2014 and 2016 separated by 2 years. Considering the average recurrence interval of 6 years, we consider that three events occurring in a series with a 2 year recurrence interval has a low probability. Furthermore, the northern Hyuga-nada Sea SSE occurred twice during the 2018-2019 event. In addition, the estimated maximum slip of 36 cm in northern Hyuga-nada Sea is too large for the plate convergence to compensate, considering the 4 to 5 cm/year slip deficit rate in this area. This may indicate that the SSE in northern Hyuga-nada Sea which was occurring shortly before the Kumamoto earthquake did not release energy fully. We cannot rule out the possibility that the estimated sequence of the northern Hyuga-nada Sea SSE has been affected by the stress perturbation of the 2016 Kumamoto earthquake. Considering the negative ΔCFS in northern Hyuga-nada Sea, we hypothesize This hypothesis remains to be examined by observing the long-term SSE sequence in this area in the future. Based on the roughly estimated back slip model using two rectangular faults, ΔCFS due to the plate convergence is estimated to be less than 100 kPa in the maximum slip area in northern Hyuga-nada Sea and Bungo channel for 1 year (see Additional file 1: Figure S9). Stress drop of the SSE is roughly 100 kPa assuming a typical length of 20 km and 30 cm slip with a rigidity of 30 GPa (Eshelby 1957). The stress budget among the plate convergence, the Kumamoto earthquake, and stress drop of the SSE remain to be solved in a future study much more quantitatively in order to examine our hypothesis. The major Bungo channel SSEs detected by the GNSS network have an approximate recurrence interval of 5 to 6 years. The recurrence interval between the 2009-2011 event and the 2018-2019 event was approximately 7 years. Thus, there is a possibility that the 2018 event was delayed In fact, we expected a Bungo channel SSE in 2016, when southwest Shikoku started undergoing aseismic slip. We cannot rule out the possibility that this delay of 1 to 2 years may have been caused by the stress perturbation of the 2016 Kumamoto earthquake, which shows a negative ΔCFS in the Bungo channel area.
Although no low-frequency earthquakes occurred at the time of the northern Hyuga-nada Sea SSE in 2018, there was an increase in the number of low-frequency earthquakes associated with the long-term Bungo channel SSE in the down-dip extension of the long-term slow slip area (see Fig. 4). This activation of low-frequency earthquakes associated with long-term slow slip events in the Bungo channel has also been observed in the previous Bungo channel long-term SSEs (e.g., Hirose and Obara 2005) (Additional file 1: Figure S10). We believe that this association is due to stress changes caused by long-term slow slip events, which show positive ΔCFS in a low-frequency earthquake area (Additional file 1: Figure S11). Since the short-term SSEs that occurred in the study period accompanying low-frequency earthquakes were very small in magnitude, they do not affect the results of our analysis.
With regard to the slip budget on the plate interface, the slip deficit rate in the SSE area is around 4 to 5 cm/ year (e.g., Yoshioka et al. 2015). The estimated aseismic   Table 1 and Additional file 1: Figure S4). Thus, the slip deficit is being mostly released near the area with maximum slip, which is the western part of the slip area, although the eastern part of the slip area does not release all the accumulated slip budget of the slip deficit. This result is consistent with the conclusion by Yoshioka et al. (2015). Northern Hyuga-nada Sea SSEs occurred in 2002-2006, 2010-2011, 2014(Ozawa 2017Takagi et al. 2019). Since the recurrence interval of northern Hyuganada Sea SSEs ranges from 2 to 5 years, the slip budget in the northern Hyuga-nada Sea SSEs remains unclear.
The effect of the 2018-2019 long-term SSE on the seismic gap shows an increase in ΔCFS of up to 10 kPa in Hyuga-nada Sea and offshore of Shikoku (Additional file 1: Figure S11). Thus, the 2018-2019 SSE increased the potential of the expected Nankai earthquake, whose source area is schematically shown in Fig. 1b. Based on a simple back slip model, the ΔCFS in the seismic gap is roughly 10 kPa for 1 year and comparable to the effect of the SSE, although more detailed back slip modeling is necessary for further discussion (Additional file 1: Figure  S9).
With regard to the Hyuga-nada area, a medium-size earthquake with Mw 6.2 occurred on May 10, 2019 in southern Hyuga-nada Sea. Figure 4d, e shows the aseismic slip in southern Hyuga-nada Sea in the down-dip extension of the 2019 Hyuga-nada earthquake before the earthquake. Although the magnitude of the estimated aseismic slip in southern Hyuga-nada Sea is very small, we cannot rule out the possibility that the aseismic slip in the down-dip extension preceded the May 10, 2019 earthquake. This hypothesis remains to be investigated in a future study.

Conclusion
Because the 2016 Kumamoto earthquake with Mw 7.6 changed the stress state on a plate interface, the 2018-2019 event is a good example of a long-term SSE affected by the stress perturbation of a large earthquake. We propose a hypothesis that the northern Hyuga-nada Sea SSE shortly before the 2016 Kumamoto earthquake ceased at the time of the earthquake and restarted 2 years later. We also suggest the possibility that the Bungo channel long-term SSE was delayed by the stress perturbation of the 2016 Kumamoto earthquake. The entire rupture propagation, which is different from those of the past Bungo channel area SSEs, may have also been caused by the perturbation of the 2016 Kumamoto earthquake. Further examination of our hypothesis is very important to understand the factors controlling the long-term SSE rupture sequence. Since a long-term SSE affects the stress state in a seismogenic zone, it is of great importance to monitor long-term SSEs. In fact, the 2018-2019 longterm SSE changed the stress state in favor of the expected Nankai earthquake at the bottom of the locked seismogenic zone. We propose the possibility that aseismic slip preceded the Mw 6.2 Hyuga-nada earthquake, which indicates the potential of assessing earthquake occurrence.