Slow slip events (SSEs) usually occur around a locked seismogenic zone in the updip, downdip, or at similar depths where the plate interface is weakly coupled (e.g., Schwartz and Rokosky, 2007; Avouac, 2015; Bürgmann, 2018). The regions where SSEs occur have time-dependent transitional friction from velocity weakening to velocity strengthening probably caused by high-pressure interface fluids (e.g., Saffer and Wallace, 2015; Obara and Kato, 2016; Behr and Bürgmann, 2021). Recent observations have shown that SSEs in the downdip of the seismogenic zone initiate immediately before seismogenic-zone earthquakes. e.g., the 2011 M_w 9.0 Tohoku earthquake (Kato et al., 2012; Ito et al., 2013), the 2014 M_w 8.1 Iquique earthquake (Ruiz et al., 2014; Kato et al., 2016; Socquet et al., 2017), the 2014 M_w 7.3 Papanoa earthquake (Radiguet et al., 2016), and the 2018 M_w 6.9 Zakynthos earthquake (Saltogianni et al., 2021). Scientists have interpreted the interaction between aseismic and seismic slip as the stress perturbations caused by the SSEs changing the state of stress in the seismogenic zone beyond the triggering stress threshold of those earthquakes. Generally, the potential of triggering such earthquakes depends on the frictional and stressed conditions at the fault interface (e.g., Mazzotti and Adams, 2004; Beeler et al., 2014; Kaneko et al., 2018). Mazzotti and Adams (2004) suggested that a large earthquake occurs immediately when the state of stress in the seismogenic zone exceeds the earthquake-triggering stress threshold caused by an SSE. Beeler et al. (2014) proposed that a large earthquake is postponed indefinitely after the stress perturbation because the seismogenic zone lacks an earthquake-triggering stress threshold, perhaps due to the earthquake nucleation process. The two hypotheses are opposite depending on whether the earthquake-triggering stress threshold represents the end of an earthquake cycle. The triggering stress threshold can explain why a large earthquake occurs after one local SSE in weeks to a few months during the SSE activity (e.g., Kato et al., 2012; Ito et al., 2013; Ruiz et al., 2014; Radiguet et al., 2016). It means the seismogenic zone is very close to the end of the earthquake cycle, and thus the earthquake can be triggered by the additional stress perturbation from the SSE. If a large earthquake occurs several months to years after the local SSE activity, the seismogenic zone is not yet very close to the end of the earthquake cycle. However, the SSE may have increased the likelihood of a large earthquake when it terminated, even if no earthquake is triggered immediately (e.g., Voss et al., 2018). In the southern Hikurangi subduction zone, New Zealand, the 2016 Kapiti SSE downdip the seismogenic zone did not trigger the large earthquake immediately. The stress perturbations caused by the SSE increased the probability of a megathrust earthquake with a maximum moment of M_w 8.6 on the seismogenic zone (Kaneko et al., 2018). The stress perturbations are at least twice that caused by the 2016 M_w 7.8 Kaikōura earthquake afterslip, which occurred approximately 100 km west of the 2016 Kapiti SSE source area. Thus, the occurrence of megathrust earthquakes linked to SSE activity cannot be ignored.

Geodetic evidence of the slip deficit rate of 8.6 cm/yr in the southernmost Ryukyu subduction zone reveals the potential of an M_w 7.5 to 8.7 earthquake in a locked region named the Ryukyu fault (Hsu et al., 2012). The occurrence of an M_w 8.15 megathrust earthquake will cause a tsunami with 7 m wave heights in the eastern Taiwan region (Sun et al., 2018). The 2002 M_w 7.1 Hualien offshore earthquake that occurred near the downdip end of the Ryukyu fault produced local destruction in the Taiwan region and a tsunami of 20 cm wave heights on Yonaguni Island (Figure 1). The Ryukyu fault extends from the Ryukyu Trench to the shallow depths of the subducting plate interface (Figure 1). The convergence rate across the Ryukyu Trench may be 125 mm/yr in the NS direction between the Philippine Sea Plate and Yonaguni Island on the Eurasian Plate (Hsu et al., 2012). The rate incorporates a back-arc rifting rate of 50 mm/yr in the NS direction from the Okinawa Trough relative to the Eurasian Plate (Nishimura et al., 2004) in a convergence rate of 80 mm/yr in the 310° direction between the two plates under Taiwan Island (Hsu et al., 2009; Chen et al., 2014; 2017). Considering the high convergence rate across the Ryukyu fault, the plate interface has not produced an earthquake with a magnitude greater than M_w 7.7 in the last 417 years, as revealed by the Taiwan historical earthquake catalog since 1,604 (Cheng and Yeh, 1989; Theunissen et al., 2010). The limitation of the earthquake size is atypical of universal observations and may imply the possibility of the seismogenic potential of a much larger event in the future. However, current knowledge of the Ryukyu fault regarding the likelihood of a megathrust event of M_w 8.7 remains unclear.

