Detection of a recent large Hyuga-nada long-term slow slip event and estimation of its spatiotemporal slip distributions

Using horizontal and vertical GNSS time series data from the GSI in sotheastern Kyushu from January 1, 2017 to June 30, 2022, we detected a recent long‑term slow slip event (L‑SSE) that occurred in the Hyuga‑nada region, southwest Japan, and estimated its spatiotemporal slip distribution. We performed such analysis considering the piecewise evolution of slip over time with time windows of 0.8 years to ensure a good signal‑to‑noise ratio in the horizontal displacements in each time window. The results showed that a slip of more than 2 cm occurred in the central part of Miyazaki Prefecture from 2018.5 to 2019.3 (the units were 0.1 years = 36.5 days). Then, the amount of slip increased and expanded slightly to the southern part of Miyazaki Prefecture from 2019.3 to 2020.1, and the amount of slip reached a maximum of 3.9 cm in the subsequent period (2020.1–2020.9). A smaller slip occurred at almost the same location in the following 0.8 years. Therefore, the duration of the L‑SSE was approximately 3.2 years from approxi‑ mately 2018.5 to 2021.7. The annual average maximum slip rate was approximately 4.9 cm/yr (3.9 cm/0.8 yr) during the period from 2020.1 to 2020.9. The maximum total slip was estimated to be approximately 12.9 cm, and the equivalent release moment was 4.9 × 10 19 Nm, corresponding to Mw7.1. Compared to previous L‑SSEs that occurred in the Hyuga‑nada region, the annual average maximum slip rate was relatively low, and the main total slip was estimated at almost the same location on the plate interface. On the other hand, the slip duration of this L‑SSE was the longest, and the release moment and moment magnitude were the greatest. The total slip area of more than approximately 10 cm in the Hyuga‑nada L‑SSE estimated in this study almost overlapped with the afterslip area of the December 3, 1996 Hyuga‑nada earthquake, with a depth range of approximately 30–40 km. Short‑term SSEs also occurred in the same depth range. These coincidental same depth ranges of interplate seismic events stem from large thermal gradients. The large thermal gradients play an important role in narrowing the depth range of frictional parameters ( a–b ) at the plate interface. In addition, the pore pressure and normal stress are closely related to the critical stiffness. In particular, contrasting values of pore pressure and normal stress at the shallow and deep sides of the mantle wedge corner may also contribute to transient aseismic slip within the same depth range.


Introduction
In the Nankai Trough subduction zone in southwest Japan, the Philippine Sea (PHS) plate, a young oceanic plate with an age of 26-15 Ma (e.g., Okino et al. 1994), is subducting at a low angle beneath the continental Amurian (AM) plate.In the Nankai Trough, M8-class megathrust earthquakes have occurred at intervals of approximately 90-150 years (e.g., Ando 1974;Kumagai 1996), and many researchers around the world have conducted research on seismic events that occur at this plate boundary.Since the discovery of tectonic tremors in the Nankai Trough subduction zone by Obara (2002) around the turn of this century, studies on slow earthquakes have been actively conducted.Megathrust earthquakes and slow earthquakes are slip phenomena that occur on plate interfaces.According to Obara and Kato (2016), shallow very-low-frequency (VLF) and tremor events, megathrust earthquakes, long-term slow slip events (hereafter referred to as L-SSEs), and deep episodic tremors and slips (ETSs) occur from shallower to deeper nonoverlapping depth segments across the Nankai Trough plate interface.
The Hyuga-nada region, the target area of this study, is located to the west of the above-mentioned subduction zone, is also an area of subduction beneath the Kyushu region, and has a significantly different tectonic environment.In this region, the PHS plate is subducting beneath the continental AM plate in the northwestward direction at convergence rates of 6.2-6.4 cm/yr (DeMets et al. 2010) (Fig. 1).Here, the Kyushu-Palau Ridge, a remnant arc associated with back-arc spreading of the Shikoku Basin, is being subducted.On its west side, the old (> 50 Ma) PHS plate is subducting at a high angle, and no M8-class megathrust earthquakes have been identified for at least the past 400 years (Earthquake Research Committee 2022).
L-SSEs with durations ranging from months to years have occurred on plate boundaries in and around the Japanese Islands.The signatures of such aseismic slow slip phenomena have been recorded in Japan by the GNSS Earth Observation Network System (GEONET), which consists of approximately 1300 GNSS sites operated by the Geospatial Information Authority of Japan (GSI) since 1996.GEONET plays an important role in monitoring crustal deformations caused mainly by earthquakes and volcanic activity.Specifically, L-SSEs have also been identified via GNSS time series data from previous studies in the Hyuga-nada region (Yarai and Ozawa 2013;Ozawa 2017;Takagi et al. 2019) (Table 1, Fig. 8).
Regarding the depth distribution of slow earthquakes at plate boundaries in the Hyuga-nada region, VLFs and tectonic tremors occur in and around shallow plate interface (e.g., Yamashita et al. 2015Yamashita et al. , 2021)), but deep tectonic tremors rarely occur.Interestingly, L-SSEs and shortterm SSEs occur at almost the same depth range (Okada et al. 2022).Furthermore, the afterslip areas of the two M6-class Hyuga-nada earthquakes that occurred in October and December 1996 (e.g., Yagi et al. 2001) also overlap with the areas where these SSEs occurred (e.g., Yarai and Ozawa 2013).Therefore, the depth distributions of the physical, frictional, and mechanical properties at the plate boundary may be quite different from those of the adjacent Nankai Trough to the east (e.g., Gao and Wang 2017), but it remains unclear what controls this depth distribution and the along-arc variation.
In this study, we analyze GNSS time series data from January 1, 2017 to June 30, 2022 and detect an L-SSE in the Hyuga-nada region that has not been previously documented.We infer the spatiotemporal slip evolution associated with such a recent L-SSE by using our own inversion method (Yoshioka et al. 2015;Seshimo and Yoshioka 2022).Based on these results, we clarify similarities or differences between these events and past L-SSEs.Furthermore, by comparing our results with those of previous studies, we attempt to deepen the understanding of the depth distribution of frictional properties at the plate boundary in the Hyuga-nada region.

