The Yaeyama Islands are located in the extreme southwestern part of Japan, situated on the Ryukyu Arc (Fig. 1a). This is the zone where the Philippine Sea plate (PH) is subducting northwestward along the Ryukyu Trench at a rate of 8.0–8.5 cm/yr (Sella et al. 2002). In addition, back-arc spreading occurs along the Okinawa Trough on the northern part of the islands (Sibuet et al. 1998) at a rate of 3.5–5.0 cm/yr (Nishimura et al. 2004). Therefore, the plate convergence rate of 12.0–13.0 cm/yr in this region, relative to the PH, is considered to be extremely fast.
Along the subducting PH interface, various types of fault slips with different time constants have been documented. Based on the Global Navigation Satellite System (GNSS) observations, slow slip events (SSEs), a slow fault slip that lasts for approximately one month, have been recurrently located at ~ 25 to 40 km depths beneath the northwestern part of Iriomote Island (Fig. 1a) (Heki and Kataoka 2008; Nishimura 2014; Tu and Heki 2017). These SSEs occur approximately every six months with a moment magnitude (Mw) of 6.6–6.7 (Kano et al. 2018b). The large slip areas of SSEs exhibit an average slip rate of ~ 13 cm/yr. These are comparable to the relative plate motion, indicating that the main slip area is entirely coupled during the inter-SSE period and the SSEs release most of the accumulated strain (Kano et al. 2018b). On the shallower side of the SSE region at depths of ~ 10 to 20 km, the seismic signatures of slow fault slips with dominant frequencies of a few Hz and a few tens of seconds, known as low-frequency earthquakes (LFEs) and very low-frequency earthquakes (VLFEs), have often been identified (Ando et al. 2012; Nakamura and Sunagawa 2015; Arai et al. 2016; Nakamura 2017). For the shallowest part of the plate interface, Nakamura (2009) indicated the possibility of occurrence of the 1771 Yaeyama earthquake close to the Ryukyu Trench (Fig. 1a). He demonstrated that the shallow thrust-type fault model with Mw of 8.0 can reproduce the tsunami height records observed at the Yaeyama Islands. In addition, tsunami deposits that recorded four tsunami events with a recurrence interval of ~ 600 years supported this shallow fault model. Based on the slip amount of the source model in Nakamura (2009), the plate convergence rate, and recurrence interval of earthquakes, the seismic coupling ratio in the region was estimated to be ~ 20% (Ando et al. 2018). Moreover, geodetic analysis indicated low interplate coupling near the trench region. Watanabe and Tabei (2004) analyzed GNSS velocity data for 1996–1999. Due to sparse GNSS stations, they combined GNSS data with the moment tensor data of shallow earthquakes to increase the spatial resolution at shallow depths and suggested that interplate coupling close to the trench in the analysis period was < 10%. Nevertheless, the source region of the 1771 Yaeyama earthquake is still debatable, with proposed mechanisms defining it as a megathrust tsunamigenic earthquake (Nakamura 2009), an intraplate earthquake with a submarine landslide (Imamura et al. 2008), an earthquake along a splay fault (Hsu et al. 2013), and the collapse of an accretionary prism near the trench (Okamura et al. 2018), as summarized in Fujiwara et al. (2020).
For further understanding of the interplate coupling along the entire southern Ryukyu subduction zone from shallow to deep depths, additional GNSS stations were installed in the 2010s on the Yaeyama Islands (Fig. 1a). Stable estimation of secular velocities and distance changes between stations is currently possible due to the continuous acquisition of data. Therefore, this study estimates the spatial distribution of the slip deficit rate (SDR) and discusses the possibility of a megathrust event in the southern Ryukyu region by analyzing the secular velocities and distance changes.
Secular velocities and distance changes between GNSS stations on the Yaeyama Islands
The GNSS daily coordinates were analyzed to estimate the secular velocities at 14 GNSS stations deployed on the Yaeyama Islands (Fig. 1a). The GNSS stations consist of four permanent GNSS Earth Observation Network System (GEONET) stations (0500, 0749, 0750, and 0751) operated by the Geospatial Authority of Japan (GSI), six stations (FNUK, INDA, KOHM, KRSM, OOHR, and SRH1) operated by Kyoto University (KT), three stations (YKNK, YOGN, and YPNR) operated by Kyushu University (KS), and one station (ISGK) operated by the National Astronomical Observatory of Japan (NAO). All GEONET stations have been operational since the 1990s, and the NAO station was installed in 2002. In addition to these stations, four of the KT stations (FNUK, KOHM, KRSM, and OOHR) were established in 2010, which contributed to improving the spatial resolution for monitoring SSEs beneath the islands (Kano et al. 2018b). The remaining KT and KS stations were established in 2017 or 2018. The GNSS data from January 2010 to February 2021 were pre-processed by Gipsy-X ver1.3/1.4 software to obtain the time series of daily coordinates. The VMF1 mapping function (Boehm et al. 2006) and FES2014b model (Lyard et al. 2021) were used for tropospheric delay modeling and correction of ocean tide loading, respectively. The precise ephemerides and translational parameters of coordinates into the International Terrestrial Reference Frame (ITRF) 2014 coordinate system were provided by the Jet Propulsion Laboratory. Finally, the coordinates were converted to east-west (EW), north-south (NS), and up-down (UD) components (Fig. 1b–1d).
The secular velocity of each component at each GNSS station was estimated using est_noise software (Langbein 2017), which assumes a temporal correlation for modeling the observation noise. Here, the noise model was adopted in combination with white noise, flicker noise, random-walk noise, and band-pass filtered noise of 0.5–2 cycles/yr with three poles, to maximize a likelihood function (Langbein 2004). The analysis modeled offsets due to antenna replacements or relocation of the antenna to the nearby building, along with larger (> 6.0) magnitude in the unified hypocenter catalog provided by the Japan Meteorological Agency (JMA) that occurred close to the Yaeyama Islands (Table 1). The obtained secular velocities in the horizontal direction exhibited a south-southeastward motion of 63.9–80.2 mm/yr with a maximum error of ~ 3.1 mm/yr in the ITRF 2014 (Fig. 2a). On the other hand, vertical velocities were estimated to range approximately from − 2.4 to 1.1 mm/yr with errors ranging from ~ 0.5 to 2.5 mm/yr, thereby exhibiting no clear uplift or subsidence (Fig. 2b).
The obtained secular velocities were assumed to be composed of the combined effect of the rigid motion of the South Ryukyu (SR) block where the Yaeyama Islands are located and that of the interplate coupling if it exists. To eliminate the effect of the rigid block motion, the horizontal secular velocities were converted to the distance changes between each pair of the stations. As a result, the distance changes ranged 2.4 mm/yr in maximum contraction and 4.7 mm/yr in maximum extension. The distances in the direction from east-west to northwest-southeast show contraction systematically and those in the direction from northeast-southwest shown extension (Fig. 2c). The former is roughly parallel to the direction of plate convergence and implies the interplate coupling. We used distance changes and vertical secular velocities in the following inversion analysis.