Tsunami deposits in the adjacent southernmost Ryukyu subduction zone have provided crucial evidence for tsunami floods from large historical earthquakes. The tsunami deposits in northern Taiwan have shown a possible recurrence interval of 100–400 years for the large earthquakes linked to the 1867 and 1,694 local events and the potential far-field events between 1,293 and 1,414, 1,090 and 1,235 (e.g., Konstantinou et al., 2013; Lin et al., 2014; Yu et al., 2020 and 2022). The two local events in 1867 and 1,694 (approximately M_w ∼7.0) in northern Taiwan are likely associated with normal faulting (e.g., Cheng et al., 2016; Sugawara et al., 2019; Yu et al., 2020) instead of the subduction-related thrusting in the study area. In the southern Ryukyu Islands, the tsunami deposits indicate a seismic recurrence interval of 150–400 years constrained by the 1771 and 1,625 local M ∼8.0 events (e.g., Araoka et al., 2013; Ando et al., 2018). The tsunami deposits also reveal a possible longer recurrence interval of 600–1,000 years that might correlate with the prehistoric tsunami in northern Taiwan (e.g., Ando et al., 2018; Yu et al., 2020). However, the age resolution of Taiwan and southern Ryukyu archives is insufficient for clear evidence of tsunami correlation between those events (Yu et al., 2022). So far, no tsunami deposits have been reported in the southernmost Ryukyu subduction zone to constrain the megathrust earthquake cycles.

The 1920 M_w 7.7 earthquake is the maximum event recorded in this subduction zone and is the only earthquake that satisfies the inferred magnitude over the Ryukyu fault. However, the hypocenter and the source region are poorly resolved and unclear. Engdahl et al. (1998) relocated this earthquake using teleseismic travel times at a location of 122.080° E, 23.813° N, and 35 km depth (Figure 1). Their relocation seems to imply a possible intraplate event of the Philippine Sea Plate. From a local travel-time phase perspective, the earthquake relocation is at 122.0° E, 24.0° N, and 20 km depth (Cheng and Yeh, 1989) and is more likely an interplate event (Figure 1). The relocations from the Central Weather Bureau of Taiwan and Theunissen et al. (2010) are a little northeastward and shallower than that of Cheng and Yeh (1989) (Figure 1). Theunissen et al. (2010) suggested that the 1920 M_w 7.7 earthquake was likely a rupture along the subducting plate interface with a possible splay fault in the Eurasian Plate. Although the debate regarding the origin of this large historical earthquake remains, the probability of it being an interplate event is high. The previous relocations reveal a possibility that the 1920 M_w 7.7 earthquake may have nucleated at the downdip end of the Ryukyu fault (Figure 1). If this earthquake is the last rupture of the Ryukyu fault, it is crucial to the timing of the megathrust earthquake cycle and the fault stress conditions.

Some studies have shown that SSEs can coexist at the edge of a locked fault region where traditional friction may occur (e.g., Ito et al., 2013; Dixon et al., 2014; Mallick et al., 2021). Chen et al. (2018) reported repeating SSEs in 2009 and 2015 that occurred deeper than the locked Ryukyu fault region at plate interface depths of 25–45 km. The SSE cumulative slip is likely complementary to the hypocenters of M_L ≥ 5.0 earthquakes in the last 20 years (Figure 1). Another SSE in 2005 may have occurred within the Ryukyu fault at the downdip end. The three observed SSEs in the southernmost Ryukyu subduction zone did not trigger a large earthquake at the Ryukyu fault thus far. The question centers around the role of the sequence of SSEs in the seismogenic potential of the Ryukyu fault. In the adjacent southern Ryukyu Islands, the SSEs recurred every 6 to 7 months with seismic moments of M_w 5.6 to 6.8 at plate interface depths of 30–50 km (e.g., Heki and Kataoka, 2008; Nishimura, 2014). The static stress changes caused by the SSEs are 0.2–1.0 kPa on the eastern-half Ryukyu fault, which may have activated the very low-frequency earthquakes (VLFEs) at the edge of the Ryukyu fault, likely even on the eastern-half plane (Nakamura and Sunagawa, 2015). Onshore Global Navigation Satellite System (GNSS) observations in northeastern Taiwan revealed that the cumulative energy released by each SSE is equivalent to the seismic moments of M_w 6.4 to 6.6 (Chen et al., 2018; Chen S. K. et al., 2022). The SSE moment may have uncertainties hampered by a lack of near-field observations for offshore SSEs (Figure 1).

Here we estimate the likelihood of an M_w ≥ 7.5 earthquake with a maximum moment of M_w 8.7 on the Ryukyu fault and examine how that likelihood is affected by the SSEs activity. We establish possible earthquake cycle scenarios for the locked Ryukyu fault with the earthquake magnitude-frequency relation (Gutenberg and Richter, 1944) constrained by large historical earthquakes. Then, we calculate the shear stress changes caused by the three SSEs with synthetic SSE moments to simulate the potential of a megathrust earthquake triggered by the SSE downdip of the Ryukyu fault. This study provides insights into how much the stress perturbations from SSEs can affect a locked seismogenic region. The stressed conditions of the Ryukyu fault and its possible rupture behaviors are investigated for the first time.