Data
In this study, we used the daily coordinate F5 solution at 26 GNSS sites (Fig. 1) resulting from the final analysis of the GSI GEONET observation data (Takamatsu et al. 2023).The three-component time series spans from January 1, 2017 to June 30, 2022.We selected the GNSS stations via the following method: since we focused only on detecting an L-SSE in the Hyuga-nada region and estimating its spatiotemporal slip evolution, avoiding GNSS data related to the 2018-2019 Bungo Channel L-SSE (e.g., Ozawa et al. 2020;Seshimo and Yoshioka 2022), we did not use GNSS stations in the northern half of the Kyushu region (see Additional file 1: Fig. S1(a)).We selected all the stations in the southern half of the Kyushu region.However, in the central part of the selected stations, postseismic deformations associated with the Mw7.0 Kumamoto earthquake on April 14, 2016 were identified in the period of analysis, showing slight northward displacements for the processed time series data (see Additional file 1: Fig. S1(b)), which are not related to the Hyuga-nada L-SSE.Similarly, the data from the GNSS stations in and around Sakurajima volcano, which erupted frequently throughout the year, were contaminated by volcanic deformation processes.Additionally, despite the consistent horizontal displacement directions between the GNSS stations in western Kagoshima Prefecture and those in eastern Kagoshima Prefecture, the amounts of former displacements are greater than those in central Kagoshima Prefecture.Such a horizontal displacement field is not consistent with a slip distribution attributable to the L-SSE in the Hyuga-nada region and may be strongly influenced by other unknown mechanisms.Therefore, we removed the data from such stations from the dataset (see gray displacements in Fig. 3a, b) and used the data from the  Baba et al. (2002), the 20-50 km isobath by Hirose et al. (2008), and the isodepth contours deeper than 60 km by Nakajima and Hasegawa (2007).The light green circle represents the approximate location of the Bungo Channel.The green bold line indicates the outer edge of the subducted Kyushu-Palau Ridge (KPR) by Yamamoto et al. (2013).The red and blue solid circles indicate the GNSS observation stations analyzed in this study whose detailed information is listed in Additional file 1: Table S1, and the blue solid squares represent the reference stations used in this study.The yellow-green box represents the horizontal projection of the source region of the assumed model on the 3D plate boundary.The red solid star represents the epicenter of the Mw7.0 Kumamoto earthquake on April 14, 2016.Processed time series data for all the stations (red and blue solid circles) with station codes are shown in Fig. 2. The gray boxed area in the inset shows our study area remaining 26 GNSS stations shown in Fig. 1 (see also Additional file 1: Table S1).

GNSS time series data analysis
As in Yoshioka et al. (2015) and Seshimo and Yoshioka (2022), we selected six reference stations in the northern part of the Chugoku region (940075, 940076, 950385, 950387, 950388, and 950407) (Fig. 1) that were not affected by the L-SSE.The analytical procedure is briefly described according to Yoshioka et al. (2015).
After steps associated with antenna exchange were taken, coseismic steps and common-mode errors were removed, and crustal displacements due to plate motion and annual and semiannual seasonal variations were estimated assuming the following model and removed from the GNSS time series: (1) where T = 1 year and the first and second terms on the right-hand side represent crustal deformation due to interplate locking.The remaining terms on the right-hand side in Eq. ( 1) represent annual and semiannual seasonal variations.Coefficients a to f are unknown parameters determined by a least-squares fit to the GNSS time series from January 1, 2017 (2017.0) to June 30, 2018June 30, (2018.5).5), during which no L-SSEs occurred (Fig. 2).Although Blewitt and Lavallée (2002) suggested that a time span of 2.5 years be adopted as a standard minimum data span for velocity solutions intended for tectonic interpretation, taking a period of this length is difficult in the Hyuganada region because the recurrence interval of L-SSEs in the region is 2-3 years (e.g., Ozawa 2017;Takagi et al. 2019).In fact, southward displacements that occurred in 2016 were identified at many stations (Fig. 2a, d); these displacements were transient signals associated with an L-SSE analyzed by Ozawa (2017).Therefore, we were not able to extend the period for estimating the trend to the preceding time period.Map plots of the linear trends for the horizontal and vertical components are shown in Additional file 1: Figs.S2(a) and (b), respectively.Spatial consistency is evident for both components, indicating that linear trends are removed properly.A curve fitting using the superposition of 20 cubic B-spline functions was performed for the corrected time series data for each component at every station, whose optimal smoothness was determined using inversion analyses based on the Akaike's Bayesian Information Criterion (ABIC) minimization principle (Akaike 1980) objectively and uniquely.We used these smooth functions to create piecewise displacement data for the spatiotemporal slip inversion of the L-SSE.By using the obtained best-fit curve (see Fig. 2) and the scattered processed time series data, we estimate a standard error of 1σ for the total displacement inferred in each time window.
Figures 4 and 5 show such horizontal and vertical displacements (arrows and bars) and related uncertainties (ellipses and bars), respectively.The standard errors were also used as weights for the GNSS displacements at each time window for the L-SSE slip inversion.Following Seshimo et al. (2023), we determined the minimum time window length so that the signal-to-noise (S/N) ratio of both the N-S and E-W components was 1.0 or greater.However, the first and last time windows, which are thought to contain almost no crustal deformation associated with the L-SSE, naturally have poor S/N ratios and do not meet these conditions.Consequently, based on the horizontal data, we selected time windows with a duration of 0.8 years (Fig. 4).Although the S/N ratio is not satisfied for vertical data when using 0.8-year time windows (Fig. 5), our main objective is to explain horizontal crustal displacements.
The period of analysis for the time series data used for the inversion analysis was set to 4.8 years, from September 13, 2017 (2017.7) to June 30, 2022 (2022.5),which includes immediately before and after the periods with little transient deformation.Using the continuous best-fit curve, each of the three-component time series at each station was divided into six equal datasets.As a result, we used all 468 datasets (= 3 (components) × 26 (stations) × 6 (0.8-year time windows)), which are shown in Figs. 4 and 5, for our inversion analysis.In the inversion analysis, the average standard errors were used as relative weights for the EW, NS, and UD components, whose values were determined to be 0.19 cm, 0.16 cm, and 0.67 cm, respectively, based on the time series analysis.

