Time-independent forecast model for large crustal earthquakes in southwest Japan using GNSS data

In this study, we developed a regional likelihood model for crustal earthquakes using geodetic strain-rate data from southwest Japan. First, the smoothed strain-rate distributions were estimated from continuous Global Navigation Satellite System (GNSS) measurements. Second, we removed the elastic strain rate attributed to interplate coupling on the subducting plate boundary, including the observed strain rate, under the assumption that it is not attributed to permanent loading on crustal faults. We then converted the geodetic strain rates to seismic moment rates and calculated the 30-year probability for M ≥ 6 earthquakes in 0.2 × 0.2° cells, using a truncated Gutenberg–Richter law and time-independent Poisson process. Likelihood models developed using different conversion equations, seismogenic thicknesses, and rigidities were validated using the epicenters and moment distribution of historical earthquakes. The average seismic moment rate of crustal earthquakes recorded during 1586–2020 was only 13–20% of the seismic moment rate converted from the geodetic data, which suggests that the observed geodetic strain rate includes considerable inelastic strain. Therefore, we introduced an empirical coefficient to calibrate the moment rate converted from geodetic data with the moment rate of the earthquakes. Several statistical scores and the Molchan diagram showed all models could predict real earthquakes better than the reference model, in which earthquakes occur uniformly in space. Models using principal horizontal strain rates exhibited better predictive skill than those using the maximum horizontal shear strain rate. There were no significant differences in predictive skill between uniform and variable distributions for seismogenic thickness and rigidity. The preferred models suggested high 30-year probability in the Niigata–Kobe Tectonic Zone and central Kyushu, exceeding 1% in more than half of the analyzed region. The model predictive skill was also verified by a prospective test using earthquakes recorded during 2010–2020. This study suggests that the proposed forecast model based on geodetic data can improve the regional likelihood model for crustal earthquakes in Japan in combination with other forecast models based on active faults and seismicity.


Introduction
The development of a Regional Earthquake Likelihood Model (RELM) is a practical product of earthquake science that can be used to assess earthquake hazards and enhance human preparedness for earthquake disasters. Most RELMs use either active faults or past seismicity to estimate long-term earthquake rates (e.g., Zechar et al. 2013). Crustal deformation measured by geodetic instruments including the Global Navigation Satellite System (GNSS) can also provide useful information on future earthquakes because earthquakes convert elastic deformation to inelastic permanent deformation according to Reid's (1911) elastic rebound theory (Rhoades et al. 2017). For example, Ward (1994) stated that geodetic data can be used to quantify the earthquake potential of poorly known faults during seismic hazard analysis in southern California. Forecast models using GNSS data have successfully forecast large earthquakes on both global (e.g.,  and regional scales (e.g., Shen et al. 2007).
The Japanese Islands are situated in a region of complex tectonics, where four tectonic plates converge (e.g., Taira 2001). Regional tectonic stress leads to numerous earthquakes, which can be categorized into two types: subduction-type earthquakes and crustal earthquakes. The Headquarters for Earthquake Research Promotion (HERP), which is an administrative branch of the Japanese government, evaluates the long-term activity of both earthquake types, and develops national seismic hazard maps of Japan . Their evaluation of crustal earthquakes is typically based on geological and geomorphological data from major active faults with the potential for M ≥ 7 earthquakes. However, several damaging crustal earthquakes, including the 2008 M7.3 Iwate-Miyagi Inland earthquake (Ohta et al.2008), have occurred on faults that are not recognized as major active faults. The HERP has recently begun conducting regional evaluations of active faults to assess M ≥ 6.8 earthquakes using short active fault and seismicity data. The HERP also considers earthquakes on unspecified faults based on seismicity data for national seismic hazard maps. However, the HERP does not currently use geodetic data, although contemporary deformation of the Japanese Islands is monitored by a permanent dense GNSS network (e.g., Sagiya 2004). Only a few pioneering studies have developed forecast models for the Japanese Islands using geodetic data (Sagiya 2015;Takahashi and Shinohara 2015;Triyoso and Shimazaki 2012). Therefore, in this study, we utilized GNSS data to develop a RELM for crustal earthquakes in southwest Japan and compare it with the distribution of actual earthquakes to verify its predictive ability.