We estimate the probability of a megathrust earthquake along the plate interface using the method of Kaneko et al. (2018) with three steps: 1) calculating stress changes on the Ryukyu fault due to SSEs, 2) establishing a synthetic shear stress earthquake cycle of megathrust events over millions of years for the Ryukyu fault, and 3) combining 1) with 2) to estimate the probability of triggering a megathrust event. First, we calculate static stress changes on the Ryukyu fault due to the observed SSEs from the two source depths in Figure 1 with synthetic moments from M_w 6.0 to 7.0. In contrast to Kaneko et al. (2018), we use the Coulomb stress criteria (e.g., King et al., 1994; Stein, 1999; Toda et al., 2011), which have been widely used to quantify the stress changes caused by SSEs (e.g., Radiguet et al., 2016; Voss et al., 2018; Cruz-Atienza et al., 2021). We calculate Coulomb stress changes from 2005, 2009, and 2015 M_w 6.4 to 6.6 SSEs and use those stress changes to scale the relative stress changes from synthetic M_w 6.0 to 7.0 SSEs. The fault dimensions and cumulative slip of SSEs on the receiver faults are scaled proportionally with the SSE moments (Michel et al., 2019; Dal Zilio et al., 2020). The moment is proportional to A^3/2, where A is the source area of each SSE and obeys the earthquake magnitude-frequency relation. The strike, dip, and rake of each SSE are the same as the fault parameters of 2005, 2009, and 2015 SSEs (Chen et al., 2018). The computations include 10 × 10 km spatial grids parallel and perpendicular to the Ryukyu Trench while projecting them on the surface (Figure 1). We determine the cumulative slip and fault dimensions of each SSE from 200 bootstrap samples in a 0.1 magnitude bin for their mean values and the resulting Coulomb stress changes (Figure 2A). Note that the fault dimensions are not assigned and are scaled only by the SSE moments (Figure 2B). We perform a checkerboard test on the interface fault to confirm that the resolution of offshore cumulative slip of SSEs is valid. Here the input slip was assumed as 0.1 m homogeneously over the checkerboards of the entire subfault plane (Supplementary Figure S3A) to examine the reproduced slip pattern. We apply the resolution matrix (e.g., Yabuki and Matsu’ura 1992; Yokota et al., 2016) as below: R = (H ^T H +α² G ^T G)⁻¹ H ^T H, where H is the static-response-function matrix; α is the hyperparameter of smoothness determined by a Bayesian information criterion (ABIC; Akaike, 1980); G is the spatial smoothness matrix. T denotes the transposed matrix. The resolution value as diagonal elements of the resolution matrix shows a poor resolution (< 0.15) of the offshore slip when the slip area is 50 km away from the coastline (Figure 3). Thus, we determine the cumulative slip and fault dimensions of each SSE within a 50 km distance.

The SSE moment greater than M_w 7.0 will produce an offshore slip 50 km away from the coastline. This would cause significant uncertainties in the following analysis, and thus we only determine the cumulative slip and fault dimensions smaller than the SSE moment of M_w 7.0. The Coulomb stress change Δα_c is Δ̰σ_c = Δτ - μΔσ _n , where Δτ is the shear stress change, Δσ_n is the normal stress change, and μ is the effective friction coefficient. Δτ and Δσ_n are estimated from the 3D strain field on the specified receiver faults derived from the SSE cumulative slip on the subducting plate interface (the source fault). The SSE cumulative slip is multiplied by the elastic stiffness to obtain the 3D strain field and Δσ_c. Poisson’s ratio (PR) and Young’s modulus (E) are 0.25 and 80 GPa, respectively. The friction coefficients are 0.4 and 0.2 for the region of the 2005 SSE and the 2009/2015 SSEs, respectively, and are related to their different source depths along the oceanic subducting plate interface (Kaneki and Hirono, 2019). Thus, the shear modulus is G = E/[2(1 + PR)] = 30 GPa.

We designed a receiver fault where each subfault plane had a specified strike, dip, and rake, and combine them into the source fault model to calculate the SSE static stress transfer. The receiver faults begin from all patches outside of the primary SSE cumulative slip region, and thus the Coulomb stress change is positive on the receiver faults. In addition, the Coulomb failure is estimated from the SSE cumulative slip theoretically for the stress drop: Δσ ≈ G D / A , where D is the average slip on the patch and A is the patch area. Since the 2005 SSE occurred at the downdip end of the Ryukyu fault, the stress changes on the Ryukyu fault caused by the synthetic SSEs may be positive or negative simultaneously on the spatial scales of the 10 × 10 km grids. A megathrust event may require homogeneous stress accumulation that is strong enough over most of the locked region. Therefore, we follow Kaneko et al. (2018) to calculate the mean stress changes in the entire fault region (Figure 1) caused by SSEs. In all bootstrap samples, we notice that the shear stress changes seem more significant than the normal stress changes caused by the observed and synthetic SSEs (Figure 4). The shear stress changes averaged over the entire fault region are distributed approximately from 0.005 to 0.03 MPa, whereas the normal stress changes are mostly lower from 0.003 to 0.01 MPa (Figure 4). A recent study regarding the effective normal stresses during the Boso SSEs in Japan suggests that the normal stress changes on the subducting plate interface were much smaller than the lithostatic pressures (Kobayashi and Sato, 2021). For simplicity, we used the shear stress component of each SSE for the resulting probability estimates and only showed the shear stress changes on and around the fault areas of SSEs in the results.