Inversion method
In this study, we estimated the spatiotemporal slip distribution using the inversion method of Yoshioka et al. (2015) that was subject to three prior constraints: (1) a smooth spatial distribution of slip, (2) a slip oriented mainly in the plate convergence direction (but without imposing a nonnegative constraint), and (3) a smooth temporal variation in slip.The slip distribution is represented by a superposition of bicubic B-spline functions, and the temporal slip variations are represented by a superposition of first-order B-spline functions.According to Yoshioka et al. (2015), the relationship between the model parameters representing the spatiotemporal slip distribution, the observed data and prior constraints is expressed by the following equation: where d is the observed displacement data vector; H is a matrix representing the relationship between the unit slip on the model fault plane and the displacement at each observation point; a is the model parameter vector; and A, B, and G are matrices representing the above-described constraints (1), (2), and (3), respectively.By minimizing the ABIC (Akaike 1980), we obtain the optimal values of the hyperparameters α = 3.7 × 10 −2 , β = 2.0 × 10 −1 , and γ = 1.6 × 10 −2 that represent the relative weights between the prior and observational constraints.We evaluated the reliability of the obtained slip distribution by using the estimation error and resolution of the estimated slip history, which were derived from the covariance matrix and resolution (2) (See figure on next page.)Fig. 2 GNSS time series data, including displacements associated with long-term slow slip events (SSEs), at all 26 stations (Fig. 1) used for the inversion analysis.The black dots indicate processed time series data after removing coseismic steps, steps associated with antenna exchange, linear trends, annual and semiannual variations, and common-mode errors.Each yellow-green curve is the best-fit curve for the processed data.The time span between the two pink vertical dashed lines is 2018.5-2021.7,during which the L-SSE was considered to have taken place.The location of each station is shown in Additional file 1: Table S1 matrix, respectively; the exact forms of these variables are described in Yoshioka et al. (2015).

Model setup
Inversion analysis was performed for crustal deformation for the aforementioned 4.8-year period by analyzing each of the subdivided 0.8-year time windows.Since an L-SSE is an aseismic slip phenomenon occurring on a plate boundary, the model source region was parametrized as a 3D nonplanar surface following the geometry of the upper plate boundary of the PHS plate by Hirose et al. (2008).Green's functions relating fault slip to surface displacements were calculated assuming a semi-infinite isotropic elastic body with linearity (Maruyama 1964), which enabled us to perform a linear inversion analysis based on Eq. ( 2).However, in reality, the deeper plate interface is a contact zone between the PHS slab and the mantle wedge, and steady slip is considered to occur there under the ductile regime caused by high temperature (Peacock and Wang 1999).Since the depth of the continental Moho is approximately 30 km in the Kyushu region (e.g., Murakoshi 2003;Abe et al. 2013) and the aseismic afterslip zone of the two 1996 Hyuganada earthquakes estimated by Yagi et al. (2001) reached 40 km depth, the maximum downdip depth of the model source region was set to approximately 65 km (Fig. 1).Therefore, the model source region was set to be 200 km in length in the strike direction and 100 km in width in the dip direction on the 3D plate interface.As GNSS data exist only on land, the resolution of slip decreases for shallower portions of the plate interface; thus, detailed, unbiased and statistically meaningful shallow slip cannot be inferred, regardless of the inversion approach (e.g., Ortega-Culaciati et al. 2021).
The slip distribution is represented by 11 × 7 bicubic B-spline basis functions distributed throughout the fault surface, where such a number of basis functions were determined following Yoshioka et al. (2004).When the number of basis functions increases, the estimation errors also tend to increase.Therefore, in this study, we determined the spatial number of basis functions appropriately, considering this tendency.Temporal changes in slip are represented by first-order B-spline basis functions, which are shaped as isosceles triangles covering two time windows.Therefore, we use seven B-spline functions for the six time windows.Slip is represented by two components, one in the strike direction and the other in the dip direction, resulting in 1078 (= 11 × 7 × 2 × 7) model parameters representing the spatiotemporal slip distributions.

