Potential of megathrust earthquakes along the southern Ryukyu Trench inferred from GNSS data

The southern part of the Ryukyu subduction zone has recorded tsunami events with a recurrence interval of several hundred years. Although their source is controversial, one model suggests that the last 1771 Yaeyama tsunami was caused by a shallow megathrust earthquake with a magnitude of 8. However, the current knowledge on interplate coupling based on recent geodetic data is limited. Here, a time series of Global Navigation Satellite System data from January 2010 to February 2021 was analyzed, including newly installed stations by Kyoto and Kyushu Universities, to obtain the distance changes between stations and vertical secular velocities. The distance changes ranged from 2.4 mm/year in contraction and to 4.7 mm/year in extension, and the vertical velocities exhibited no clear uplift or subsidence, with − 2.4 to 1.1 mm/year. The back slip inversion results indicated a slip deficit of 17–47 mm/year to the south of the Yaeyama Islands. The large slip deficit area is complementarily intervened between the shallower source area of low-frequency earthquakes and the deeper slow slip region, suggesting the spatial heterogeneity of frictional properties along the plate interface. If the large slip deficit area accumulates stress in the same rate since the last 1771 earthquake, it could result in a megathrust event with a moment magnitude greater than 7.5. Because the limited onshore data cannot resolve the slip deficit on the shallow plate interface, seafloor geodetic observations are essential to clarify the detailed spatial distribution of the slip deficit and discuss its earthquake and tsunami potential.