To evaluate the probability of a megathrust earthquake on the Ryukyu fault in a given time interval, we need to know the timing of the stressed conditions in the earthquake cycle. Establishing earthquake cycles requires the timing and size of large historical earthquakes. This modeling approach addresses the evolution of shear stress as earthquake cycles of many megathrust events averaged over the locked region. The total seismic moment released by all megathrust events equals the geodetically estimated moment deficit (Kaneko et al., 2018). Thus, stress changes caused by SSEs can be regarded as the amplitudes of stress perturbation relative to the stress accumulation and stress decreases in large earthquakes. This simple approach has been used for potential seismic hazard assessments in the Nankai, Himalayas, and Hikurangi subduction zones, where the records of historical earthquakes are incomplete (Parsons et al., 2012; Stevens and Avouac, 2016; Kaneko et al., 2018). The establishment of a synthetic shear stress earthquake cycle on the Ryukyu fault requires three parameters: the earthquake stress drop (Δσ), earthquake recurrence intervals (tr), and the background tectonic stress rate (τ). We used three assumptions for parameters to ensure that the synthetic earthquake cycles reproduced the natural characteristics. First, Δσ during the earthquake occurrences is a log-normal distribution with a standard deviation log₁₀Δσ of 0.4 and a mean Δσ of 2.0 MPa, which is a general estimate of Δσ from global subduction earthquakes (e.g., Allmann and Shearer, 2009; Kaneko and Shearer, 2015; Courboulex et al., 2016). Δσ is randomly distributed over time and increases with the earthquake size, which means that the rupture dimension of each earthquake was not assigned and was scaled only by the earthquake size and stress decrease.

Second, tr increases with earthquake size following the GR Law: log10N = a - bM (Gutenberg and Richter, 1944). The earthquake magnitude-frequency relation means a megathrust earthquake has a longer tr and recurs less frequently than a large earthquake in many earthquake cycles. Kaneko et al. (2018) assumed the b value to be empirically 1.0 to establish t _r for M_w ≥ 7.8 earthquakes in the southern Hikurangi-locked region. It has been shown that b values can change in different subduction zones from 0.7 to 1.4 and possibly vary with the stress state (Nishikawa and Ide, 2014; Scholz, 2015; Petruccelli et al., 2019). b values in the subduction locked region have been estimated to be lower than 1.0 in the Nankai Trough due to a higher stress state (Nanjo and Yoshida, 2018). However, the subduction zone offshore typically lacks constraints by large historical earthquakes for constraining b values in their megathrust earthquake cycles. We construct tr for geodetically estimated M_w ≥ 7.5 earthquakes on the locked Ryukyu fault (Hsu et al., 2012) using the past earthquakes in this subduction zone. There are sixteen M_w ≥ 7.0 earthquakes that have occurred in the last 120 years and have been relocated by two local earthquake catalogs: 1) Cheng and Yeh (1989) and Theunissen et al. (2010) and 2) Chang et al. (2016). The historical earthquake data provide preliminary constraints on the characteristics of the b value for the megathrust earthquake cycle. Figure 5 shows the hypocenters of the past M_w ≥ 7.0 earthquakes relocated from Cheng and Yeh (1989) and Theunissen et al. (2010). Most hypocenters are located along the subducting plate interface with a depth difference from the plate interface of fewer than 15 km (Figure 5). Since the difference is within the uncertainty of relocated hypocenters (Cheng and Yeh, 1989; Theunissen et al., 2010), we assume that the earthquakes are almost interplate events that can be used to estimate the b value along the subducting plate interface.

The b value is estimated by the maximum-likelihood method (Aki, 1965), and we quantify the uncertainty following Shi and Bolt (1982). The b value of the southernmost Ryukyu subduction zone is approximately 1.20–1.25, constrained by the historical earthquakes from M_w 7.0 to 7.4 in the two catalogs (Figure 5). The regression line of the b value is no longer constrained if we involve M_w ≥ 7.5 earthquakes in the b value estimation in both catalogs (Figure 5). We note that the two M_w 8.0 earthquakes in Chang et al. (2016) do not seem to follow the same magnitude-frequency relation as M_w 7.0 to 7.4 (Figure 5). Since no earthquakes are estimated within the range between M_w 7.4 to 8.0, we presume the magnitudes of the two M_w 8.0 earthquakes from Chang et al. (2016) may be overestimated. The regression line of the b value determined from the Cheng and Yeh, (1989) and Theunissen et al. (2010) catalog can explain the data in M_w 7.7 (Figure 5). The estimated b value of 1.2 is higher than the average of 1.0 in the entire Ryukyu subduction zone estimated from magnitudes of M_w 5.0 to 7.0 (e.g., Petruccelli et al., 2019).