Crustal deformations associated with the 2018-2021
Hyuga-nada L-SSE Figure 2 shows the processed time series data of the three components at the 26 GNSS stations used for the inversion analysis (Fig. 1, Additional file 1: Table S1).Most of the horizontal time series indicate that southward and eastward displacements began to be observed at approximately 2018.5 and ended at approximately 2021.7.During this period, the trend changed relatively monotonously.If there was a period of displacement whose slope of the time series was almost zero during the period of 3.2 years, there was little or no slippage of the L-SSE.However, we did not identify such signals during this period.Therefore, transient aseismic slow slip continued for the 3.2-year period, resulting in the longest duration ever reported.
Figure 3a shows the horizontal displacement field during the period from 2018.5 to 2021.7.This figure shows the dominance of displacements toward the east and southeast directions.Vertical displacements during the same period (Fig. 3b) were characterized by uplift along the eastern coastal region of Miyazaki Prefecture, reaching 5.2 cm.The systematic east-southeast-to southeastoriented horizontal displacements, which were in the almost opposite direction to that of plate convergence, and uplift along the eastern coast provided strong evidence of L-SSE occurrence.
We divided the data in Fig. 3a, b into equal time windows of 0.8 years in duration.In Fig. 4, we show the horizontal displacement field for each time window, including the windows immediately before and after the L-SSE period.Small displacements at all stations can be observed in Fig. 4a.This is because this time window was included in the detrended period from 2017.0 to 2018.5.In the following periods from 2018.5 to 2021.7, eastsoutheast to southeastward displacements were clearly identified, especially at stations along the eastern coast (Fig. 4b-e).In the subsequent time window from 2021.7 to 2022.5 (Fig. 4f ), small displacements reappeared at nearly all stations.Therefore, L-SSE was considered to have taken place during the period from 2018.5 to 2021.7.
The vertical displacement fields at the 0.8-year time window are shown in Fig. 5.The standard errors of the vertical displacements are larger than those of the horizontal displacements (Fig. 4).Since the period of the first time window (2017.7-2018.5) is fully included in the detrended period (2017.0-2018.5), the vertical displacements in Fig. 5a tend to be much smaller than the vertical displacements during the following time windows (Fig. 5b-e).However, because of the larger fluctuations in the vertical displacement data than in the horizontal displacement data, deviations in the vertical displacements from the 1.5-year linear trend appear for the 0.8-year time window (Fig. 5a).However, it should be noted that the standard errors of 1σ are larger than the vertical displacements at nearly all the stations.In the subsequent time windows from 2018.5 to 2021.7, uplifts were observed in the eastern coastal areas (Fig. 5b-e).In the following time window, this trend changed and turned into slight subsidence, which was also mostly below the standard error of 1σ (Fig. 5f ), indicating the end of the L-SSE by 2021.7.

Spatiotemporal slip distribution associated with the 2018-2021 Hyuga-nada L-SSE
The estimated spatiotemporal slip distributions are shown in Fig. 6 (see also Additional file 1: Fig. S3 for the distributions of estimation errors of 1σ).By comparing the observed displacements and slip inferences in the six time windows, we observed the following temporal change.In Fig. 6a, as expected from Figs. 4a and 5a, the slip amounts were less than the estimation errors of 1σ in most of the model source region.In the subsequent four time windows, statistically meaningful slips were inferred almost everywhere in the model source region (Fig. 6be).A clear, large slip of more than 2 cm was obtained in the central part of Miyazaki Prefecture, ranging in depth from approximately 20 to 40 km during the period of 2018.5-2019.3 (Fig. 6b).The overall slip directions in the third time window (Fig. 6c) were similar to those in the previous time window (Fig. 6b), but the slip amount increased, and the slip slightly expanded to the southern part of Miyazaki Prefecture.In the subsequent period of 2020.1-2020.9 (Fig. 6d), the maximum inferred slip reached 3.9 cm, and the main slip area was located at the same location and slip directions were almost the same as those in the previous time window (Fig. 6c).In the subsequent period, the overall slip amount decreased (Fig. 6e).
As was the case in Fig. 6a, in Fig. 6f, the slip amounts were less than the estimation errors of 1σ in most of the model source region.An annual average maximum slip rate of approximately 4.9 cm/yr (3.9 cm/0.8 yr) was estimated during the period from 2020.1 to 2020.9 (Fig. 6d).
It is difficult to directly connect the S/N ratios of the observed horizontal and vertical displacements with statistically meaningful inverted spatiotemporal slip distributions.In this study, we determined the length of the minimum time window so that the average S/N ratios of both the N-S and E-W components could be greater than 1.0 (Seshimo et al. 2023).However, the procedure resulted in S/N ratios smaller than 1.0 for the vertical data and hence in slip distributions constrained mostly by the horizontal observations.Nonetheless, this procedure for determining the duration of the time window proves useful for obtaining statistically meaningful spatiotemporal slip distributions where most of the slip amounts are larger than the estimation errors of 1σ.By summing the slip amounts inferred at the four time windows (Fig. 6b-e), we obtained the total slip distribution associated with the L-SSE (Fig. 7) (see also Additional file 1: Fig. S4 for the distribution of the estimation error of 1σ).The maximum total slip was estimated to be approximately 12.9 cm, and the equivalent release moment of the L-SSE was estimated to be 4.9 × 10 19 Nm, assuming a rigidity of 30 GPa, corresponding to Mw7.1.These release moment and moment magnitude were larger than those of previous L-SSEs in the region.
Figure 4 shows a comparison between the observed surface displacements at the GNSS stations and those predicted by the slip distributions (Fig. 6).Overall, the observed horizontal displacements at the stations were well explained by the calculations for the six time windows.For the vertical displacements, except for those in the eastern coastal region, the fit of the calculations to the observations was not good (Fig. 5).The main reason for the poor fit is that the weight of the vertical data was lower than that of the horizontal data at the time of the inversion analysis and because of the low (< 1.0) S/N ratio of the vertical data.The other reason is the spatial inconsistency of the vertical displacement field.On the other hand, the observed uplift patterns at stations along the eastern coastal region in southeastern Kyushu can be well explained by the calculations (Fig. 5b-e), which are useful for determining the depth extent of the slip distributions.