Introduction
The Yaeyama Islands (which is the generic collective term for this group of 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/year relative to the South China block (Sella et al. 2002). In addition, back-arc spreading occurs along the Okinawa Trough to the north of the islands (Sibuet et al. 1998) at a rate of 3.5-5.0 cm/ year (Nishimura et al. 2004). Therefore, the plate convergence rate of 12.0-13.0 cm/year in this region, relative to the PH, is considered to be fast.
Along the subducting PH interface, various types of fault slips with different time constants have been documented. Based on Global Navigation Satellite System (GNSS) observations, slow slip events (SSEs), slow fault slips that last for approximately 1 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 6 months with a moment magnitude (Mw) of 6.6-6.7 (Kano et al. 2018b). Kano et al. (2018b) investigated the slip history at the plate interface for ~ 2.5 years, including five SSEs. They estimated an average slip rate of ~ 13 cm/year at the large slip area of SSEs by dividing the total slip for the five SSEs by the analysis period, including both SSE and inter-SSE periods. This average slip rate is 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 (Fig. 1a) (Ando et al. 2012;Nakamura and Sunagawa 2015;Arai et al. 2016;Nakamura 2017). The magnitudes of LFEs and VLFEs were estimated to be 2.8-3.4 and 3.4-4.8, respectively (Nakamura and Sunagawa 2015;Nakamura 2017). For the shallowest part of the plate interface, Nakamura (2009) indicated the possibility of the occurrence of the 1771 Yaeyama earthquake close to the Ryukyu Trench (Fig. 1a). He demonstrated that a shallow thrust-type fault model with an Mw of 8.0 could 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 (Ando et al. 2018). Based on the slip amount of the source model in Nakamura (2009), the plate convergence rate, and the 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 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 and related tsunami is still under debate, with proposed mechanisms of a megathrust tsunamigenic earthquake (Nakamura 2009), an intraplate earthquake with a submarine Fig. 1 Tectonic setting and GNSS time series. a Study area. The black rectangle in the inset is the location of the Yaeyama Islands where the Philippine Sea plate (PH) is subducting beneath the islands at a rate of 8.0-8.5 cm/year relative to the South China block (SC) (Sella et al. 2002). The Yaeyama Islands are considered to be located on the South Ryukyu block (SR), which are roughly indicated by the dotted blue rectangle (Nishimura et al. 2004). The squares are the locations of the GNSS stations with four-digit station codes. The purple squares indicate the stations newly installed after 2010. The red-, blue-, and orange-colored regions correspond to the source areas of slow slip events (SSEs) (Heki and Kataoka 2008), low-frequency earthquakes (LFEs) (Arai et al. 2016;Nakamura 2017), and the estimated source area of the 1771 Yaeyama earthquake (Nakamura 2009). The white arrow represents the motion of the back-arc spreading of the Okinawa Trough (OT) (Nishimura et al. 2004). The location of OT is indicated in the inset. The solid and dashed lines are the Ryukyu Trench and the depth contours of the upper surface of the subducting PH with an interval of 10 km (Slab 2) (Hayes et al. 2018). b-d GNSS time series sorted from westernmost (top) to easternmost (bottom) stations in b east-west, c north-south, and d up-down components. The red and blue vertical lines represent the timings of earthquakes considered in this study and the occurrence of SSEs. The timings of SSEs were determined by Tu and Heki (2017) before 2016 and by visual inspection after 2017. The orange arrows indicate the timings of antenna replacements, and green arrows indicate an unknown offset appearing in the up-down component, which was eliminated in estimating velocities. GNSS Global Navigation Satellite System 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 southern Ryukyu subduction zone, 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.

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 Information 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 five KT and KS stations were established in 2017 or 2018. Through the installation of these new KT and KS stations (highlighted by purple squares in Fig. 1a), we expected to improve the spatial resolution for monitoring the crustal deformation on the Yaeyama Islands, especially across the western parts. The GNSS data from January 2010 to February 2021 were preprocessed by the 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-d). We used all data from the period from January 2010 to February 2021 at each site to estimate the secular velocities in the following analysis.
The secular velocity of each component at each GNSS station was estimated using the est_noise software (Langbein 2017), which assumes a temporal correlation for modeling observation noise. Here, the noise model was adopted in combination with white noise, flicker noise, random-walk noise, and bandpass-filtered noise of 0.5-2 cycles/year with three poles, to maximize a likelihood function (Langbein 2004). The analysis modeled offsets due to antenna replacements, relocation of the antenna to the nearby building, or earthquakes, and apparent offsets with unclear origin (Fig. 1). Here, we considered earthquakes 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, Fig. 1 and Additional file 1: Fig. S1). The obtained secular velocities in the horizontal direction exhibited a south-southeastward motion of 63.9-80.2 mm/year with a maximum error of ~ 3.1 mm/ year in the ITRF 2014 ( Fig. 2a). On the other hand, vertical velocities were estimated to range approximately from − 2.4 to 1.1 mm/year with errors ranging from ~ 0.5 to 2.5 mm/year, thereby exhibiting no clear uplift or subsidence (Fig. 2c). We also estimated the secular velocities from January 2015 to February 2021 (Additional file 1:  Fig. S2), resulting in a difference in secular velocities of < 1.4 mm/year and < 2.8 mm/year in horizontal and vertical components, respectively. However, this difference did not affect the location of the large SDR region and the amplitudes shown in the following sections ( Fig. 3 and Additional file 1: Fig. S3). Since this study focuses on long-term average crustal movements, we did not model the occurrence of SSEs (although signals of SSEs were repeatedly observed). This means that the secular velocities estimated here include the effects of stress release during the SSE period as well as stress accumulation during the inter-SSE period. 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 interplate coupling, if existing. To eliminate the effect of rigid block motion, the horizontal secular velocities were converted to the distance changes between each pair of stations. As a result, the distance changes ranged from 2.4 mm/year for maximum contraction to 4.7 mm/year for maximum extension. The distances in the direction from east-west to northwest-southeast showed systematical contraction, and those in the direction from northeast-southwest showed extension (Fig. 2d). The former was roughly parallel to the direction of plate convergence and implied interplate coupling. We used distance changes and vertical secular velocities in the following inversion analysis.

Inversion analysis
Assuming that the distance changes and the vertical secular velocities can be caused by interplate coupling, we conducted inversion to estimate SDR by implementing the backslip model (Savage 1983) in the following procedure. The plate interface of the subducting PH was modeled using the plate geometry of Slab 2 (Hayes et al. 2018) with planar subfaults (Fig. 3a). The entire region was divided into 9 × 10 subfaults in the strike and dip directions, respectively. The size of each subfault was 20 × 20 km 2 when the subfault was projected on the surface. We tested additional cases with a fault discretization of 15 × 15 km 2 and 10 × 10 km 2 subfaults and found that this fault discretization did not affect the main results (Additional file 1: Fig. S4). The surface response due to the fault slip rate in each subfault was calculated assuming a homogeneous elastic half-space (Okada 1992). The slip direction in each subfault was fixed to N45W, which is the approximately parallel direction to that of the subduction of the PH relative to the SR (Sella et al. 2002;Nishimura et al. 2004). We constrained non-negative fault slips and assumed zero fault slip in the additional subfaults surrounding the entire fault region, which are not explicitly represented in each figure, except for Fig. 3a. The smoothness constraint of the SDR in the neighboring subfaults was introduced to stabilize the inversion. Inversion was conducted to minimize the following cost function s(m): where d is the N-dimensional observed vector containing the distance changes and the vertical velocities, G is a matrix consisting of an elastic response on the ground surface, m is the M-dimensional model parameter vector including the amounts of SDR, E is the observation error covariance matrix, and L indicates a Laplacian matrix having a rank of P. In this calculation, the cost function s(m) monotonically decreases with an increasing value of hyperparameter α 2 in our problem, and thus, the hyperparameter ( α 2 = 200) is adopted when the changing rate of a value in the cost function is small (Additional file 1: Fig. S5). The estimation error of the SDR is calculated from the corresponding diagonal element of the error covariance matrix C using the optimum parameter vector m. (1)