Some studies have pointed out a change in the b values for larger earthquake magnitudes (e.g., Utsu, 1999; Wiemer and Wyss, 2000). Chang et al. (2016) have shown that b values are approximately 1.0 in the earthquake magnitudes from M_w 5.0 to 7.0 in their relocated historical catalog. The b values increase to over 1.2 in the earthquake magnitudes greater than M_w 7.0. It is consistent with our b value estimation for the past M_w ≥ 7.0 earthquakes that indicate a possible difference in earthquake magnitude-frequency relation across the magnitudes of M_w 7.0 in the southernmost Ryukyu subduction zone. The increase in the b value for larger earthquake magnitudes may result from either the short history of instrumental records or the nucleation process of large earthquakes. A complete historical earthquake catalog in the future will resolve whether the b values increase with larger earthquake magnitudes. At this stage, we use the current observations of the 1.2 b value to establish possible megathrust earthquake cycles on the Ryukyu fault. We assume the locked Ryukyu fault region shares the same earthquake magnitude-frequency relation as we observe in the southernmost Ryukyu subduction zone to establish t _r for the M_w ≥ 7.5 earthquakes. M_w 7.5 is the minimum earthquake size of the Ryukyu fault in our synthetic earthquake cycles. This result is consistent with the geodetically estimated moment deficit for a minimum of M_w 7.5 in this locked region (Hsu et al., 2012). To evaluate the maximum M_w 8.7 event as suggested (Hsu et al., 2012), we assume that earthquakes with magnitudes from M_w 7.5 to 8.7 follow the same estimated b values (dotted line in Figure 5). Considering the uncertainties of the b value are approximately ± 0.05 to 0.06, the b value could range from 1.15 to 1.31 in both catalogs (Figure 5). However, if the two M_w 8.0 earthquakes from Chang et al. (2016) are not overestimated, the b value will be close to 1.0 with a 0.2 difference compared to our estimates. Thus, we tested the b values between 1.0 and 1.4 with an 0.2 interval to establish the t _r on the Ryukyu fault.

After determining the b values, the t _r and τ can thus be designed in the earthquake cycles. Here t _r is the same as Δσ and is randomly distributed in the earthquake cycles, resulting in unpredictability using time and slip. In this sense, τ on the Ryukyu fault is the ratio of Δσ to t _r , and τ ˙ = ∑ i = 1 N Δ σ i ∕ ∑ i = 1 N t r i , where N is the total number of simulated earthquakes, assumed to be 50,000 (Kaneko et al., 2018). Their seismic moments are balanced by the geodetically estimated moment deficit (Ṁ _o ). Ṁ _o = µV _pl A, where Ṁ _o is the moment deficit rate, µ is the shear modulus, and V _pl is the geodetic slip deficit rate of the Ryukyu fault estimated from the fault dimension (A; the dotted rectangle in Figure 1). We use Vpl of 8.6 cm/yr (Hsu et al., 2012) to estimate the seismogenic potential. However, the Vpl may be overestimated by the onshore GNSS observations. It could also be accommodated by aseismic processes, such as the 2005 SSE, in the northwest-downdip region of the Ryukyu fault (Figure 1). Seafloor GNSS-acoustic observations in the Nankai Trough also supported the overestimation of 10–30 percent of the V _pl value, determined from only the onshore GNSS data (Yokota et al., 2016). Therefore, we also use a lower, suitable V _pl value of 6.0 cm/yr for the Ryukyu fault to estimate the seismogenic potential. Under the framework of a 1.2 b value, we use the constraints of 50 years recurrence intervals from M_w 7.0 historical events (Cheng and Yeh, 1989; Theunissen et al., 2010) and the geodetically inferred recurrence interval of 500 years for an M_w 8.7 event (Hsu et al., 2012) to scale the corresponding recurrence intervals from M_w 7.5 to 8.7 events. We assume the observed recurrence interval of 150–400 years in M ∼8.0 events at the adjacent southern Ryukyu Islands (e.g., Araoka et al., 2013; Ando et al., 2018) is compatible with our study area to help scale the t _r. Thus, the average t _r is approximately 180 years for M_w ≥ 7.5 earthquakes and 280 years for M_w ≥ 8.0 earthquakes. The tested frameworks of 1.0 and 1.4 b values are scaled similarly as in the 1.2 b value in the earthquake cycles but with different magnitude frequencies.

As the timing of the downdip SSEs is unknown in the Ryukyu fault earthquake cycle, we calculate the stress change caused by those SSEs in the synthetic earthquake cycles every t _pert years, where t _pert is the recurrence interval of the SSEs (Kaneko et al., 2018). It has been shown that M_w 6.0 SSEs recurred every 2 years, and M_w 7.0 SSEs recurred every 6 years in the Nankai, Hikurangi, and Guerrero subduction zones (Wallace et al., 2012; Obara and Kato, 2016; Radiguet et al., 2016). Since the recurrence interval of SSEs remains unclear in the study area (Chen S. K. et al., 2022), we assume that t _pert is proportional to the moment of SSEs, as observed in those subduction zones. Note that the resulting probability of a large subduction earthquake is independent of the choice of t _pert when the number of simulated megathrust events is designed to be large enough (Kaneko et al., 2018). As the observed duration of M_w ∼6.4 to 6.6 SSEs is two to 4 months (Chen et al., 2018), the M_w 6.0 to 7.0 SSEs are assumed to last 1 week to 6 months to agree with that observed in global subduction zones (e.g., Wallace et al., 2012; Obara and Kato, 2016; Cruz-Atienza et al., 2021). Taking t _pert and SSE duration into consideration for the time resolution of probability estimation, we calculate the probability of an M_w ≥ 7.5 earthquake with a maximum moment of M_w 8.7 on the Ryukyu fault every year.