Comparison with past L-SSEs in the Hyuga-nada region
L-SSEs have occurred in the Hyuga-nada region in the past.Yarai and Ozawa (2013) estimated the spatiotemporal evolution of aseismic slow slip from January 1, 1997 to July 20, 2010 constrained by GNSS data (Fig. 8a, Table 1).They reported that three L-SSEs have occurred since 2004: the first occurred between January 1, 2005 and January 1, 2006 (2005 L-SSE), the second occurred between March 1996 to September 2017 (Fig. 8c) using a newly developed method for GNSS time series data.They also showed that the total slip amount was at a maximum in the southern Hyuga-nada region, reaching 79 cm among the four investigated segments in the subduction zone from southeastern Miyazaki Prefecture to western Shikoku.They also suggested that the L-SSE recurrence intervals were 2-3 years, which is consistent with the findings of Yarai and Ozawa (2013) and Ozawa (2017).
The annual average maximum slip rate for each L-SSE reported in Yarai and Ozawa (2013) was estimated to be 5-6 cm/yr for the 2005 L-SSE, 7-8 cm/yr for the 2007 L-SSE, and 4-5 cm/yr for the 2009 L-SSE (Table 1).Ozawa (2017) obtained approximately 8 cm/yr and 10 cm/yr annual average maximum slip rates for the first and second L-SSEs, respectively, in the southern part of Miyazaki Prefecture.Thus, a combined maximum slip of approximately 16 cm developed at that location between January 1, 2013 and April 10, 2016.Among the 12 L-SSEs detected in the Hyuga-nada region by Takagi et al. (2019), the annual maximum slip rate was estimated to be 15.5 cm/yr.In this study, the annual average maximum slip rate and maximum total slip were estimated to be approximately 4.9 cm/yr and 12.9 cm, respectively.Thus, these values are close to the values of the 2009 L-SSEs reported by Yarai and Ozawa (2013) and the combined maximum slip reported by Ozawa (2017), respectively.Compared with other L-SSEs, the annual maximum slip rate obtained in this study is relatively low, indicating gradual strain release.
On the other hand, the duration of all L-SSEs that occurred in the Hyuga-nada region was less than 1.7 years, and their moment magnitudes ranged from Mw6.0 to Mw7.0 (Table 1) (Yarai and Ozawa 2013;Ozawa 2017;Takagi et al. 2019).In this study, the estimated duration and moment magnitude of the L-SSE were approximately 3.2 years and Mw7.1, respectively.The duration was the longest, and the moment magnitude was the largest among the L-SSEs that occurred in the Hyuga-nada region.Ozawa (2017) reported that the L-SSE in the northern part of Miyazaki Prefecture occurred at the same time as or after that in the southern part of Miyazaki Prefecture; however, in this study, the slip patch in the central part of Miyazaki Prefecture (Fig. 6b) started slightly earlier than the southern slip patch (Fig. 6c-e).