Spatial distribution of slip deficit rate
The inversion results indicated a SDR distribution in a wide area of the fault region, with a value of up to 47 mm/year to the south of Iriomote Island (Fig. 3a). The SDR distributions reproduced the observed distance changes and the vertical secular velocities well ( Fig. 3d-g). However, their estimation errors were large as well, ranging from 17 to 45 mm/year (Fig. 3b), which were comparable or larger than the SDR, especially in the deep or eastern subfaults as well as the shallow subfaults. Therefore, hereafter, we discuss only the subfaults that had SDRs larger than their errors (Fig. 3c). The SDRs in these subfaults contributed to 56% of the residual reduction compared to the value calculated using the whole subfaults. Vertical velocities and distance changes calculated due to these subfaults (Additional file 1: Fig. S6) showed a similar spatial pattern compared to those due to the whole subfaults (Fig. 3d, f ). However, there were differences between these calculated vertical displacements at stations located on the eastern and western side of the network, resulting in a worse fitting to the observation if we only considered selected subfaults. Similarly, the distance changes showing northwest-southeast contraction resulted in overestimation. This indicates that SDRs of deep subfaults smaller than the estimation errors would be necessary to further explain the data. However, it is difficult to show a significant existence of deep SDRs from the limited observations. To investigate the robustness of the results, we conducted five synthetic tests assuming the SDRs in a checkerboard pattern (Additional file 1: Fig. S7), in a single region where large SDRs were estimated (Additional file 1: Fig. S8), or in shallow subfaults (Additional file 1: Fig. S9), as summarized in Additional file 1: Text S1. The checkerboard test indicated that SDRs deeper than ~ 20 km were well resolved, while those at shallower depths (< 15 km) were not (Additional file 1: Fig. S7). The locations of deep large SDR subfaults were not biased (Additional file 1: Fig. S8), although the smoothing constraint resulted in smoothed SDR distributions. On the other hand, the SDRs assumed at shallower depths were estimated to be deeper than their true location (Additional file 1: Fig. S9). In addition, the estimated SDRs were far too smoothed to show the existence of shallow SDRs. These synthetic tests indicated that the SDRs deeper than ~ 20 km were reliably estimated; however, we cannot discuss the SDRs at shallow depths.
Most of the subfaults, with SDRs larger than their errors, were located slightly to the south of the islands at depths of ~ 20 to 30 km (Fig. 3c). These SDRs were Fig. 3 Estimation of slip deficit rates (SDR) and comparison with observations. a Spatial distribution of SDRs and b their estimation error. c Same as a but only for subfaults whose SDR is greater than the estimation error shown. Note that the additional subfaults imposing zero SDR are shown only in a. The solid and dashed lines are the Ryukyu Trench and the depth contours of the upper surface of the subducting PH with an interval of 10 km (Slab 2) (Hayes et al. 2018). d Calculated distance changes and e residual distance changes. f Calculated and g residual vertical velocities considered to be reliable according to the synthetic tests (Additional file 1: Figs. S7, S8, Text S1). In contrast, the SDRs in the shallow subfaults were rarely resolved (Additional file 1: Fig. S9), so that limited onshore GNSS solely cannot clarify the existence of the SDRs in the shallow plate interface where Nakamura (2009) determined the source model of the 1771 Yaeyama earthquake (Fig. 1a). We confirmed that these characteristics of the spatial distributions of SDRs were common irrespective of the choice of the analysis period (Additional file 1: Fig. S2), fault discretization (Additional file 1: Fig. S3), and selection of data (Additional file 1: Fig. S10, Text S2) used in the inversion. Figure 4 compares the spatial distributions of SDRs with the epicenters of LFEs (Arai et al. 2016;Nakamura 2017) and the slip rates of cumulative SSEs from July 2010 to February 2013 (Kano et al. 2018b). Notably, our analysis suggests the existence of a large SDR region at depths of 20-25 km that is complementarily distributed with LFE epicenters on the shallower side and the large slip area of SSE on the deeper side. This may reflect that frictional properties vary with depth along the plate interface, although it is uncertain whether the shallow plate interface where the LFEs occur is coupled or not due to the low resolution (Additional file 1: Fig. S1b). The depth boundary between large SDR and SSE areas was ~ 25 km, which corresponds to the structural boundary between the crust and mantle wedge in the overriding plate, as investigated by a seismic reflection survey (Figure 3 in Arai et al. 2016); the large SDR areas are located below the arc crust, while the large SSE areas are located below the mantle wedge. Although the profile in the seismic survey (indicated by the purple line in Fig. 4) was located ~ 50 km east of Ishigaki Island, the plate geometry does not change significantly in the along-strike direction, and thus, the fault slip style along the plate interface may be related to the structure in the overriding plate. A similar relation between the structure in the overriding plate and fault slip style was indicated by seismic tomography (Yamamoto et al. 2018). Yamamoto et al. (2018) discussed the high Vp/Vs of oceanic crust inferred below the SSE region, indicating the existence of fluid that contributes to SSE activities. In addition, their result may suggest an along-dip variation of Vp/Vs, that is, lower Vp/Vs values near the large SDR region compared to the SSE region. However, considering the spatial resolution of tomography (~ 30 km in the horizontal direction), further investigation will be essential to discuss an along-dip variation of the velocity structure and its relation to fault slip style. Figures 3a and 4 show small SDRs in the main slip area of the SSEs. As already mentioned, the secular velocities are the long-term average crustal movements that include the effects of both stress release and accumulation in the SSE region. Therefore, small SDRs in the SSE region indicate that the accumulated stress during the inter-SSE period are fully released by the repeated SSEs, as suggested by Kano et al. (2018b).