For the timing of the earthquake cycle of the Ryukyu fault, there was no paleoseismological evidence for a magnitude greater than our assumed minimum size of M_w 7.5. Only the 1920 M_w 7.7 earthquake occurred in proximity to the locked region (e.g., Theunissen et al., 2010). Since the hypocenter of historical earthquakes is poorly resolved, the 1920 earthquake cannot be ruled out as an interface slip event. Here, we assume that the 1920 earthquake was the last rupture on the Ryukyu fault. The timing of the first SSE in 2005 is at least 85 years into the earthquake cycle to obtain a specific resulting probability. For computing efficiency, we exclude the next 85 years following each Δσ from the probability calculations. We note the possibility that the 1920 earthquake was not the last rupture of the Ryukyu fault and the interface remained locked since 1,604. Thus, we also model the timing of the 2005 SSE that was 400 years into the earthquake cycle. Since the three SSEs have shown a time-dependent t _pert, we evaluated whether the stress level of each t _pert exceeds the triggering stress threshold of the next earthquake in a time interval of one to 2 years immediately after each t _pert. We apply three possible scenarios that consider each t _pert as positive stress in the timing advance. The closer the timing of the stress perturbation caused by the SSE is to the triggering stress threshold of the next event, the higher the seismogenic probability is. For the details and calculation of the resulting probability and uncertainties, refer to Section 2.3 in Kaneko et al. (2018).

Figure 6 shows shear stress changes of M_w 6.0 to 7.0 SSEs at the two source depths on the subducting plate interface. The regions of the shear stress decrease, i.e., stress drop, in each subfigure approximately correspond to the SSE cumulative slip regions. The size of the stress drop area is proportional to the moment of SSE. The shear stress increase may have occurred around the stress drop region, and the maximum stress increase appears at the edge of the stress drop region. For SSEs that occur downdip of the Ryukyu fault, a shear stress increase of 0.05 ± 0.02 MPa by an M_w 6.0 SSE is estimated within the downdip region (Figure 6A). There is no shear stress increase on the updip side since the slip dimension of the M_w 6.0 SSE is insufficient for a broader influence on the shear stress increase. When the moment of the SSE reaches M_w 6.5, the shear stress increase of 0.2 ± 0.03 MPa is visible within much of the fault region at northwest-downdip depths (Figure 6B). Note that shear stress decreases begin to appear at the downdip depth limit of the Ryukyu fault. This result is not surprising because the slip dimension of the SSE is proportional to the moment of the SSE, which changes the spatial distribution of the shear stress increase and decrease. The shear stress increase in an M_w 7.0 SSE reaches 0.3 ± 0.05 MPa within much of the fault region at northwest-downdip depths (Figure 6C). The above shear stress increases are less than one-tenth of the estimated values when averaged over the entire fault region (Figure 4). The results are thus not greater than 0.005–0.03 MPa, which is broadly consistent with the stress changes on locked faults caused by similar SSE sizes in large subduction zones (e.g., Nakamura and Sunagawa, 2015; Radiguet et al., 2016; Voss et al., 2018; Saltogianni et al., 2021). The shear stress increase is insignificant in each scenario where SSEs occur within much of the Ryukyu fault downdip region. We observed only a shear stress increase of 0.05 ± 0.04 MPa, estimated from an M_w 6.0 SSE (Figure 6D). Such slight stress increases of 0.10 ± 0.05 MPa and 0.15 ± 0.06 MPa are determined from M_w 6.5 and 7.0 SSEs, respectively (Figures 6E,F). The shear stress increase in the entire locked region does not seem proportional to the moment of the SSE.

Figure 7 shows the annual probability as a function of time since the first 2005 SSE based on each shear stress increase in Figure 6. Here, we demonstrate the probabilities of an M_w ≥ 7.5 earthquake with a maximum moment of M_w 8.7 and a magnitude frequency relation of a b value of 1.2 (blue curves, Figures 7A–C). The annual probability in the first year following the SSE in 2005 only increases slightly. The probability change is within the uncertainty of background probability in the 10 year sequence. Note that the resulting probability uncertainty is always the largest since the influence of both shear stress decrease and increase on the Ryukyu fault. We observed that the 2009 and 2015 SSEs show a significant probability increase over the uncertainty of background probability (Figures 7A–C). In each scenario, the annual probabilities in the first year following the SSEs are approximately 1.1–1.2 times the background probability. As expected, the ratio of the probability increase is governed by the amount of shear stress increase over the first year (Kaneko et al., 2018). Note that the probability increase looks small because the stress perturbations from the SSE in Figure 6 are averaged over the entire fault region in Figure 4. The annual probability reduces to the background level of 0.45% after 1 year in each case because we assumed the durations of all SSEs are less than half of a year. Thus, the probability does not change when the external stress perturbation is absent, relying on the failure criterion that a large earthquake occurs immediately at a stress threshold (e.g., Mazzotti and Adams, 2004).