Depth profile of slow earthquakes in the Hyuga-nada
region and their possible friction properties Yagi et al. (2001) estimated the afterslip areas of the October 19, 1996 and December 3, 1996 Hyuga-nada earthquakes (Mw6.8 and Mw6.7, respectively) using GNSS data from Kyushu.The characteristic decay times (= final slip/initial slope) were 15 and 100 days, respectively.Figure 7 shows the spatial distribution of cumulative slip over the 3.2 years (2018.5-2021.7)identified in this study,  together with the slip areas of the two 1996 Hyuga-nada earthquakes and their afterslip areas on the plate boundary (Yagi et al. 2001).The slip area that had a slip greater than approximately 10 cm on the Hyuga-nada L-SSE estimated in this study overlaps with the afterslip area of the Hyuga-nada earthquake that occurred on December 3, 1996.Since the slip area of the L-SSE inferred in this study almost coincides with the occurrence locations of the past L-SSEs in the Hyuga-nada region, it can be concluded that the past L-SSEs in the Hyuga-nada region also occurred in the afterslip areas of the 1996 Hyuga-nada earthquakes.This result is consistent with the findings of Yarai and Ozawa (2013), who inferred that the afterslip of the 1996 earthquakes and previous L-SSEs can occur in the same place.Thus, a similar frictional constitutive law is suggested for the occurrence of both slip types (L-SSE and afterslip) under the same temperature and pressure conditions.However, the duration of the L-SSE estimated in this study was approximately 3.2 years, which is much longer than the characteristic time of the afterslip of the December 3, 1996 Hyuga-nada earthquake.Asano et al. (2015) detected shallow very-lowfrequency earthquakes (VLFEs) off the coast of the Miyazaki Prefecture from January 1, 2010 to March 31, 2010 (Fig. 7).Yamashita et al. (2015) and Yamashita et al. (2021) determined the locations of epicenters of shallow tectonic tremors from April 2013 to July 2013 and from March 2014 to February 2017, respectively, which were mostly overlapped with the epicenters of the VLFEs.Since these observation periods were different from  2013).The gray area represents the area with an average resolution of 0.08 or less, as shown in Fig. 6.The gray lines represent isodepth contour lines of the upper surface of the PHS plate.The pink line represents the location of the mantle wedge corner (MWC), which coincides with the isodepth contour of 30 km.The yellow-green line represents the 350 °C isotherm of the upper surface of the subducted PHS plate by Ji et al. (2016).The solid blue area and the blue contour lines with 2 cm intervals represent the coseismic slip area and afterslip area of the October 19, 1996 Hyuga-nada earthquake, respectively (Yagi et al. 2001).The solid red and red contour lines with 2 cm intervals represent the coseismic slip area and afterslip area, respectively, of the December 3, 1996 Hyuga-nada earthquake (Yagi et al. 2001).The yellow-green solid square indicates the epicenter of the Mw6.2 Hyuga-nada earthquake on May 10, 2019 together with its CMT solution.The small solid orange stars represent the epicenters of the shallow very long frequency earthquakes (VLFEs) reported by Asano et al. (2015) from January 1, 2010 to March 31, 2010.The small gray and red solid circles represent the epicenters of shallow tectonic tremors reported by Yamashita et al. (2015) from April 2013 to July 2013 and Yamashita et al. (2021) from March 2014 to February 2017, respectively our period of analysis, it would be difficult to discuss the causal relationships among the L-SSE, VLFEs, and shallow tectonic tremors.However, based on Fig. 7, the L-SSE appeared to have formed at a much greater depth than the shallow VLFEs and shallow tectonic tremors.
Based on numerical simulations of 2D thermal structures along an across-arc profile passing through western Shikoku (Fig. 1), Gao and Wang (2017) attempted to explain the depth-dependent separation of slow earthquakes in the Nankai subduction zone (Obara and Kato 2016), which was described in Sect. 1. Gao and Wang (2017) suggested that the shear stress at which an L-SSE occurs is controlled by either semifrictional or viscous behavior depending on the thermal regime.They also noted that an L-SSE requires that the frictional behavior be very mildly velocity-weakening, as represented by a slightly negative value of the friction parameter (a -b), which is the velocity dependence of steady-state friction, according to laboratory experiments on rock friction (e.g., Ruina 1983;Dieterich 1994;Scholz 1998).
Unlike in Shikoku, in offshore Kyushu, L-SSEs and short-term SSEs (S-SSEs) occur at depths ranging from approximately 30 to 40 km (Okada et al. 2022) and are not separated.Furthermore, in the Hyuga-nada region, the afterslip of the 1996 Hyuga-nada earthquakes also occurred within the same depth range (Yagi et al. 2001;Yarai and Ozawa 2013) (Fig. 7).In the shallow region, significant VLFEs and shallow tremors occurred (Asano et al. 2015;Yamashita et al. 2015Yamashita et al. , 2021)).Therefore, in Hyuga-nada, the depth dependence of friction properties inferred from the conceptual model of Gao and Wang (2017) does not hold.Scholz (1998) focused on the depth dependence of the critical stiffness and argued that the above-mentioned transient aseismic slip occurs when the critical stiffness is conditionally stable and thus has a small positive value.Such regions exist within a depth range limited to the shallower and deeper sides of the so-called seismogenic zone of large interplate earthquakes.Applying this finding to the Hyuga-nada region, it is possible that the significant VLFE and tectonic tremors occurring in the shallow area are conditionally stable on the shallow side.On the other hand, regions where L-SSE, S-SSE, and afterslip occur may correspond to conditionally stable conditions on the deep side.
Regarding condition (1), namely, (a -b), since the dip angle in the deeper zone in the Hyuga-nada region is greater than that in the Nankai subduction zone, this effect contributes to producing a greater thermal gradient in the former than in the latter in the across-arc direction.Suenaga et al. (2018) performed 2D box-type time-dependent thermal modeling associated with the subduction of the PHS plate along an across-arc profile passing through the Hyuga-nada region to investigate the generation mechanisms of interplate seismic events.They found that the temperature range of the upper surface of the PHS plate, where the afterslip of the 1996 Hyuga-nada earthquake occurred, reached approximately 300-350 °C.According to Hyndman et al. (1997), this temperature range is located at the downdip end of the coseismic slip area, corresponding to the brittle-ductile and unstable to stable sliding transition zone.The depth of the continental Moho beneath Kyushu is estimated to be approximately 30 km based on receiver function analysis (Abe et al. 2013).Therefore, the mantle wedge corner (MWC) is located approximately parallel to and offshore from the east coast of Kyushu (Figs. 1  and 7).According to a three-dimensional thermal model associated with subduction of the PHS plate in southwest Japan by Ji et al. (2016), the downdip limit of the seismogenic zone in the Hyuga-nada region is controlled not by temperature but by the MWC because the temperature at the MWC is less than 350 °C (Fig. 7) and because the width of the seismogenic zone tends to narrow toward southern offshore Kyushu.This means that the temperature difference at the plate boundary between the shallow and deep sides of the MWC increases toward the south.Therefore, these effects play an important role in narrowing (a-b) and thus favoring conditionally stable zones.
Regarding condition (2), namely, pore pressure, based on the results of seismic tomography and receiver function analysis, a region with a high Poisson's ratio exists at depths greater than 30 km in the MWC, indicating the presence of serpentinized mantle and suggesting its relationship with the occurrence of transient aseismic slip (Tahara et al. 2008;Abe et al. 2013).
Regarding condition (3), namely, normal stress, Hyuganada is characterized by a strong negative Bouguer gravity anomaly reaching -140 mgal, and numerical modeling has shown that approximately 60% of this anomaly can be explained by the buoyancy of the subducted KPR (Tahara et al. 2008).This means that the buoyancy of the subducted KPR acts to strengthen the normal stress at the plate boundary in its exit region.
Therefore, (2) pore pressure and (3) normal stress, both of which are parameters of critical stiffness, increase at the plate boundary on the deep and shallow sides of the MWC, respectively.Therefore, the critical stiffness increases rapidly on the shallow side of the MWC and decreases on the deep side, contributing to a sharp depth contrast of the critical stiffness at the MWC.
Furthermore, the Kyushu-Palau Ridge (KPR) is subducted in the Hyuga-nada region, and the outer edge of the subducted KPR has been obtained by analyzing active and passive seismic sources (Yamamoto et al. 2013) (Fig. 7).Interestingly, this northwestern margin is approximately at 30 km depth and approximately coincident with the MWC, and the coseismic slips of the two 1996 Hyuga-nada earthquakes occurred within the KPR, whereas their afterslip occurred at the deeper plate boundary outside of the KPR.
Additionally, during the period of analysis, on May 10, 2019, a low-angle thrust earthquake with a magnitude of Mw6.2 occurred on the plate interface in the Hyuga-nada region (Fig. 7).The coseismic steps associated with the earthquake were corrected by calculating the average displacements for the 10 days before and after the earthquake and subtracting them from the original time series data of the three components.This event also occurred within the KPR and was located updip of the main slip area of the L-SSE, although the causal relationship between them is not clear.The earthquake occurred at the beginning of the third time window (2019.3-2020.1),during which clear large slips occurred in the deeper portion.Unlike the two 1996 Hyuga-nada earthquakes, an afterslip associated with the 2019 event with temporal decay did not seem to have taken place.The reasons might be that the moment magnitude of the earthquake was smaller than that of the 1996 event, that the earthquake occurred far from the GNSS stations on land and that the earthquake occurred within the unstable sliding zone.
Therefore, in the Hyuga-nada region, (1) rapid temperature changes at the plate boundary between the shallow and deep sides of the MWC, which may also be related to a high dip angle, play an important role in narrowing the conditionally stable zone.Additionally, (2) high pore pressure and serpentinization on the deep side of the MWC and (3) the increase in normal stress due to the buoyancy of the KPR at the shallow side of the MWC, which is coincident with the deepest northwestern rim of the subducted KPR, may contribute to the sharpness of the critical stiffness at the MWC.This sharp contrast in these factors at the shallow and deep sides of the MWC may limit the narrow depth range within which conditionally stable slip can be achieved.For these reasons, it is inferred that deep transient aseismic slips (L-SSE, S-SSE, and afterslip) occur within the same depth range.