Discussion
Our results reliably indicate SDRs of 17-47 mm/year to the south of the Yaeyama Islands. However, if we consider a 1 − σ estimation error, the SDRs on the significant patches range from 0.09 to 71 mm/year. In the following discussion, to roughly investigate the magnitude of possible fault slips, we used the representative SDR value of 33 ± 24 mm/year calculated by averaging the SDRs for the subfaults to the south of the Yaeyama Islands (the black rectangle in Fig. 4). This value corresponds to 26 ± 19% of the relative plate motion (~ 13 cm/year) of the PH relative to the SR. This result indicates that SDRs exist despite the low coupling ratio, suggesting the accumulation of elastic strain for future possible earthquakes.
One possibility to release the slip deficit could be a seismic slip. As there are no records of magnitude 8 class earthquakes around the Yaeyama Islands for the 250 years following 1771 (Ando et al. 2018), it can be stated that if this region accumulates a slip deficit at the same rate in the 2010s analyzed in this study as for the Fig. 4 Summary of spatial distributions of various fault slips in the southern Ryukyu subduction zone. The red-colored squares indicate the SDR estimated in this study (Fig. 3c). The blue contours are the average slip rate in the SSE region, which were calculated from the cumulative slip from July 2010 to February 2013 (Kano et al. 2018b) with a contour interval of 40 mm/year. The orange dots represent the epicenters of LFEs (Arai et al. 2016;Nakamura 2017). Note that only the relocated LFEs in Nakamura (2017) using the S-P time delays, as well as the cross-correlation of the envelope waveforms, are plotted. The orange-and yellow-colored squares are the source models of the 1771 Yaeyama earthquake estimated by Nakamura (2009) and Imamura et al. (2001), respectively. The purple line indicates the location of the seismic survey (Arai et al. 2016) 250-year period, the expected slip could be 8.4 ± 6.1 m. If we assume a rigidity of 30 GPa and the source region indicated by the black rectangle in Fig. 4 with a uniform slip of 8.4 m (we refer to this fault model as model I), this leads to an earthquake with a seismic moment of 8.1 ± 5.9 × 10 20 Nm corresponding to an Mw of 7.5-8.0. As mentioned before, the resolution of the subfaults shallower than the large SDRs is low (Additional file 1: Fig. S1b), thus it may be possible that the source region extends to shallower depths. If we extend the source region (~ 120 km wide) near the trench axis (we refer to this fault model as model II as indicated in Fig. 4) to reach the location where Nakamura (2009) estimated the source model of the 1771 Yaeyema earthquake, the Mw becomes larger, up to ~ 8.3, assuming the same uniform slip as in model I. Usami (2010) reported from historical documents that Ishigaki Island experienced a ground shaking corresponding to the JMA seismic intensity of IV in the 1771 Yaeyama earthquake as summarized in Ando et al. (2018). Supposing that either models I or II caused this earthquake, both models could result in the large ground shaking in the Ishigaki Island as reported by Usami (2010). Tsunami deposits indicated that the 1771 earthquake generated large tsunamis leading to maximum runups of ~ 27 m on the eastern coast of Ishigaki Island (e.g., Ando et al. 2018). As mentioned above, model II includes a part of the source region of Nakamura (2009), which was determined to explain the recorded tsunami heights by numerical simulation. Therefore, model II can generate tsunamis as well although the amplitude of shallow slip needs to be further considered for the quantitative comparison of tsunami run-up heights. On the other hand, model I is unlikely to generate large tsunamis because it is an interplate earthquake at depths of 20-25 km. Thus, a possible scenario to cause such a large tsunami may be that a submarine landslide occurred in addition to the earthquake itself. Imamura et al. (2001) proposed a model in which the submarine landslide triggered by an intraplate earthquake caused the large 1771 tsunami. Although the location of this source model is ~ 40 km away from the eastern edge of the large SDR area (indicated by a yellow square in Fig. 4), if the interplate earthquake in model I resulted in a submarine landslide in the same location, the large tsunami run-up height could be explained. Other models of the 1771 tsunami were proposed (Hsu et al. 2013;Okamura et al. 2018). Based on seismic reflection surveys, Hsu et al. (2013) proposed that the tsunami was related to a thrust fault along the ~ 450 km long mega-splay fault in the trench parallel direction east of 125.5° E close to the Ryukyu Trench. On the other hand, Okamura et al. (2018) discovered a seafloor depression on the seafloor located above the source model of Nakamura (2009). This depression was interpreted to be caused by the collapse of the accretionary prism and the resulting rotational slide. They showed that the collapse can explain the 1771 tsunami regardless of the coupling state of the shallow plate interface. Both models indicated a tsunami source close to the Ryukyu Trench; however, our results cannot judge these possible scenarios, as we only considered the SDRs along the subducting plate interface, and onland GNSS data cannot resolve the coupling of splay faults located close to the trench.
Another possibility through which SDRs are released is by aseismic slip, that is, through SSEs. Although SSEs were recurrently identified with an Mw greater than 6.2, particularly on the deeper side of the large SDR areas (Heki and Kataoka 2008;Nishimura 2014), SSEs including the large SDR region as the source area have not been detected at least since the establishment of GEONET. It may be possible that such a transient slip would occur in the edge of the large SDR region as inferred in the Nankai subduction zone (Kano et al. 2019). Considering the high resolution in the large SDR region, the GNSS observations would be able to detect similar magnitudes of SSEs as those of the deeper ones if an SSE were to occur in the large SDR region. There are three possibilities. One possibility is that quite a large SSE would occur in this region, which has, however, not occurred at least in recent two decades. Supposing that such a large SSE would occur with a recurrence interval of 30 years and release all the slip deficits of 33 ± 24 mm/year, the SSE will have a seismic moment of 9.5 ± 6.9 × 10 19 Nm, corresponding to an Mw of 6.9-7.4, which would be certainly detected by the GNSS. In this case, the SDRs can be explained by aseismic slip only. However, large SSEs with such a long recurrence interval have not been observed because of a lack of long-term GNSS observations. Therefore, this possibility is less plausible, although we cannot entirely rule it out. The second possibility is that a long-lived SSE with durations longer than a few decades is occurring; however, this would be difficult to identify from recent GNSS observations. The third possibility is that SSEs would have occurred within the large SDR region, but were too small to be detected by the GNSS. However, it would be difficult to explain the large SDR based on such a series of relatively small SSEs. In this case, the large slip deficit would not be released by an aseismic slip alone but seismically. Although the observations, which are limited in both time and space, cannot resolve which of the three possibilities discussed above is true, continuous GNSS observations as well as seafloor geodetic observations will reveal more about the potential of megathrust earthquakes in the southern Ryukyu region.
Our inversion results indicated significant SDR regions of 17-47 mm/year that would be possibly released by seismic slip. These SDR values correspond to the low coupling ratio of 13-36% calculated from the relative plate motion (~ 13 cm/year) of the PH relative to the SR. Although these values include estimation errors, they suggest that these regions may mostly show a steady sliding, but occasionally result in a dynamic slip. Numerical simulations on earthquake cycles suggested that constant sliding and dynamic slip will occur in the same fault during complex cycles (e.g., Hori and Miyazaki 2011;Nakata et al. 2012;Noda and Lapusta 2013). For example, Hori and Miyazaki (2011) reported that aseismic stable sliding with a speed of ~ 34-62% of the plate subduction rate occurred after a nearby magnitude 7 class earthquake, while the same area indicated a dynamic rupture during a magnitude 9 class earthquake. This complex slip behavior may not be the result of a single source fault but that of a complex interaction between nearby multiple source regions. Although this model is not directly applicable to our study area, coexistence of both fast slip and steady sliding within a low coupling area may imply the existence of another potential seismic source located, for example, in the shallow plate interface or at similar depths on the western or eastern side of the large SDR region revealed in this study. The latter region may be too far from the onshore GNSS stations to investigate plate coupling. Seafloor geodetic observations will be important from this point of view as well.
A possible modification of the proposed SDR model would involve the viscoelastic effect in the asthenosphere, while our modeling is based on the elastic model alone. Ignoring the viscoelastic effect was shown to lead to the overestimation of the lower limit of the locking depth and the amount of SDR (e.g., Li et al. 2015Li et al. , 2018Itoh et al. 2021). This is because the viscous relaxation lends itself to a long-wavelength deformation, resulting in movement at more distant stations compared to the elastic model. These results were based on inversion results using secular velocities, while we used distance changes in addition to vertical velocities. Therefore, further analysis is necessary to examine the effect of viscoelasticity on our results. Introducing viscoelasticity is also important while interpreting the temporal variation of vertical velocities. When recurrence intervals of earthquakes are considerably longer than the asthenospheric viscous relaxation time, as in the case of the southern Ryukyu area, decadal-scale temporal variation of interseismic deformation becomes more significant in the later stage of the interseismic period (Sagiya 2015;Hashima and Sato 2017). Hashima and Sato (2017) numerically calculated the vertical deformation associated with the 2011 Tohoku-oki earthquake using a simple two-layered lithosphere-asthenosphere model. Their results indicated that the viscoelastic effect due to coseismic slip is more significant than the effect of interplate locking in the early stage of the interseismic period, while the latter is predominant in the later interseismic period. As a result, temporal variation of vertical velocities at onland stations > 50 km from the interplate coupling area became more significant. In addition, the spatial pattern of the vertical deformation was different from the assumed slip models whether the coseismic slip reaches the asthenosphere (Hashima and Sato 2017). Therefore, how viscoelasticity affects SDR estimation depends on the assumed model. In any case, continuous monitoring of vertical velocities and their decadal-scale temporal change will lead to a more realistic modeling of interplate coupling with a viscoelastic model. In addition, although we only consider the motion due to plate subduction, viscoelastic effects due to back-arc spreading to the north of the Yaeyama Islands will affect the secular velocities. Therefore, this should be taken into account when considering viscoelasticity, which will be subject to future study.

Conclusion
We estimated the spatial distributions of SDRs in the southern Ryukyu subduction zone based on GNSS observations during the last 10 years. The inversion results indicated a reliable estimate of 17-47 mm/year SDRs to the south of the Yaeyama Islands. The large SDR areas are located at depths of 20-25 km, intervened by LFE epicenters on the shallower side and SSE regions on the deeper side of the plate interface. If this region accumulates stress in the same rate as in the recent decade since the last 1771 Yaeyama earthquake, this may result in an Mw of 7.5 or larger megathrust earthquakes. Although limited onshore observations cannot determine whether the large SDR regions extend to the shallow fault region near the trench axis, as estimated by Nakamura (2009), and the true tsunami source in 1771, current geodetic data indicate that the southern side of the Yaeyama Islands has the potential for megathrust earthquakes. Seafloor geodetic observations as well as continuous onland GNSS observations will be essential to reveal more about the earthquake potential in the southern Ryukyu region including the detailed spatial distribution of SDRs especially at shallower depths.