We also consider the annual probability for b values of 1.0 and 1.4 for the earthquake cycles. In both scenarios, the potential of a large earthquake is the highest over the first year of each SSE, namely the result of b values of 1.2. We find that the background probability is affected by the b value used for t _r in synthetic shear stress earthquake cycles, which was found empirically in Kaneko et al. (2018). In the case of the b value of 1.0 (Figures 7D–F), the background probability increases to 0.60% and then slightly increases with the moments of SSEs. When the b value increases to 1.4, we observed a decrease in background probability to 0.30% and then slightly increases with the same ratios as Figures 7A–C. The background probability is dependent on the choice of the b value because τ is determined from ∑ i = 1 N Δ σ i ∕ ∑ i = 1 N t r i . Additionally, the background probability is independent of the assumed Δσ because τ also decreases linearly with Δσ (Kaneko et al., 2018). If so, it implies that a higher b value is used in the t _r, and a lower frequency of megathrust events is a release of the total Δσ.

We find that the background probability of an M_w ≥ 7.5 subduction earthquake is small on the Ryukyu fault, between 0.30% and 0.60%, compared to the predictions in the Nankai and southern Hikurangi subduction zones (Kaneko et al., 2018). Kaenko et al. (2018) assumed that the total seismic moments in the two subduction zones were released by assumed M_w ≥ 7.8 earthquakes with an empirical b value of 1.0 in their earthquake cycles. The total Δσ is thus balanced by a higher frequency of megathrust events, which leads to a higher background probability than our results. Our results are more robust based on the geodetically estimated minimum earthquakes of M_w ≥ 7.5 on the Ryukyu fault (Hsu et al., 2012) and the historical earthquake b values. It has been shown that the background probability of a large subduction megathrust is not affected much by the V _pl value (Kaneko et al., 2018). Although we assumed a higher V _pl value than that assumed in Kaneko et al. (2018), our results show that in the case of the V _pl = 8.6 cm/yr (green curves, Figure 7): the background probability will be higher than V _pl = 6.0 cm/yr of only about 0.3%. We find the effect of V _pl is secondary to the probability increase compared to that of the timing of SSE into the earthquake cycle (red curves, Figure 7). It could be the case if the V _pl is very large enough to produce a similar probability increase in the case of 400 years into the earthquake cycle. However, the present observations in this subduction zone do not support this kind of scenario. Thus, V _pl is unlikely to be a critical parameter that results in the difference in background probability.

The seismogenic potential of an M_w ≥ 7.5 earthquake on the Ryukyu fault is estimated to be very small when considering the fault zone can rupture entirely as an M_w 8.7 event. The low probability relies on the 1920 M_w 7.7 earthquake, which was assumed to be the last rupture of the Ryukyu fault. Thus, the timing of the first SSE in 2005 is 85 years in the early stage of the earthquake cycle relative to a 500 years recurrence of M_w 8.7 events on the Ryukyu fault. There is a possibility that the 1920 M_w 7.7 earthquake did not occur on the fault interface, which may imply that the Ryukyu fault has remained locked for at least the last 417 years because there are no records of a geodetically estimated minimum magnitude of M_w 7.5 during that period. The stress level of the Ryukyu fault is higher than this assumption, and the 2005 SSE may have changed into the late stage of the earthquake cycle, much closer to the end of the 500 years recurrence of M_w 8.7 events. In this scenario, we obtained a higher annual probability of the megathrust earthquake with a maximum of M_w 8.7, perhaps close to 1.2%, with the possibility of > 1.2% over the next decade (Figure 7). This scenario highlights the potential of a coming megathrust earthquake of M_w 7.5 to 8.7 in the next 100 years. Thus, we discuss three possible rupture scenarios of the Ryukyu fault (Figure 8): (A) regular and repeated ruptures with M_w 7.0 to 7.5 events on several segments of the entire locked region. (B) The Ryukyu fault ruptures separately in two to three main sections, with M_w 7.5 to 8.2 for each earthquake size in the locked area. (C) The entire Ryukyu fault ruptures coherently as an M_w 8.7 event (Hsu et al., 2012), even if we lack paleoseismological evidence to confirm the existence of such a large M_w 8.7 event in the past. Types (A) and (C) are the end-member scenarios, and type (B) is a mixture between them.

If we consider only the hypocenters of the M_w ≥ 7.0 historical earthquakes in the last 120 years, the M_w 7.0 to 7.2 events occurred almost entirely in the western half of the Ryukyu fault (Figure 5). It may imply that the fault interface is locked on the western-half plane but not fully locked over the entire fault interface. In contrast, the fault interface of the eastern-half plane may be weakly locked since there were no frequent M_w ≥ 7.0 earthquakes but was replaced by very low frequency earthquakes (VLFEs) (Nakamura and Sunagawa, 2015). However, considering the uncertainties in the locations, those events in Figure 5 may not all have occurred on the slip interface linked to the Ryukyu fault. If they were, it could imply that the frictional properties on the fault interface exhibited a time-dependent difference in that period between the western and eastern half. Time-dependent frictional variations have been inferred from geodetically observed SSEs along the strike direction of the plate interface at depths of almost 25–45 km (Chen et al., 2018).