Conclusions
In this study, we used GNSS time series data from January 1, 2017 to June 30, 2022 to detect a long-term slow slip event that occurred in the Hyuga-nada region, and estimated the spatiotemporal slip distribution.For this analysis, we used the inversion method proposed by Yoshioka et al. (2015).The significant results revealed in this study are summarized as follows: 1.An L-SSE was identified in the southern part of Miyazaki Prefecture for approximately 3.2 years around from 2018.5 to 2021.7, which was the longest duration for an L-SSE ever reported in the Hyuganada region.The maximum total slip was estimated to be approximately 12.9 cm.The release moment of the event was estimated to be 4.9 × 10 19 Nm, corresponding to Mw7.1, which was the largest L-SSE that occurred in this region.2. For the temporal evolution of the inferred slip, a slip of more than 2 cm was present in the central part of Miyazaki Prefecture during the period of 2018.5-2019.3.In the subsequent period from 2019.3 to 2020.1, the amount of slip increased and expanded to the southern part of Miyazaki Prefecture.Slip was the greatest in the subsequent period (2020.1-2020.9),reaching 3.9 cm.The slip area was almost the same in the following 0.8 years, and the overall slip amount decreased.The annual average maximum slip rate was approximately 4.9 cm/yr (3.9 cm/0.8 yr) during the period from 2020.1 to 2020.9. 3. The slip area of more than approximately 10 cm in the Hyuga-nada L-SSE region estimated in this study almost overlapped with the afterslip area of the December 3, 1996 Hyuga-nada earthquake, which had a depth range of approximately 30-40 km.S-SSEs also occurred within the same depth range.These coincidental same depth ranges of occurrence for interplate seismic events originate from large thermal gradients.The large thermal gradients play an important role in narrowing the depth range of (a-b) at the plate interface.In addition, the pore pressure and normal stress are closely related to the critical stiffness.
A sharp contrast in pore pressure and normal stress between the shallow and deep sides of the mantle wedge corner may also contribute to transient aseismic slip within the same depth range.

Fig. 1
Fig. 1 Tectonic map in and around the Hyuga-nada region.The thick barbed black line denotes the Nankai Trough.The black dashed lines in Kyushu are the boundaries of the prefectures.The thin black solid lines represent isodepth contours of the geometry of the upper surface of the PHS plate subducting beneath Kyushu: the 10 km isodepth by Baba et al. (2002), the 20-50 km isobath by Hirose et al. (2008), and the isodepth contours deeper than 60 km by Nakajima and Hasegawa (2007).The light green circle represents the approximate location of the Bungo Channel.The green bold line indicates the outer edge of the subducted Kyushu-Palau Ridge (KPR) by Yamamoto et al. (2013).The red and blue solid circles indicate the GNSS observation stations analyzed in this study whose detailed information is listed in Additional file 1: Table S1, and the blue solid squares represent the reference stations used in this study.The yellow-green box represents the horizontal projection of the source region of the assumed model on the 3D plate boundary.The red solid star represents the epicenter of the Mw7.0 Kumamoto earthquake on April 14, 2016.Processed time series data for all the stations (red and blue solid circles) with station codes are shown in Fig. 2. The gray boxed area in the inset shows our study area