GNSS data
We used velocity data from permanent GNSS stations in southwest Japan, as estimated by Nishimura (2017). The secular velocity was estimated using the daily coordinates from April 2005 to December 2009. Nishimura (2017) removed the elastic deformation attributed to interplate coupling between the continental plate and the subducting Philippine Sea plate along the Nankai Trough from the observed velocity from the GNSS stations. They noted that the distribution of crustal earthquakes is correlated with that of the strain rates after removal. Here, we assume that the stress causing crustal earthquakes in southwest Japan is not directly loaded by coupling of the subduction megathrust but by permanent inelastic deformation below seismogenic faults (e.g., Iio and Kobayashi 2002). In other words, elastic deformation attributed to interplate coupling is fully released by interplate megathrust earthquakes; therefore, it can modulate the timing of crustal earthquakes but does not permanently load crustal faults (e.g., Hori and Oike 1999). This assumption is valid for a linear viscoelastic medium as well as an elastic medium (Mitogawa and Nishimura 2020). Therefore, we used the corrected velocity after removal of elastic deformation and estimated the strain rate distribution using the method proposed by Okazaki et al. (2021), in which the velocity field is objectively optimized by the expansion of the bicubic B-spline basis function. Figure 1 shows the estimated strain-rate distribution at 0.2 × 0.2° cells. The plotted strain component is the maximum of |ε h1 | , |ε h2 | , and |ε h1 +ε h2 | , where ε h1 and ε h2 are the maximum and minimum principal strain rates in the horizontal plane, respectively, which are used in the conversion to seismic moment rates in the next subsection. The other strain components are shown in Additional file 1: Fig. S1. These generally indicate that many M ≥ 6 crustal earthquakes occur in cells with high GNSS strain rates. We regard this strain-rate distribution as a proxy for the long-term average value of strain accumulation causing crustal earthquakes.

Method for calculating earthquake probability
First, we converted the observed strain rates to seismic moment rates, ṁ 0 . Several equations for this conversion have been proposed, of which the following were tested (Savage and Simpson, 1997;Ward, 1994): where µ , H , and A are the rigidity, thickness of the seismogenic layer, and surface area, respectively. ε h1 , ε h2 , and (1) γ are the maximum and minimum principal strain rates, and the maximum engineering shear strain rate in the horizontal plane, respectively. Here, we add the empirical correction coefficient β , which is like the seismic coupling coefficient introduced by Bird and Kagan (2004), and to the original equations proposed by Ward (1994) and Savage and Simpson (1997). This coefficient calibrates the moment rates converted from the geodetic strain rates to those of actual crustal earthquakes. The most important reason for introducing the empirical coefficient is that the geodetic strain rates were much higher than the seismic and geological strain rates in Japan, as reported in previous studies (e.g., Shen-Tu et al.1995). A part of this discrepancy can be explained by elastic deformation attributed to interplate coupling on the subduction megathrust fault (e.g., Ikeda 2014). The discrepancy remains after we removed the elastic strain attributed to interplate coupling, which suggests that the observed strain includes considerable inelastic strain. The inelastic strain is attributed to ongoing permanent aseismic processes, including fault creep, slow slip events, active folding, and numerous small earthquakes occurring in the seismogenic layer. We estimate β by comparing the moment rates converted from geodetic strain rates (hereafter, geodetic moment rates) with those of past earthquakes in the following subsection. Bird and Liu (2007) proposed another conversion equation where ε 1 , ε 2 , and ε 3 are the principal strain rates in three dimensions. Because GNSS data can provide only horizontal strain rates, the three principal strain rates are derived by assuming the conservation of volume ( ε 1 +ε 2 +ε 3 = 0 ) and a vertical orientation of one principal strain rate. Using these assumptions, Eq. (4) is equivalent to Eq. (3), with β = 1.
The rigidity ( µ ) and seismogenic thickness ( H ) in Eqs. (1)-(3) generally exhibit regional variations. Therefore, we used D90, the depth above which 90% of earthquakes occurred, as a proxy for the seismogenic thickness (Fig. 2a). The D90 distribution of Omuralieva et al. (2012) was used after interpolation and extrapolation with adjustable tension continuous curvature splines (Wessel et al. 2013). We calculated the average rigidity of the seismogenic thickness based on the deep subsurface structure model of the Japan Seismic Hazard Information Station (J-SHIS) (National Research Institute for Earth Science and Disaster Resilience 2019) (Fig. 2b). Hereafter, we refer to the model that uses variable µ and H as the VARI model. The VARI models developed using (4) m 0 = 2µHAε 1 ,ε 2 < 0 −2µHAε 3 ,ε 2 ≥ 0 , Eq. (1)-(3) are termed VARI1, VARI2, and VARI3. The product of µ and H , which controls the geodetic moment rates, ranges between 51 GPa km and 979 GPa km, with an average of 351 GPa km. Moreover, we selected 30 GPa and 12 km as uniform µ and H values in the UNI model. The UNI models developed using Eq. (1)-(3) are termed UNI1, UNI2, and UNI3. The earthquake magnitude used in this study was the Japan Meteorological Agency (JMA) scale. The empirical relationship between the JMA magnitude ( M ) and seismic moment ( m 0 ) in Nm was used for shallow crustal earthquakes (Takemura, 1990): We then used a truncated Gutenberg-Richter law for the relationship between the magnitude ( M ) and number of earthquakes ( n): where M max is the maximum magnitude in the analyzed region. The number of earthquakes N (M) with magnitude M and larger is given by: The total seismic moment of earthquakes per unit time, M max −∞ m 0 (M)n(M)dM , following a truncated Gutenberg-Richter law, is assumed to be equal to the geodetic moment rate, ṁ 0 , converted from the strain rate (Eq. (1)-(3)). Then, a in Eq. (6) and (7) can be given by: Although the other parameters in Eq. (6) may be dependent on the region (e.g., Hirose and Maeda, 2011), we used 8.0 and 0.9 for M max and b as spatially uniform typical values for crustal earthquakes, respectively. M max was derived from the magnitude of the 1891 Nobi earthquake.
Assuming a stationary Poisson process, the probability of p(l, t) having l events in period t is given by: where Ṅ is the rate of earthquakes in the cell, which is a temporal derivative of the number of earthquakes. Then, the probability of one or more events can be expressed as: (5) logm 0 = 1.17M + 10.72.
We used Eqs. (7) and (10) to obtain the probability for a certain magnitude threshold and a certain time period.