The evidence of VLFEs activity on the eastern-half plane only (Nakamura and Sunagawa, 2015) also supports the difference in the interface frictional behaviors. Recent seafloor geodetic evidence indicates that the convergence rate is likely to increase along the strike of the Ryukyu fault offshore from 92 to 123 mm/yr (Chen H.Y. et al., 2022). The new observations may reveal a strain partitioning between the western/eastern-half planes, complementary to the locations of M_w ≥ 7.0 historical earthquakes and VLFEs. If lateral frictional variation exists in the entire locked region, types (A) and (B) are more similar to the possible rupture behaviors of the Ryukyu fault as fault segmentation. The idea of fault segmentation is reasonable for the higher observed b value of 1.2 for the M_w 7.0 to 7.7 historical earthquakes, rather than the empirical 1.0 value. Seafloor geodetic measurements across the entire fault region will help us understand the frictional interface properties regarding possible fault segmentation. We propose that future M_w 7.0 to 7.5 earthquakes could occur frequently offshore eastern Taiwan, the same as the past events very close to the coastline (Figure 5). The close epicenter distances to the Taiwan region would lead to local destruction and a potential tsunami threat.

A large earthquake triggered by SSE activity is still an open question and has been extensively investigated in global subduction zones. The first observed cases regarding the triggering or delayed triggering of the large earthquakes were rare. The 2014 M_w 7.3 Papanoa earthquake was estimated to be triggered by a 0.04 MPa ongoing stress perturbation from the downdip SSE 2 months before the M_w 7.3 event (Radiguet et al., 2016). The 2018 M_w 6.9 Zakynthos earthquake was likely triggered by the downdip SSE with a 0.025 MPa stress perturbation, which terminated 2 months before the M_w 6.9 event (Saltogianni et al., 2021). In the Costa Rica subduction zone, the 2012 SSE that occurred downdip the 2012 M_w 7.6 Nicoya earthquake was delayed by at least 6 months (Voss et al., 2018). The stress perturbation of 0.01 MPa caused by the SSE on the entire Nicoya locked region may have been insufficient for the triggering stress threshold of the Nicoya earthquake, as could the delay (Voss et al., 2018). The common point is that they were very close to the end of their earthquake cycles. It means the stress levels of the locked faults were high enough, and thus, fault ruptures could have been triggered by small stress perturbations caused by SSEs. In the southernmost Ryukyu subduction zone, our estimated shear stress increase on the Ryukyu fault caused by the 2005, 2009, and 2015 SSEs is similar to those cases when averaged over the entire locked region (less than 0.005–0.04 MPa).

A megathrust earthquake has not yet been triggered on the Ryukyu fault after a sequence of SSEs; it is therefore very likely that the end of the earthquake cycle is not close as in those cases. However, this is not true if the SSEs are very close to the end of the earthquake cycle and there has been no large earthquake at the fault interface in the past 450 years, as in the case of the 2011 M_w 9.0 Tohoku earthquake (Kato et al., 2012; Ito et al., 2013). Our results suggest that the 2009 and 2015 SSEs can produce a significant probability increase in the seismogenic potential of the Ryukyu fault. Since the two SSEs occurred on the same patch, a future SSE recurring there with a similar moment can be expected. Note that our results only show the influence of M_w 6.0 to 7.0 SSEs on the triggering of the megathrust earthquake. If a future M_w > 7.0 SSE occurs down there and the offshore slip cannot be effectively detected (Figure 3), such SSE could produce more stress perturbations than our scenarios. We alert the Taiwan region to monitor the local SSE activity and notice the potential of a large earthquake triggered by the SSEs downdip the locked Ryukyu fault. In addition, the 2005 SSE located at the downdip depth limit of the locked Ryukyu fault seems to increase the potential of a large earthquake slightly. Time-dependent frictional variations at the fault interface could be a reason for SSE generation, as observed in the Japan Trench and Costa Rica subduction zones (e.g., Ito et al., 2013; Dixon et al., 2014). If so, it remains a subject of future work hampered by the limited observations in the southernmost Ryukyu subduction zone for the time-dependent frictional variations.

We report the likelihood of a megathrust earthquake on the Ryukyu fault in the southernmost Ryukyu subduction zone for the first time. The annual probability of an M_w ≥ 7.5 earthquake with a maximum earthquake size of M_w 8.7 is approximately 0.3 to 0.6%, constrained by the observed local b value. The b value for historical earthquakes with magnitudes from M_w 7.0 to 7.7 is 1.2, higher than the average of 1.0 b value derived from M_w 5.0 to 7.0 earthquakes in the entire Ryukyu subduction zone. It implies that the recurrence of megathrust earthquakes is longer than previously thought on the Ryukyu fault if the megathrust events follow the observed magnitude-frequency relation. We also estimate the stress changes on the Ryukyu fault caused by the SSEs in 2005, 2009, and 2015. The shear stress increases averaged over the entire locked zone for a megathrust event range from 0.005 to 0.04 MPa, broadly consistent with the previous observations. Such stress increases have caused the probabilities of a megathrust earthquake triggering over the 2009 and 2015 SSEs 1.1 to 1.2 times relative to the background probability. A future SSE recurring downdip the locked Ryukyu fault should be monitored for the triggering of the large earthquake. It is very likely the earthquake cycle is not close to the end; thus far, there has been no triggered megathrust earthquake after a sequence of SSEs. This work provides unique observation for the critical issue of slow slip-triggered large earthquakes.