Fig. 3
Fig. 3 Displacement fields associated with the Hyuga-nada L-SSE during the period from 2018.5 to 2021.7.The solid green circle indicates the observation station that exhibits a maximum displacement of 2.3 cm in the figure.The white solid circles represent stations strongly affected by the 2018-2019 Bungo Channel L-SSE and the 2016 Kumamoto earthquake (940090: Bungo Channel L-SSE, 960702: Kumamoto earthquake).(a) Horizontal displacement field.The red arrows indicate the stations used in the inversion, and the gray arrows indicate the stations not used in the inversion.The ellipse at the tip of each arrow represents the standard error of 1σ.(b) Vertical displacement field.The short horizontal line denotes the station location, and the red and blue vertical bars represent uplift and subsidence, respectively.The gray bars indicate the stations not used in the inversion.The gray thin vertical line at the tip of the observed bar represents a standard error of 1σ

Fig. 4
Fig.4Spatial distributions of horizontal displacements in southeastern Kyushu during the period from 2017.7 to 2022.5 with a 0.8-year time window.The red arrows represent observed processed displacements, and the blue arrows represent displacements calculated from the spatiotemporal slip distributions in Fig.6.The ellipse at the tip of each red arrow represents the standard error of 1σ. a 2017.7-2018.5.b  2018.5-2019.3.c 2019.3-2020.1.d 2020.1-2020.9. e 2020.9-2021.7.f 2021.7-2022.5

Fig. 5
Fig. 5 Spatial distributions of vertical displacements in the southeastern Kyushu during the period from 2017.7 to 2022.5 with a 0.8-year time window.The short black horizontal line represents the location of the station, and the bars above and below the horizontal line indicate uplift and subsidence, respectively.The red and blue bars represent the observed uplift and subsidence, respectively.The orange thin vertical line at the tip of the observed bar represents a standard error of 1σ.The yellow and light blue bars represent uplift and subsidence, respectively, calculated from the spatiotemporal slip distributions in Fig. 6.The periods of a-f are the same as those in Fig. 4

Fig. 6
Fig. 6 Slip distributions in the Hyuga-nada region with 0.8-year time windows from 2017.7 to 2022.5.The gray lines represent isodepth contour lines of the upper surface of the PHS plate.The black dashed lines in southeastern Kyushu are the boundaries of the prefectures.The arrows indicate the directions and amounts of slip on the model source region, and the arrow-tip circles represent an estimation error of 1σ.The solid black contours represent slip amounts with 1 cm intervals.The gray area denotes the area with a resolution of 0.08 or less.The periods of a-f are the same as those in Fig. 4

Fig. 7
Fig. 7 Cumulative slip distribution during the period from 2018.5 to 2021.7 at 2 cm contour intervals.The green bold line indicates the outer edge of the subducted Kyushu-Palau Ridge (KPR) byYamamoto et al. (2013).The gray area represents the area with an average resolution of 0.08 or less, as shown in Fig.6.The gray lines represent isodepth contour lines of the upper surface of the PHS plate.The pink line represents the location of the mantle wedge corner (MWC), which coincides with the isodepth contour of 30 km.The yellow-green line represents the 350 °C isotherm of the upper surface of the subducted PHS plate byJi et al. (2016).The solid blue area and the blue contour lines with 2 cm intervals represent the coseismic slip area and afterslip area of the October 19, 1996 Hyuga-nada earthquake, respectively(Yagi et al. 2001).The solid red and red contour lines with 2 cm intervals represent the coseismic slip area and afterslip area, respectively, of the December 3, 1996 Hyuga-nada earthquake(Yagi et al. 2001).The yellow-green solid square indicates the epicenter of the Mw6.2 Hyuga-nada earthquake on May 10, 2019 together with its CMT solution.The small solid orange stars represent the epicenters of the shallow very long frequency earthquakes (VLFEs) reported byAsano et al. (2015) from January 1, 2010 to March 31, 2010.The small gray and red solid circles represent the epicenters of shallow tectonic tremors reported byYamashita et al. (2015) from April 2013 to July 2013 and Yamashita et al. (2021) from March 2014 to February 2017, respectively

Fig. 8
Fig.8The slip distributions of the Hyuga-nada L-SSEs in previous studies, which are also summarized in Table1.The pink line represents the location of the mantle wedge corner (MWC), which coincides with the isodepth contour of 30 km.(a) Slip distribution byYarai and Ozawa (2013).The dark red line represents the slip distribution of an L-SSE from January 1, 2005 to January 1, 2006, the red line represents the slip distribution of an L-SSE from January 1, 2007 to January 1, 2008, and the bright red line represents the slip distribution of an L-SSE from January 1, 2009 to July 20, 2010.The contour intervals of the slip amounts are 2 cm.(b) Slip distribution of L-SSEs by Ozawa (2017).The dark blue line represents the slip distribution of an L-SSE from October 1, 2013 to May 1, 2014, the blue line represents the slip distribution of an L-SSE from July 1, 2015 to February 1, 2016, and the light blue line represents the slip distribution of an L-SSE from February 1, 2016 to April 10, 2016.The contour intervals of the slip amounts are 4 cm.(c) Locations of L-SSEs represented by rectangular faults, which were detected by Takagi et al. (2019).The dashed line of each rectangular fault represents its shallower side.The black, dark blue, dark green, dark red, blue, green, red, orange, purple, yellow, light blue, and yellow-green rectangles indicate the estimated locations of rectangular fault planes for the L-SSEs whose onset dates wereSeptember 30, 1996, December 23, 1997, November 2, 1998, February 12, 2000, January 28, 2002, March 28, 2002, February 21, 2005, October 6,  2006, May 22, 2008, August 12, 2009, April 10, 2011, and July 15, 2013, respectively

Table 1
L-SSEs and large earthquakes detected in the Hyuga-nada region in previous studies