Historical earthquake data
For calibration and validation of the RELM, we used a catalog of historical earthquakes before 1919 (Utsu 1990(Utsu , 2002. We also used an earthquake catalog determined using seismograms instrumentally recorded by the JMA in 1919, with a focus on large crustal earthquakes with depths of ≤ 20 km and M ≥ 6. The depth of historical earthquakes before 1919 is poorly understood. Therefore, we used earthquakes denoted as very shallow or shallow during 1885-1918 and all earthquakes before 1885 in Utsu's catalog. We manually excluded six M 6.0-7.4 historical earthquakes that the Earthquake Research Committee (2004) regarded as intraslab earthquakes similar to the 2001 M w 6.8 Geiyo earthquake (Kakehi 2004). A total of 151 earthquakes (indicated by red circles in Fig. 1 The number of earthquakes was counted using the number of epicenters in each cell. Therefore, only a single cell is responsible for an earthquake, although the source area of a large earthquake can occupy multiple cells. Because detailed source models for most historical earthquakes are not available, we simply distributed the seismic moment of an earthquake to surrounding cells to ensure that it was proportional to the area of a circle with a diameter of D in each cell. D in kilometers was calculated from the length of the fault in the empirical scaling law of Takemura (1998) as follows: The epicenters and moment distribution of historical earthquakes ( Fig. 3a) suggest that historical earthquakes are concentrated in areas of high strain rate, including the Niigata-Kobe Tectonic Zone (Sagiya et al. 2000), the San-in Shear Zone (Nishimura and Takada 2017), and central Kyushu , which have all been identified as strain concentration zones in previous studies. Temporal evolution of the total seismic moment in the analyzed area shows that the moment rate increased with time (Fig. 3b). The total moment was 8.68 × 10 20 Nm and 6.46 × 10 20 Nm during 599-2020 and 1586-2020, respectively. It is preferable to take a longer time to ensure a stable estimate of the average moment rate. However, the much higher rate in the later period suggests that many M ≥ 6 earthquakes were not documented before the sixteenth century. The moment rate, including two largest earthquakes during 1586-2020, was approximately equal to that during 1919-2020, i.e., the era of instrumental recorded seismograms, in which no M ≥ 6 earthquakes were likely missed (e.g., Ohta et al. 2002). Therefore, we used the period between 1586 and 2020 to calculate the average moment rate of long-term seismicity (1.49 × 10 18 Nm/yr). We determined the correction coefficient β to balance the average moment rate of historical earthquakes with the geodetic moment rates derived from Eqs. (1)-(3), as listed in Table 1.
The derived value of β ranges between 0.12 and 0.2, suggesting that most strain observed by GNSS measurements is not seismogenic one. It may be underestimated by the insufficient length of the earthquake catalog. The average moment rate can be almost double (2.73 × 10 18 Nm/year) if we adopt the period for the average between 1891 and 2020, just including the 1891 M8.0 Nobi earthquake which is the largest inland earthquake in history. It is also notable that the average moment rate for a certain period strongly depends on the uncertainties of the 0 5x10 20 1x10 21 500 1 000 1 500 2 000  599-1918599- (Utsu., 19902002) and instrumental seismograms from 1919 (Japan Meteorological Agency 2022) Table 1 Correction coefficient β used to balance the seismic moment calculated from geodetic strain rates to that of past large earthquakeṡ moment for the largest earthquakes, which were large before the era of instrumental recorded seismograms. However, the derived value of β is consistent with the previous studies that the observed geodetic strain rates are much larger than the seismological and geological ones (e.g., Shen-Tu et al. 1995;Sagiya et al. 2000).

Estimated probability distribution of crustal earthquakes
We calculated the probability of M ≥ 6 crustal earthquakes within 30 years for a total of six models using three conversion equations and two structures ( Fig. 4 and Additional file 1: Fig. S2). The highest probability in each model ranged from 2.9 to 3.7% and was found in Kyushu. The lowest probability was between 0.027 and 0.11%. More than half of the cells had a probability of ≥ 1% in all models. The calculated probability distribution for the UNI models (Fig. 4b) reflected the strain-rate distribution ( Fig. 1). Although the effect of variable thickness and rigidity on the probability was rather minor, we observed an increase in the probability due to a thick seismogenic layer in Hiroshima and Okayama Prefectures (132.5-134°E and 34-35°N) (Figs. 2a and 4a).

Retrospective testing of the forecast models
We calculated several statistical scores to examine and validate the performance of the forecast models (Table 2). All M ≥ 6 crustal earthquakes during 599-2020 (Figs. 1 and 3a) were used for the analysis. The first score, which was the correlation coefficient between the geodetic moment rate and the moment of past earthquakes, suggested weak positive correlations for all models. Scatter plots (Additional file 1: Fig. S3) showed no large moment release of past earthquakes in cells with low geodetic   moment rates, and the maximum moment release of past earthquakes increased with the geodetic moment rates. Although our forecast model is based on the proportionality between earthquake moment rates and geodetic moment rates on a long-term average, the period of the earthquakes employed in the model was apparently shorter than the average recurrence intervals of crustal earthquakes (more than hundreds to thousands of years), as manifested by the null moment of past earthquakes in many cells. The VARI models exhibited slightly higher correlation coefficients than the UNI models. Moreover, the models with maximum shear strain rates (i.e., VARI1 and UNI1) had lower correlation coefficients than those with principal strain components (i.e., VARI2, VARI3, UNI2, and UNI3). The second score was the probability gain per earthquake ( G ) for the models relative to a reference model in which earthquakes are uniformly distributed (e.g., Helmstetter et al. 2007): where L and L ref are the log-likelihood of the present model and reference model, respectively, and N is the number of earthquakes. The log-likelihood of the model is given by: where p i (l, t) is the probability of the i th cell, and K is the total number of cells.
The probability gains per earthquake exceeded 1.13 for the models with principal strain-rate components (i.e., VARI2, VARI3, UNI2, and UNI3), but were approximately 1.05 for models with the maximum shear strain rate (i.e., VARI1 and UNI1). This is attributed to the various focal mechanisms of crustal earthquakes in southwest Japan, including normal, strike-slip, and reverse faulting (Terakawa and Matsu'ura 2010). These gains may be less than our expectation, when compared with the values of 1.17 to 7.08 obtained in a seismicity-based forecast model in California (Helmstetter et al. 2007). However, all models were significantly better than the uniform density model because the probability gains of all 151 earthquakes exceeded 1,800, even for models that used the maximum shear strain rate.
The Molchan error diagram (Molchan 1997) is often used to demonstrate the performance of an alarm-based forecast. It represents the relationship between the miss rate, that is, the proportion of unforecasted earthquakes that does not satisfy a certain forecast threshold, and the proportion of forecasted space/time. In our case, the horizontal axis of Fig. 5 is the cell number in descending order of the geodetic moment rate, normalized by the total number of cells. We can arbitrarily choose a forecast threshold, such as 1% of the 30 year probability and the probability at the cells with the upper 25% probability, and alert the cells satisfying the threshold. The vertical axis represents the number and moment of earthquakes in the cells without the alert, normalized by the total number and moment of earthquakes, respectively. We also plotted the normalized geodetic moment rates in the corresponding cells. If there is no correlation between the model forecast and past earthquakes, the error curve is expected to be a diagonal line from the upper-left corner to the lower-right corner in the Molchan diagram. This indicates that the forecast model does not have a predictive skill. In contrast, a model with predictive skill has an error curve to the lower left of the diagonal line, indicating that earthquakes are concentrated in the cells with higher earthquake probability. We plotted the 95% confidence level of a non-correlation model simulated by a random cell number. The error curves for the number and moment of past earthquakes in all forecast models traced to the lower left of the 95% confidence level of the non-correlation model, which suggests a significant positive correlation between the geodetic moment rate and past earthquakes (Fig. 5).
Additionally, we used the proportion of forecasted earthquakes in the number and the moment at x = 0.25 in the Molchan diagram as a third score (e.g., Shen et al., 2007). The rest (e.g., the proportion of forecasted earthquakes) was obtained by subtracting the proportion of unforecasted earthquakes, as indicated by the error curves in Fig. 5, from one. Accordingly, 36-46% of the epicenters of past earthquakes were located in cells of the upper 25% of high probability/geodetic moment rates (Table 2). This proportion increased to 46 to 64% for the moment of past earthquakes. The UNI1 model had a smaller proportion of past earthquake concentrations than the other models. The fourth score was the area skill score (ASS), that is, the area above the error curve inside the unit square of the Molchan diagram. The ASS was used to indicate the overall performance of the covariate in earthquake forecasting models (Zechar and Jordan, 2010). The diagonal line of the Molchan diagram had ASS = 0.5 and the predictive model gave ASS > 0.5. The ASS for the epicenters of past earthquakes ranged from 0.62 to 0.65, whereas that for the moment ranged from 0.69 to 0.75. The ASS of models with the maximum shear strain rate was smaller than models with principal strain rate components for both UNI and VARI models.
These scores indicated past earthquakes were concentrated in a high-strain rate zone observed by GNSS and that the forecast model based on GNSS observations exhibited predictive ability for earthquakes. The models using the strain-rate component of Max(|ε h1 |, |ε h2 |) or Max(|ε h1 |, |ε h2 |, |ε h1 +ε h2 |) (i.e., UNI2, UNI3, VARI2, and VARI3) exhibited better performance than models using the maximum shear strain rate. There was no substantial difference between models with variable or uniform distributions of rigidity and seismogenic thickness. Considering that Eq. (3), which was proposed by Savage and Simpson (1997), uses Max(|ε h1 |, |ε h2 |, |ε h1 +ε h2 |) as an improved expression between the seismic moment and the surface strain rate after examining Eqs. (1) and (2)  used the UNI3 and VARI3 models as representative models for further analysis. Although all the statistical scores used in this section are related to the relative spatial distribution of geodetic strain rates and past earthquakes, the consistency of forecasted probability with the frequency of past earthquakes is also an important factor. Because we introduce β in Eqs. (1)-(3), the predicted seismic moment is equal to that of the past earthquakes. An annual rate of M ≥ 6.0 earthquakes predicted by our models is 0.23. It corresponds to ~ 100 and ~ 328 earthquakes for 435 and 1422 years, respectively. They are roughly equal to the actual number of earthquakes during 1586-2020 (i.e., 122) but much larger than that during 599-2020 (i.e., 155). It supports our inference that many M ≥ 6 earthquakes have not been documented before the sixteenth century.
Triyoso and Shimazaki (2012) constructed a forecast model using GNSS data and compared it with the previous 400-year-long seismicity. They showed that their GNSS-based model did not effectively reproduce the distribution of past crustal earthquakes and that its predictive skill was approximately equal to that of the uniform density model. Conversely, our study demonstrated that forecast models using GNSS data have better predictive skill. These contrasting findings can be attributed to notable differences between data analysis procedures. That is, we used a high-resolution distribution of the GNSS strain-rate derived from a larger number of GNSS stations, optimal smoothing of the strain rate distribution (Okazaki et al. 2021), and a finer coupling model on the subducting plate interface . Furthermore, the region analyzed was limited to southwest Japan in this study, but covers most of Japan in the work of Triyoso and Shimazaki (2012). The minimum magnitudes of verified earthquakes were 6.0 and 6.8 in this study and in the work of Triyoso and Shimazaki (2012), respectively. A comparison of the corrected strain-rate distribution suggests that the different coupling models had a large impact on the results.

Prospective performance of the forecast models
In the previous subsection, most earthquakes occurred before the geodetic data period because the data periods for seismicity and geodetic data were 599 to 2020 and 2005 to 2009, respectively. Although we assume time independence of seismicity, it remains unclear whether the validation results are applicable to future earthquakes. Triyoso and Shimazaki (2012) suggested that the strain concentration zone observed by GNSS was not a steady long-term feature, but a temporal one caused by past large earthquakes. A simple simulation of postseismic deformation due to viscoelastic relaxation demonstrates significant temporal variation in the strain rates around a fault (e.g., Wang et al. 2021). Therefore, one may doubt our assumption that the strain-rate distribution used is a proxy for the long-term average value of strain accumulation.
We consider the assumption is valid as the first order because of the following two reasons, although we need to clarify how significant the distribution of the strain rate changes with the longtime scale in future studies. First, the GNSS velocity observed over the past two decades in southwest Japan is rather stable except for in the period of large slow slip events and postseismic deformation after large earthquakes (Nishimura et al. 2014). Meneses-Gutierrez et al. (2022) also suggest that crustal deformation rates in northeastern Japan were relatively constant over a centennial time scale before the 2011 Tohoku-oki earthquake using conventional geodetic data.  Nishimura et al. (2006) Second, the simulation of postseismic deformation suggests that the postseismic strain rates included in the used data are small in most parts of the analyzed area. We calculated postseismic velocity due to viscoelastic relaxation of 12 large earthquakes (Table 3) in our GNSS data period using the PSGRN/PSCMP software (Wang et al. 2006). These earthquakes whose fault models were studied in previous studies occurred from 1891 to 2005. The assumed viscoelastic structure comprises a 35-kmthick elastic layer overlying a viscoelastic half-space with a viscosity of 10 19 Pas. The predicted horizontal velocity in 2007 shows large-scale trenchward movements up to 4 mm/year in southwest Japan (Fig. 6a). Postseismic strain rates are smaller than 5 × 10 -8 /year in the inland of southwest Japan except for the region along the Pacific coast where large contraction rates of up to 2 × 10 -7 /year are calculated. In addition, it is notable that the predicted postseismic displacement has a long wavelength and can be explained by coupling and slip on the subducting plate interface. It is likely that the postseismic deformation is removed as the elastic deformation due to interplate coupling in the correction process of GNSS velocity. As such, we used earthquakes that occurred only after the geodetic data period to examine the prospective performance of the forecast models. Only eight M ≥ 6 crustal earthquakes occurred in the area analyzed between 2010 and 2020. Three of them were foreshocks and the mainshock of the 2016 Kumamoto earthquake (Kato et al., 2016). Because a small number of tested earthquakes lead to uncertainty and instability in the statistical test, we lowered the magnitude threshold to examine a total of 51 M ≥ 5 earthquakes with a depth of ≤ 20 km, as shown in Fig. 4.
The VARI3 and UNI3 models forecast approximately 21 M ≥ 5 earthquakes over 11 years, with a total moment of ~ 1.3 × 10 19 Nm. These numbers and moments are approximately one-third of the actual values (i.e., 51 and  3.42 × 10 19 Nm), which are mainly attributed to clustered seismic activity of the M7.3 Kumamoto earthquake. Statistical scores are listed in Table 4. The probability gains per earthquake were 1.15 and 1.33 for the VARI3 and UNI3 models, respectively. The error curves in the Molchan diagram show the concentration of epicenters and earthquake moments in both models (Fig. 7a, b). They also indicate that the UNI3 model had better predictive ability than VARI3, as they were close to the 95% confidence limits of the no-correlation model. We also tested our models with the declustered earthquake catalog by removing clustered seismicity including aftershocks. We applied the window method (Gardner and Knopoff 1974)  The temporal windows of declustering are the same as Gardner and Knopoff (1974), but the spatial ones are given by Eq. (11). The number and moments of earthquakes in the declustered catalog are 31 and 2.87 × 10 19 Nm, respectively, which are closer to those predicted than those in the original catalog. The probability gains per earthquake and ASS for the declustered catalog show less predictive capacity than those for the original catalog (Table 4). However, these scores suggest that our models still provide better forecasts than the reference model and the error curves in the Molchan diagram show the concentration of epicenters and earthquake moments in geodetic high-strain-rate cells (Figs. 7c,d). Although the declustering of the earthquakes in the retrospective test reduces the number of earthquakes from 151 to 137, the difference of statistical scores is smaller than that in the prospective test.

Implications of forecast model validation
In both retrospective and prospective examinations, the statistical scores for the moment distribution of past earthquakes were better than those for the number of earthquakes (Table 2), which means that our models forecast the source areas of earthquakes better than their epicenters. It is also notable that the seismic moment of past earthquakes is more concentrated in geodetic highstrain-rate cells than the number of past earthquakes in the Molchan diagram (Figs. 5 and 7). These may be reasonable because geodetic measurements observe the interseismic strain that accumulates around the earthquake source region, not only its epicenter. Because the dimensions of the cell used in this study (0.2°) are approximately equal to the fault length of an M7.2 earthquake, the source area of a large earthquake must occupy multiple cells. Thus, using only the epicenter of large earthquakes can bias the validation test for a high-resolution forecast model with a small cell size. It may be apparently unconcordant with previous studies that seismicity rates, rather than seismic moment rates, are geodetic moment rates (e.g., Kagan 1999;Kreemer et al. 2002), because the seismic moment rates are dominated by the occurrence of rare huge earthquakes in a time window. The large difference in spatiotemporal scale can lead to different conclusions.
Comparisons of error curves between the moments of past earthquakes (red curve) and the geodetic moment rate (blue curve) in the Molchan diagram (Figs. 5 and 7) indicate that the moment distribution of large earthquakes observed geodetically was more concentrated in cells with a high-strain-rate. For example, 64% of the seismic moment was concentrated in the upper quarter of high-strain rate cells for the VARI3 model, whereas 43% of the geodetic moment rate was observed in the same cells (Fig. 5c). Several factors probably contributed to this apparent concentration. First, surface deformation in the interseismic period was more distributed than the coseismic deformation concentrated immediately above the source fault. The elastic dislocation theory demonstrates that interseismic strain rates are distributed around the fault, depending on the down-dip limit of the locked part of the fault (cf. Segall, 2010). Second, the distribution of the GNSS stations was too sparse to capture the real heterogeneous distribution of the interseismic strain rates. For example, a high-strain-rate zone along the coast of the Japan Sea, called the San-in Shear Zone (Fig. 1), was revealed by the velocity data after significant enhancement of the GNSS network in 2002 (Nishimura and Takada 2017). A regional dense GNSS array (e.g., Meneses-Gutierrez and Nishimura 2020; Sagiya et al. 2004) and InSAR observations (e.g., Takada et al. 2018) near active faults are useful for improving the spatial resolution of interseismic deformation. Third, the noise included in the GNSS velocity data increases the "observed" strain rates from those of tectonic deformation. Thus, we speculate that the real tectonic strain rates may be reduced in cells with low strain rates and concentrated in cells with higher strain rates.

Comparison with active faults and seismicity based on forecast models
The HERP started regional evaluations of M ≥ 6.8 crustal earthquakes using the data of active faults and seismicity to publish a 30-year probability of the earthquakes (Earthquake Research Committee 2017). It may be worth comparing our model with the model of the HERP using different datasets. We calculate a 30-year probability of M ≥ 6.8 earthquakes in cells included in the HERP's seven sub-regions of southwest Japan (Fig. 8) using Eqs. (7)-(10) to compare that of the HERP (Table 5). We note that the probability difference among the models is not large except for the northern Chugoku region, where the seismicity-based model gives one order of magnitude higher probability. The HERP seismicity-based model in this region may be biased from the long-term average because several M ≥ 7 earthquakes occurred there over the past 80 years (e.g., Gutscher and Lallemand 1999). On the other hand, the probabilities of our geodetic model and the seismicity-based model are roughly twice as large as those of the active-fault-based model in southern Kyushu. It implies that the earthquake activity on a short timescale of geophysical observations is higher than that in a long timescale of geological observations there.
Previous studies have shown that a model that combines different types of data can improve forecast performance on both global (e.g.,  and regional scales (e.g., Rhoades et al. 2017). GNSS data have the advantage of a relatively uniform spatial resolution throughout Japan. Active fault data are crucial for timedependent forecast models because they provide a history of past earthquakes. We propose that a combination of multidisciplinary data, including GNSS, seismicity, and active faults, is essential to improve RELMs for longterm earthquake forecasting in Japan.

Conclusions and prospects
In this study, we developed time-independent forecast models for crustal earthquakes using GNSS data from southwest Japan. In the forecast models, optimally smoothed GNSS strain rates from 2005 to 2009 were employed after removing the elastic deformation attributed to interplate coupling on the megathrust fault along the Nankai Trough. Then, we converted the strain rates to geodetic moment rates using three equations proposed in previous studies. Finally, we constructed uniform and N a n k a i T r o u g h  Shikoku (G) 10* 10* 13.0 11.8 variable structure models for different values of seismogenic thickness and rigidity in the conversion equations. The derived geodetic moment rate was ~ 5-8 times larger than the seismic moment rate of historical earthquakes during 1586-2020. Therefore, we introduced an empirical factor to calibrate the geodetic moment rate with the seismic moment rate. Finally, we calculated the probability of M ≥ 6 crustal earthquakes over 30 years, assuming a stationary Poisson process and a truncated Gutenberg-Richter law with values of 0.9 and 8.0 as the uniform b value and maximum earthquake magnitude, respectively. The estimated probability exceeded 1% in more than half of the analyzed area and exhibited significant a spatial variation of two orders of magnitude. Furthermore, we validated the forecast models using retrospective and prospective earthquake catalogs with several score indexes and the Molchan diagram. These analyses suggested that the models that use the absolute value of principal horizontal strain rates and areal strain rates better forecast the real seismicity than the models that use the maximum shear strain rates. Models using the absolute value of only the principal horizontal strain rates exhibited the same predictive ability. There was no substantial difference in the forecast performance between models using uniform or variable seismogenic thickness and rigidity. All models explained the real seismicity significantly better than the uniform earthquake distribution model. This study demonstrates that RELMs based on GNSS data can provide enhanced preparedness for earthquake disasters.
Although parameters b , M max , and β , which control the earthquake probability, were assumed to be spatially uniform in this study, they clearly exhibit regional variation in nature. Many studies (e.g., Hirose and Maeda 2011) have estimated the heterogeneous b-value distribution of the Gutenberg-Richter law for the Japanese Islands. However, it remains practically difficult to estimate stable values because they can change not only spatially but also temporally (e.g., Ogata et al. 2018). Spatial variation of M max in Japan are proposed (e.g., Kudo et al. 2009), and it is probably related to seismogenic thickness H . No significant improvement in the performance of the VARI models relative to the UNI models may imply the importance of regional variation of these parameters as well as seismogenic thickness and rigidity. A small β value is expected in a region with ongoing inelastic deformation, including active folding and fault creep. It is also difficult to distinguish between inelastic strain and elastic strain in the observed insufficient data, although a few methods have been proposed (e.g., Noda and Matsu'ura 2010). Estimating these variable parameters is a challenge for future research.
Abbreviations ASS: Area skill score; GNSS: Global Navigation Satellite System; HERP: Headquarters for Earthquake Research Promotion; JMA: Japan Meteorological Agency; J-SHIS: Japan Seismic Hazard Information Station; RELM: Regional Earthquake Likelihood Model.