Kawabata et al.
Progress in Earth and Planetary Science(2024) 11:24
https://doi.org/10.1186/s40645-024-00624-z

Progress in Earth andPlanetary Science

RESEARCH ARTICLE

Open Access

Detection of a recent large Hyuga‑nadalong‑term slow slip event and estimation of its
spatiotemporal slip distributions
Hiroki Kawabata1* , Yukinari Seshimo1, Shoichi Yoshioka1,2 and Francisco Ortega‑Culaciati3,4

Abstract
Using horizontal and vertical GNSS time series data from the GSI in sotheastern Kyushu from January 1, 2017 to June30, 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 partof 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 × ­1019 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 morethan 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 rangeof 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 deepsides of the mantle wedge corner may also contribute to transient aseismic slip within the same depth range.
Keywords L-SSE, Hyuga-nada, Inversion analysis, Spatiotemporal slip distribution, Critical stiffness, Conditionallystable*Correspondence:
Hiroki Kawabata222s410s@gsuite.kobe-u.ac.jp
Full list of author information is available at the end of the article© The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 2 of 17

1 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-classmegathrust earthquakes have occurred at intervals of
approximately 90–150 years (e.g., Ando 1974; Kumagai1996), and many researchers around the world have conducted research on seismic events that occur at this plateboundary. Since the discovery of tectonic tremors in the
Nankai Trough subduction zone by Obara (2002) aroundthe turn of this century, studies on slow earthquakes havebeen actively conducted. Megathrust earthquakes andslow earthquakes are slip phenomena that occur on plateinterfaces. According to Obara and Kato (2016), shallow
very-low-frequency (VLF) and tremor events, megathrust earthquakes, long-term slow slip events (hereafterreferred to as L-SSEs), and deep episodic tremors andslips (ETSs) occur from shallower to deeper nonoverlapping depth segments across the Nankai Trough plateinterface.
The Hyuga-nada region, the target area of this study, islocated 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, theold (> 50 Ma) PHS plate is subducting at a high angle, andno M8-class megathrust earthquakes have been identified for at least the past 400 years (Earthquake ResearchCommittee 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 andtectonic tremors occur in and around shallow plate interface (e.g., Yamashita et al. 2015, 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 fromthose 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 fromJanuary 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 andYoshioka 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 plateboundary in the Hyuga-nada region.

2 DataIn 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 spansfrom January 1, 2017 to June 30, 2022. We selected
the GNSS stations via the following method: since wefocused only on detecting an L-SSE in the Hyuga-nadaregion 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 stationsin the northern half of the Kyushu region (see Additional file 1: Fig. S1(a)). We selected all the stations inthe southern half of the Kyushu region. However, inthe 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

Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 3 of 17

Fig. 1 Tectonic map in and around the Hyuga-nada region. The thick barbed black line denotes the Nankai Trough. The black dashedlines in Kyushu are the boundaries of the prefectures. The thin black solid lines represent isodepth contours of the geometry of the uppersurface 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 locationof 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 projectionof 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

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 byother unknown mechanisms. Therefore, we removed
the data from such stations from the dataset (see graydisplacements in Fig. 3a, b) and used the data from theKawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 4 of 17

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

L-SSE occurrence
period andnumber of L-SSE
occurrences

Duration (year) Maximum
Mw
slip rate (cm/yr)

This study

From January 1,
2017 to June 30,
2022

From 2018.5
to 2021.7 (theunits are 0.1 years(= 36.5 days))

3.2

4.9

7.1

Yarai and Ozawa(2013)

From January 1,
1997 to July 20,
2010

From January 1,
2005 to January 1,
2006

1.0

5–6

6.7 to 6.8

From January 1,
2007 to January 1,
2008

1.0

7–8

From January 1,
2009 to January 1,
2010

1.0

4–5

From October 1,
2013 to May 1,
2014; two L-SSEs

0.58

8

6.4, 6.5

From July 1, 2015
to April 10, 2016;
two L-SSEs

0.78

10

6.0, 6.4

0.40–1.67

15.5a

6.4–7.0

Ozawa (2017)

From January 1,
2013 to April 10,
2016

Takagi et al. (2019)

From January 1,
From Septem‑
1996 to December ber 30, 1996
31, 2017
to November 8,
2014; 12 L-SSEs

Yagi et al. (2001)

Earthquake
occurrence date

Maximum amountof slip of anearthquake

October 19, 1996

6 cm

December 2, 1996

13 cm

a

The maximum slip rate was determined by dividing the uniform slip amount of each L-SSE on its rectangular fault by its duration, and the largest value among the 12
L-SSEs is listed here

remaining 26 GNSS stations shown in Fig. 1 (see also
Additional file 1: Table S1).

3 Methods
3.1 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 notaffected by the L-SSE. The analytical procedure is briefly
described according to Yoshioka et al. (2015).
After steps associated with antenna exchange weretaken, 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 fromthe GNSS time series:y(t) = a + bt + c sin+ e sin

4π t
T

2π t
T

+ d cos+ f cos

4πt
T

2π t
T

(1)

where T = 1 year and the first and second terms on theright-hand side represent crustal deformation due tointerplate locking. The remaining terms on the right-handside in Eq. (1) represent annual and semiannual seasonalvariations. 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, 2018 (2018.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 spanfor 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 in2016 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 notable to extend the period for estimating the trend to thepreceding 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

Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 5 of 17

consistency is evident for both components, indicatingthat 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 (Akaike1980) objectively and uniquely. We used these smoothfunctions to create piecewise displacement data for the
spatiotemporal slip inversion of the L-SSE. By using theobtained best-fit curve (see Fig. 2) and the scattered processed time series data, we estimate a standard errorof 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 toexplain 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 timeseries 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 tobe 0.19 cm, 0.16 cm, and 0.67 cm, respectively, basedon the time series analysis.
3.2 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 asuperposition 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:
  

H
d
 0   α̂A 
[a]
0=
(2)
β̂B
0
γ̂ G
where d is the observed displacement data vector; H isa 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 theabove-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−2that represent the relative weights between the prior andobservational constraints. We evaluated the reliabilityof the obtained slip distribution by using the estimation
error and resolution of the estimated slip history, whichwere derived from the covariance matrix and resolution(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 curvefor 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 havetaken place. The location of each station is shown in Additional file 1: Table S1. a–c Time series data at the GNSS stations indicated by the red solid
circles in Fig. 1. d–f Time series data at the GNSS stations indicated by the blue solid circles in Fig. 1. a, d North‒south component. b, e East‒west
component. c, f Vertical component

Kawabata et al. Progress in Earth and Planetary Science

80

(a)

960711

NS95047870(2024) 11:24

80

(b)

96071160

Displacement (cm)5095047930960715

80

2022

NS

950481

70

0

80

2016

2018

(e)

60

Displacement (cm)50

30

Fig. 2 (See legend on previous page.)(f)095048350
02P211021090409504843094200496071620

940099021093

2016

2018

2020

Time (year)UD950482021093109504912022

2022

95048160

940099102020960714960716950491Time (year)

2018

Time (year)

02108794200420021093

2020

80

2016

950484940099

2018

0

70

021090409420042016

EW

02P21196071602022

950482950484100210889504830210902096071510021087

02P21130960713021086950481950483400210852020209607147095048250

30Time (year)

021087609504800210882020

96071495047940960715

10

Time (year)(d)0210840210860210882018

50

960713200210862016

960712021085960713

0

02P122950480021085109504786095047940

UD94009402108495048020960711960712021084

30

(c)

70

02P1229607124080940094

02P12250EW95047870

94009460

Page 6 of 17950491
2022

0

2016

2018

2020

Time (year)2022

Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 7 of 17

matrix, respectively; the exact forms of these variables
are described in Yoshioka et al. (2015).

4 Results

3.3 Model setup

Figure 2 shows the processed time series data of the threecomponents at the 26 GNSS stations used for the inversion analysis (Fig. 1, Additional file 1: Table S1). Mostof the horizontal time series indicate that southwardand 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 slopeof 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 thisperiod. 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 andsoutheast directions. Vertical displacements during thesame 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 thealmost 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 thehorizontal displacement field for each time window,
including the windows immediately before and after the
L-SSE period. Small displacements at all stations can beobserved 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, east‒
southeast 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 atnearly all stations. Therefore, L-SSE was considered tohave 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 thevertical displacements are larger than those of the horizontal displacements (Fig. 4). Since the period of the firsttime window (2017.7–2018.5) is fully included in thedetrended 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

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. Sincean 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 asemi-infinite isotropic elastic body with linearity (Maruyama 1964), which enabled us to perform a linear inversion analysis based on Eq. (2). However, in reality, thedeeper 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 hightemperature (Peacock and Wang 1999). Since the depthof 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 forshallower 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 thefault surface, where such a number of basis functions
were determined following Yoshioka et al. (2004). Whenthe 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 inslip 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 otherin the dip direction, resulting in 1078 (= 11 × 7 × 2 × 7)
model parameters representing the spatiotemporal slip
distributions.

4.1 Crustal deformations associated with the 2018–2021
Hyuga‑nada L‑SSE

Kawabata et al. Progress in Earth and Planetary Science

130˚E
33˚N

(2024) 11:24

131˚E

132˚E
940090(a)960702

32˚N

1 cm
31˚N

2018.5-2021.7

50 km

33˚N

940090960702(b)

32˚N

1 cm (uplift)
31˚N
130˚E

1 cm (subsidence)

2018.5-2021.7

50 km
131˚E

132˚E

from the 1.5-year linear trend appear for the 0.8-yeartime window (Fig. 5a). However, it should be notedthat 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 wereobserved in the eastern coastal areas (Fig. 5b–e). In thefollowing time window, this trend changed and turned
into slight subsidence, which was also mostly below thestandard error of 1σ (Fig. 5f ), indicating the end of the
L-SSE by 2021.7.

Page 8 of 17

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 errorof 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 tipof the observed bar represents a standard error of 1σ

4.2 Spatiotemporal slip distribution associatedwith the 2018–2021 Hyuga‑nada L‑SSE

The estimated spatiotemporal slip distributions areshown 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 thesix time windows, we observed the following temporalchange. In Fig. 6a, as expected from Figs. 4a and 5a, theslip amounts were less than the estimation errors of 1σ inmost of the model source region. In the subsequent fourtime windows, statistically meaningful slips were inferredalmost everywhere in the model source region (Fig. 6b–e). 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 amountincreased, and the slip slightly expanded to the southernpart of Miyazaki Prefecture. In the subsequent periodof 2020.1–2020.9 (Fig. 6d), the maximum inferred slip
reached 3.9 cm, and the main slip area was located at thesame 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 themodel 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 theobserved horizontal and vertical displacements with statistically meaningful inverted spatiotemporal slip distributions. In this study, we determined the length of theminimum time window so that the average S/N ratiosof both the N–S and E–W components could be greater

Kawabata et al. Progress in Earth and Planetary Science

130˚

131˚

132˚

130˚

(2024) 11:24

131˚Page 9 of 17

132˚

130˚

131˚

(b)

(a)

132˚

(c)

32˚

32˚
0.5 cm
obs.
cal.

31˚

2017.7-2018.5

2018.5-2019.3

2019.3-2020.1

(d)

(e)

(f)

32˚

31˚
130˚

31˚

32˚

2020.1-2020.9
131˚

132˚

2020.9-2021.7
130˚

131˚

132˚

2021.7-2022.5
130˚

131˚

132˚

31˚

Fig. 4 Spatial distributions of horizontal displacements in southeastern Kyushu during the period from 2017.7 to 2022.5 with a 0.8-yeartime 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

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 procedurefor determining the duration of the time window proves
useful for obtaining statistically meaningful spatiotemporal slip distributions where most of the slip amounts arelarger than the estimation errors of 1σ.
By summing the slip amounts inferred at the fourtime windows (Fig. 6b–e), we obtained the total slip
distribution associated with the L-SSE (Fig. 7) (seealso 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 theequivalent release moment of the L-SSE was estimatedto be 4.9 × ­1019 Nm, assuming a rigidity of 30 GPa,corresponding to Mw7.1. These release moment andmoment 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 thosepredicted by the slip distributions (Fig. 6). Overall, theobserved 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 waslower than that of the horizontal data at the time of theinversion analysis and because of the low (< 1.0) S/N ratioof 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 theeastern coastal region in southeastern Kyushu can be well
explained by the calculations (Fig. 5b–e), which are usefulfor determining the depth extent of the slip distributions.

5 Discussion
5.1 Comparison with past L‑SSEs in the Hyuga‑nada region

L-SSEs have occurred in the Hyuga-nada region in thepast. 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

Kawabata et al. Progress in Earth and Planetary Science

130˚

131˚

132˚

130˚

(2024) 11:24

131˚

(a)

Page 10 of 17

132˚

130˚

131˚

132˚

(c)

(b)

32˚

32˚
0.3 cm (obs.)
0.3 cm (cal.)31˚

2017.7-2018.5

2018.5-2019.3

2019.3-2020.1

(d)

(e)

(f)

32˚

31˚
130˚

31˚

32˚

131˚

132˚

2021.7-2022.5

2020.9-2021.7

2020.1-2020.9
130˚

131˚

132˚

130˚

131˚

132˚

31˚

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 upliftand 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

January 1, 2007 and January 1, 2008 (2007 L-SSE), andthe third occurred between January 1, 2009 and July 20,
2010 (2009 L-SSE). Ozawa (2017) estimated the spatiotemporal aseismic slip distribution from January 1, 2013
to April 10, 2016 using GNSS data and reported two different rupture processes in relation to the L-SSEs that
occurred in the northern part of Miyazaki Prefecture andin the southern part of Miyazaki Prefecture (Fig. 8b andTable 1). The first rupture process was a simultaneous
L-SSE in the northern and southern parts of Miyazaki
Prefecture spanning from October 1, 2013 to May 1,
2015, along with the 2014 Bungo Channel L-SSE. The
second rupture process began on approximately July 1,2015 in the southern part of Miyazaki Prefecture, followed by the L-SSE in the northern part of Miyazaki
Prefecture and the Bungo Channel L-SSE in 2016. The
main slip area of the 2018–2021 L-SSE identified in thisstudy almost coincided with that identified by Yarai andOzawa (2013) and the slip patches of the southern partof Miyazaki Prefecture obtained by Ozawa (2017). Takagi
et al. (2019) detected 12 L-SSEs, assuming uniform slip
on a rectangular fault, in the Hyuga-nada region fromMarch 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 inthe 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 thefindings 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 and10 cm/yr annual average maximum slip rates for the firstand second L-SSEs, respectively, in the southern part ofMiyazaki Prefecture. Thus, a combined maximum slip ofapproximately 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

Kawabata et al. Progress in Earth and Planetary Science

131˚

132˚

133˚

(a)

(2024) 11:24

131˚Page 11 of 17

132˚

133˚

131˚

132˚

133˚

(c)

(b)

1

1

2

32˚

1

32˚
2

3

1

3.0 cm
31˚

50 40

30

2017.7-2018.5

20
3

(d)

2018.5-2019.3(e)2019.3-2020.1

31˚

(f)

1

1

32˚
1

32˚
1

2

3

2

1

1

1

31˚

2020.1-2020.9
131˚

132˚

133˚

2020.9-2021.7
131˚

132˚
0

2

4

slip[cm]

133˚

2021.7-2022.5
131˚

132˚

31˚

133˚

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 sameas those in Fig. 4

slip rate and maximum total slip were estimated to beapproximately 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 occurredin 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).
5.2 Depth profile of slow earthquakes in the Hyuga‑nadaregion 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,

Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 12 of 17

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 edgeof the subducted Kyushu-Palau Ridge (KPR) by Yamamoto 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 locationof the mantle wedge corner (MWC), which coincides with the isodepth contour of 30 km. The yellow-green line represents the 350 °C isothermof 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

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 greaterthan approximately 10 cm on the Hyuga-nada L-SSE
estimated in this study overlaps with the afterslip area ofthe Hyuga-nada earthquake that occurred on December
3, 1996. Since the slip area of the L-SSE inferred in thisstudy 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 inthe 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 estimatedin 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 fromMarch 2014 to February 2017, respectively, which were
mostly overlapped with the epicenters of the VLFEs.
Since these observation periods were different from

Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 13 of 17

Fig. 8 The slip distributions of the Hyuga-nada L-SSEs in previous studies, which are also summarized in Table 1. The pink line represents
the location of the mantle wedge corner (MWC), which coincides with the isodepth contour of 30 km. (a) Slip distribution by Yarai 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 were
September 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

our period of analysis, it would be difficult to discuss
the causal relationships among the L-SSE, VLFEs, andshallow 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 aslightly 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 andshort-term SSEs (S-SSEs) occur at depths ranging fromapproximately 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 (Asanoet al. 2015; Yamashita et al. 2015, 2021). Therefore, inHyuga-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 thecritical 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 zoneof 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.
As described before, L-SSEs and deep episodic tremorsand slips (ETSs) occur from shallower to deeper nonoverlapping depth segments across the Nankai Troughplate interface (e.g., Obara and Kato 2016). We consider

Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 14 of 17

why, unlike in the Nankai Trough, offshore Kyushu, transient aseismic slips, such as those in the Nankai Trough,occur within the same depth range in terms of the depth
dependency of parameters related to the critical stiffness:(1) (a – b), (2) pore pressure, and (3) normal stress.
Regarding condition (1), namely, (a – b), since the dip
angle in the deeper zone in the Hyuga-nada region isgreater than that in the Nankai subduction zone, thiseffect contributes to producing a greater thermal gradientin 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 theHyuga-nada region to investigate the generation mechanisms of interplate seismic events. They found that thetemperature 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 toHyndman 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 onreceiver function analysis (Abe et al. 2013). Therefore, themantle 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 bytemperature but by the MWC because the temperature at
the MWC is less than 350 °C (Fig. 7) and because the widthof the seismogenic zone tends to narrow toward southern
offshore Kyushu. This means that the temperature difference at the plate boundary between the shallow and deepsides of the MWC increases toward the south. Therefore,
these effects play an important role in narrowing (a–b) andthus favoring conditionally stable zones.
Regarding condition (2), namely, pore pressure, basedon 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 thepresence 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 modelinghas 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 theplate boundary in its exit region.
Therefore, (2) pore pressure and (3) normal stress, bothof which are parameters of critical stiffness, increaseat the plate boundary on the deep and shallow sides ofthe MWC, respectively. Therefore, the critical stiffness
increases rapidly on the shallow side of the MWC anddecreases 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 edgeof the subducted KPR has been obtained by analyzing active and passive seismic sources (Yamamoto et al.
2013) (Fig. 7). Interestingly, this northwestern margin isapproximately at 30 km depth and approximately coincident with the MWC, and the coseismic slips of thetwo 1996 Hyuga-nada earthquakes occurred within the
KPR, whereas their afterslip occurred at the deeper plateboundary 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-nadaregion (Fig. 7). The coseismic steps associated with theearthquake were corrected by calculating the average displacements for the 10 days before and after the earthquakeand subtracting them from the original time series data ofthe 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 themis not clear. The earthquake occurred at the beginningof the third time window (2019.3–2020.1), during whichclear large slips occurred in the deeper portion. Unlike thetwo 1996 Hyuga-nada earthquakes, an afterslip associatedwith the 2019 event with temporal decay did not seem tohave taken place. The reasons might be that the moment
magnitude of the earthquake was smaller than that of the1996 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 temperaturechanges at the plate boundary between the shallow and deepsides of the MWC, which may also be related to a high dip
angle, play an important role in narrowing the conditionallystable zone. Additionally, (2) high pore pressure and serpentinization on the deep side of the MWC and (3) the increasein 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Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 15 of 17

inferred that deep transient aseismic slips (L-SSE, S-SSE, andafterslip) occur within the same depth range.

Abbreviations
GNSS	Global navigation satellite system
GSI	Geospatial Information Authority of Japan
L-SSE	Long-term slow slip event
GEONET	GNSS earth observation network system
NIED	
National Research Institute for earth science and disasterprevention

6 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, andestimated the spatiotemporal slip distribution. For thisanalysis, we used the inversion method proposed by
Yoshioka et al. (2015). The significant results revealed inthis study are summarized as follows:
1. An L-SSE was identified in the southern part ofMiyazaki 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 estimatedto be approximately 12.9 cm. The release moment of
the event was estimated to be 4.9 × ­1019 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 slipof more than 2 cm was present in the central part ofMiyazaki Prefecture during the period of 2018.5–
2019.3. In the subsequent period from 2019.3 to2020.1, the amount of slip increased and expanded tothe southern part of Miyazaki Prefecture. Slip was thegreatest in the subsequent period (2020.1–2020.9),
reaching 3.9 cm. The slip area was almost the samein the following 0.8 years, and the overall slip amountdecreased. 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 thisstudy 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 forinterplate seismic events originate from large thermal
gradients. The large thermal gradients play an important role in narrowing the depth range of (a–b) at theplate interface. In addition, the pore pressure and normal stress are closely related to the critical stiffness.
A sharp contrast in pore pressure and normal stressbetween the shallow and deep sides of the mantle
wedge corner may also contribute to transient aseismic slip within the same depth range.

Supplementary Information
The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s40645-​024-​00624-z.
Additional file 1. Table S1. GNSS stations used in the inversion analysis
in this study. Figure S1. Time series data at GNSS stations affected by tran‑sient slip caused by the 2018-2019 Bungo Channel LSSE and the postseis‑mic deformation associated with the 2016 Kumamoto earthquake. Figure
S2. Map plots of the linear trend obtained by Eq. (1) at each GNSS station.
Figure S3. Distributions of the estimation error of 1σ for the inverted slip
distributions shown in Fig. 6. Figure S4. Distributions of estimation errorsof 1σ for the cumulative slip distribution shown in Fig. 7.
Acknowledgements
We thank T. Yabuki for sharing the source code of geodetic data inversion
used in Yabuki and Matsu’ura (1992). We also thank the Editor Y. Yagi, twoanonymous reviewers, R. Burgman and H. Hirose for their constructive com‑ments to improve the manuscript. We also thank Geospatial Information
Authority of Japan for allowing us to use daily GNSS time series data. Figures
were created by using the Generic Mapping Tools (GMT) (version: GMT 4.5.7,
URL link: https://​www.​gener​ic-​mappi​ng-​tools.​org/​downl​oad/) (Wessel andSmith 1998). F.O-C. acknowledges support from Data Observatory Foundation,
ANID Technology Center No. DO210001.
Author contributions
H.K. carried out the data processing of GNSS time series and carried out theinversion analysis, created the figures, and prepared some of the figures andsubmitted the paper. Y.S. performed preliminary data processing of GNSS
time series and the inversion analysis and wrote a draft. S.Y. organized andinstructed this study, pointing out possible problems, advised the solutions tothose in this study and updated the manuscript. F.O.-C. pointed out possible
problems with the study and advised about solutions, contributed to the writ‑ing and revision of the manuscript.
Funding
This work was partly supported by a MEXT KAKENHI grant [grant number
21H05203] and a Kobe University Strategic International Collaborative
Research Grant (Type B Fostering Joint Research).
Availability of data and materials
The original GNSS time series data used in this study are available at https://​www.​gsi.​go.​jp/​ENGLI​SH/​geonet_​engli​sh.​html

Declarations
Competing interests
The authors declare that they have no competing interest.
Author details
1
Department of Planetology, Graduate School of Science, Kobe University,
Rokkodai‑Cho 1‑1, Nada Ward, Kobe 657‑8501, Japan. 2 Research Centerfor Urban Safety and Security, Kobe University, Rokkodai‑Cho 1‑1, Nada
Ward, Kobe 657‑8501, Japan. 3 Department of Geophysics, Faculty of Physicaland Mathematical Sciences, University of Chile, Santiago, Chile. 4 Data Observa‑
tory Foundation, ANID Technology Center No. DO210001, Santiago, Chile.
Received: 28 September 2023 Accepted: 31 March 2024

Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Page 16 of 17

References
Abe Y, Ohkura T, Hirahara K, Shibutani T (2013) Along-arc variation in water
distribution in the uppermost mantle beneath Kyushu, Japan, as derived
from receiver function analyses. J Geophys Res Solid Earth 118:3540–
3556. https://​doi.​org/​10.​1002/​jgrb.​50257
Akaike H (1980) Likelihood and the Bayes procedure. In: Bernardo JM, DeGroot
MH, Lindley FV, Smith AFM (eds) Bayesian Statistics. University Press,
Valencia, pp 143–166
Ando M (1974) Source mechanisms and tectonic significance of historicalearthquakes along the Nankai Trough. Japan Tectonophys 27(2):119–140.
https://​doi.​org/​10.​1016/​0040-​1951(75)​90102-X
Asano Y, Obara K, Matsuzawa T et al (2015) Possible shallow slow slip eventsin Hyuga-nada, Nankai subduction zone, inferred from migration of very
low frequency earthquakes. Geophys Res Lett 42:331–338. https://​doi.​
org/​10.​1002/​2014G​L0621​65
Baba T, Tanioka Y, Cummins PR, Uhira K (2002) The slip distribution of the 1946
Nankai earthquake estimated from tsunami inversion using a new platemodel. Phys Earth Planet Inter 132:59–73. https://​doi.​org/​10.​1016/​S0031-​
9201(02)​00044-4
Blewitt G, Lavallée D (2002) Effect of annual signals on geodetic velocity. J
Geophys Res Solid Earth 107:ETG 9-1-ETG 9-11. https://​doi.​org/​10.​1029/​
2001J​B0005​70
DeMets C, Gordon RG, Argus DF (2010) Geologically current plate motions.
Geophys J Int 181:1–80
Dieterich J (1994) A constitutive law for rate of earthquake production andits application to earthquake clustering. J Geophys Res Solid Earth
99:2601–2618. https://​doi.​org/​10.​1029/​93JB0​2581
Earthquake Research Committee (2022) Long-term evaluation of seismic
activity around the Hyuga-nada and Nansei Islands Trench (2nd edition)(in Japanese). https://​www.​jishin.​go.​jp/​main/​chousa/​kaikou_​pdf/​hyuga​nada_2.​pdf
Gao X, Wang K (2017) Rheological separation of the megathrust seismogenic
zone and episodic tremor and slip. Nature 543:416–419. https://​doi.​org/​
10.​1038/​natur​e21389
Hirose F, Nakajima J, Hasegawa A (2008) Three-dimensional seismic velocity
structure and configuration of the Philippine Sea slab in SouthwesternJapan estimated by double-difference tomography. J Geophys Res Earth.
https://​doi.​org/​10.​1029/​2007J​B0052​74
Hyndman RD, Yamano M, Oleskevich DA (1997) The seismogenic zone of
subduction thrust faults. Island Arc 6:244–260. https://​doi.​org/​10.​1111/j.​
1440-​1738.​1997.​tb001​75.x
Ji Y, Yoshioka S, Matsumoto T (2016) Three-dimensional numerical modelingof temperature and mantle flow fields associated with subduction of
the Philippine Sea plate, southwest Japan. J Geophys Res Solid Earth
121:4458–4482. https://​doi.​org/​10.​1002/​2016J​B0129​12
Kumagai H (1996) Time sequence and the recurrence models for large earth‑quakes along the Nankai Trough revisited. Geophys Res Lett 23(10):1139–
1142. https://​doi.​org/​10.​1029/​96GL0​1037
Maruyama T (1964) Statical elastic dislocations in an infinite and semi-infinite
medium. Bull Earthq Res Inst Tokyo Univ 42:289–368
Murakoshi T (2003) Seismic imaging of the crust and uppermost mantle
structure beneath Kyushu by receiver function analysis. Dissertation,
Kyushu University.
Nakajima J, Hasegawa A (2007) Subduction of the Philippine Sea plate
beneath southwestern Japan: Slab geometry and its relationship to arc
magmatism. J Geophys Res 112:B08306. https://​doi.​org/​10.​1029/​2006J​
B0047​70
Obara K (2002) Nonvolcanic deep tremor associated with subduction inSouthwest Japan. Science 296:1679–1681. https://​doi.​org/​10.​1126/​scien​
ce.​10703​78
Obara K, Kato A (2016) Connecting slow earthquakes to huge earthquakes.
Science 353:6296. https://​doi.​org/​10.​1126/​scien​ce.​aaf15​12
Okada Y, Nishimura T, Tabei T, Matsushima T, Hirose H (2022) Development ofa detection method for short-term slow slip events using GNSS data andits application to the Nankai subduction zone. Earth, Planets Space 74:18.
https://​doi.​org/​10.​1186/​s40623-​022-​01576-8

Okino K, Shimakawa Y, Nagaoka S (1994) Evolution of the Shikoku Basin. J Geo‑
magnetism Geoelectr 46(6):463–479. https://​doi.​org/​10.​5636/​jgg.​46.​463
Ortega-Culaciati F, Simons M, Ruiz J, Rivera L, Díaz-Salazar N (2021) An EPIC
Tikhonov regularization: application to quasi-static fault slip inversion. J
Geophys Res: Solid Earth 126:e2020JB021141. https://​doi.​org/​10.​1029/​
2020J​B0211​41
Ozawa S (2017) Long-term slow slip events along the Nankai trough subduc‑tion zone after the 2011 Tohoku earthquake in Japan. Earth, Planets Sp
69:56. https://​doi.​org/​10.​1186/​s40623-​017-​0640-4
Ozawa S, Kawabata R, Kokado K, Yarai H (2020) Long-term slow slip
events along the Nankai trough delayed by the 2016 Kumamoto
earthquake, Japan. Earth, Planets Sp 72:61. https://​doi.​org/​10.​1186/​
s40623-​020-​01189-z
Peacock SM, Wang K (1999) Seismic consequences of warm versus cool sub‑duction metamorphism: examples from southwest and northeast japan.
Science 286(5441):937–939. https://​doi.​org/​10.​1126/​scien​ce.​286.​5441.​937
Ruina A (1983) Slip instability and state variable friction laws. J Geophys Res
Solid Earth 88:10359–10370. https://​doi.​org/​10.​1029/​JB088​iB12p​10359
Scholz CH (1998) Earthquakes and friction laws. Nature 391:37–42. https://​doi.​
org/​10.​1038/​34097
Seshimo Y, Yoshioka S (2022) Spatiotemporal slip distributions associated withthe 2018–2019 Bungo Channel long-term slow slip event inverted from
GNSS data. Sci Rep. https://​doi.​org/​10.​1038/​s41598-​021-​03982-6
Seshimo Y, Kawabata H, Yoshioka S, Ortega-Culaciati F (2023) Spatiotemporal
slip distribution associated with the 2012–2016 Tokai long-term slow
slip event inverted from GNSS data. Geosci Lett 10:61. https://​doi.​org/​10.​
1186/​s40562-​023-​00316-4
Suenaga N, Yoshioka S, Matsumoto T, Ji Y (2018) Two-dimensional thermal
modeling associated with subduction of the Philippine Sea plate insouthern Kyushu, Japan. Tectonophysics 723:288–296. https://​doi.​org/​10.​1016/j.​tecto.​2017.​12.​017
Tahara M, Uehira K, Shimizu H, Nakada M, Yamada T, Mochizuki K, Shinohara
M, Nishino M, Hino R, Yakiwara H, Miyamachi H, Umakoshi K, Goda M,
Matsuwo N, Kanazawa T (2008) Seismic velocity structure around theHyuganada region, Southwest Japan, derived from seismic tomographyusing land and OBS data and its implications for interplate coupling andvertical crustal uplift. Phys Earth Planet Inter 167:19–33. https://​doi.​org/​
10.​1016/j.​pepi.​2008.​02.​001
Takagi R, Uchida N, Obara K (2019) Along-strike variation and migration of
long-term slow slip events in the Western Nankai Subduction Zone,
Japan. J Geophys Res Solid Earth 124:3853–3880. https://​doi.​org/​10.​1029/​
2018J​B0167​38
Takamatsu N, Muramatsu H, Abe S, Hatanaka Y, Furuya T, Kakiage Y, Ohashi K,
Kato C, Ohno K, Kawamoto S (2023) New GEONET analysis strategy at GSI:
daily coordinates of over 1300 GNSS CORS in Japan throughout the lastquarter century. Earth Planets Space 75(1):49. https://​doi.​org/​10.​1186/​
s40623-​023-​01787-7
Wessel P, Smith WHF (1998) New, improved version of generic mapping tools
released. Eos Trans Am Geophys Union 79:579. https://​doi.​org/​10.​1029/​
98eo0​0426
Yabuki T, Matsu’ura M (1992) Geodetic data inversion using a Bayesian
information criterion for spatial distribution of fault slip. Geophys J Int
109:363–375. https://​doi.​org/​10.​1111/j.​1365-​246X.​1992.​tb001​02.x
Yagi Y, Kikuchi M, Sagiya T (2001) Co-seismic slip, post-seismic slip, and after‑
shocks associated with two large earthquakes in 1996 in Hyuga-nada,
Japan. Earth Planets Sp 53:793–803. https://​doi.​org/​10.​1186/​BF033​51677
Yamamoto Y, Obana K, Takahashi T, Nakanishi A, Kodaira S, Kaneda Y (2013)
Imaging of the subducted Kyushu-Palau Ridge in the Hyuga-nada region,
western Nankai Trough subduction zone. Tectonophysics 589:90–102.
https://​doi.​org/​10.​1016/j.​tecto.​2012.​12.​028
Yamashita Y, Yakiwara H, Asano Y et al (2015) Migrating tremor off southern
Kyushu as evidence for slow slip of a shallow subduction interface. Sci‑
ence 348:676–679. https://​doi.​org/​10.​1126/​scien​ce.​aaa42​42
Yamashita Y, Shinohara M, Yamada T (2021) Shallow tectonic tremor activities
in Hyuga-nada, Nankai subduction zone, based on long-term broadband
ocean bottom seismic observations. Earth Planets Sp 73:196. https://​doi.​
org/​10.​1186/​s40623-​021-​01533-x
Yarai H, Ozawa S (2013) Quasi-periodic slow slip events in the afterslip areaof the 1996 Hyuga-nada earthquakes, Japan. J Geophys Res Solid Earth
118:2512–2527. https://​doi.​org/​10.​1002/​jgrb.​50161

Kawabata et al. Progress in Earth and Planetary Science(2024) 11:24

Yoshioka S, Mikumo T, Kostoglodov V, Larson KM, Lowry AR, Singh SK (2004)
Interplate coupling and a recent aseismic slow slip event in the Guerrero
seismic gap of the Mexican subduction zone, as deduced from GPS data
inversion using a Bayesian information criterion. Phys Earth Planet Inter
146:513–530. https://​doi.​org/​10.​1016/j.​pepi.​2004.​05.​006
Yoshioka S, Matsuoka Y, Ide S (2015) Spatiotemporal slip distributions of three
long-term slow slip events beneath the Bungo Channel, southwest
Japan, inferred from inversion analyses of GPS data. Geophys J Int
201:1437–1455. https://​doi.​org/​10.​1093/​gji/​ggv022

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in pub‑lished maps and institutional affiliations.

Page 17 of 17