Efficiency of earthquake forecast models based on earth tidal correlation with background seismicity along the Tonga–Kermadec trench

The correlation between Earth’s tides and background seismicity has been suggested to become stronger before great earthquakes and weaker after. However, previous studies have only retrospectively analyzed this correlation after individual large earthquakes; it thus remains vague (i) whether such variations might be expected preceding future large earthquakes, and (ii) the strength of the tidal correlation during interseismic periods. Therefore, we retrospectively investigated whether significant temporal variations of the tidal correlation precede large interplate earthquakes along the Tonga–Kermadec trench, where Mw 7-class earthquakes frequently occurred from 1977 to 31 December 2020. We evaluated a forecast model based on the temporal variations of the tidal correlation via Molchan’s error diagram, using the tidal correlation value itself as well as its rate of change as threshold values. For Mw ≥ 7.0 earthquakes, this model was as ineffective as random guessing. For Mw ≥ 6.5, 6.0, or 5.5 earthquakes, the forecast model performed better than random guessing in some cases, but even the best forecast only had a probability gain of about 1.7. Therefore, the practicality of this model alone is poor, at least in this region. These results suggest that changes of the tidal correlation are not reliable indicators of large earthquakes along the Tonga–Kermadec trench.


Introduction
The Schuster test (Schuster 1897) is commonly used to evaluate the relation between seismicity and Earth's tides (e.g., Tsuruoka et al. 1995;Emter 1997); the smaller the index p-value becomes, the higher the correlation between tides and seismicity (see Sect. 3 for a detailed definition). Tanaka (2012) found that the p-value around the source area decreased before and increased after the 2011 Tohoku earthquake (M w 9.0). Similar p-value variations were also documented for the 2004 Sumatra earthquake (M w 9.0) and its larger aftershocks (M w 8.5 and 8.6; Tanaka 2010), and the 19 December 1982 South Tonga earthquake (M w 7.5; Tanaka et al. 2002b). Furthermore, Tanaka (2012) reported that the precursory time scaled with the mainshock magnitude; thus, the p-value of the tidal correlation has been considered a promising tool for forecasting large earthquakes. However, each of these studies only reported retrospective analyses of individual large earthquakes, and it remains vague whether a decrease in p-value variations might precede future large earthquakes and whether p-values remain high during interseismic periods.
In contrast, other researchers reported negative correlations between spatiotemporal p-value variations and large earthquakes (e.g., Heaton 1982;Wang and Shearer 2015). After analyzing a large number of M ≥ 3.0 events off the Tohoku coast of Japan, Wang and Shearer (2015) argued that the low p-value preceding the 2011 Tohoku earthquake reported by Tanaka (2012) was due to chance because the p-value was sensitive to their selection of spatiotemporal grids.
The Tonga-Kermadec trench is an area of very high seismicity and the site of frequent M 7-class earthquakes. Since the 1982 South Tonga earthquake (Tanaka et al. 2002b), large interplate earthquakes of M w ≥ 7.0 have occurred often in this area (Fig. 1). In our previous study (Hirose et al. 2019), we investigated the correlation between Earth's tides and 661 interplate earthquakes of M w ≥ 5.5 that occurred during a 40-year span along the Tonga-Kermadec trench and found that tidal normal stresses controlled the occurrence of earthquakes. However, we did not analyze temporal variations of the p-value. To address this issue, we here report a retrospective investigation of whether p-value variations preceded any of 729 interplate earthquakes of M w ≥ 5.5 that occurred during a 44-year span along the Tonga-Kermadec trench. Fig. 1 Interplate (thrust) earthquakes (focal depth ≤ 70 km, M w ≥ 5.5) that occurred along the Tonga-Kermadec trench during the period 1 January 1977 to 31 December 2020 (N = 729). a Epicenter distribution. The purple lines denote plate boundaries at Earth's surface (Bird 2003). Colored, dashed contours indicate the depth (km) to the upper surface of the subducting Pacific slab (Hayes 2018). White and black arrows indicate the movement of the Pacific plate relative to the Tonga plate (TO) and Kermadec plate (KE), respectively. b Time-space plot of earthquakes. c Time-series plot of earthquake magnitudes (red bars indicate events with M w ≥ 7.0) and the cumulative number of earthquakes (black trend). Events with letters * and + indicate events excluded from the forecast target because of insufficient data to calculate a p-value. Events with letters A to J indicate the forecast targets. These events with letters correspond to red circles in d. d Relation between phase angle and tidal normal stress on the fault at the time of earthquake occurrence; red symbols indicate events with M w ≥ 7.0. (e) Frequency distribution of tidal phase angles estimated using the normal stress Δσ on the fault. The

Data
We extracted 729 interplate-type earthquakes (strike, 150-230°; rake angle, 55-125°; 0-70 km depth; M w ≥ 5.5) from the Global Centroid Moment Tensor (GCMT) catalog Ekström et al. 2012) that occurred from 1 January 1977 through 31 December 2020 at 15-35° S latitude along the Tonga-Kermadec trench (Fig. 1). Our event selection criterion was M w ≥ 5.5 because this magnitude is above the detection threshold based on the magnitude-frequency distribution (Hirose et al. 2019). We used the phase angle of tidal normal stress (Δσ) at the origin time on the assumed fault plane based on the extracted GCMT solution because seismicity is known to depend on Δσ in this region (Tanaka et al. 2002b;Hirose et al. 2019).

Theoretical tidal response
We calculated theoretical tidal responses on the faults using "TidalStrain.2" software (https:// mri-2. mri-jma. go. jp/ owncl oud/s/ tjqx7 HfK8b D3KQf ), which separately estimates the effects of solid and ocean tidal loadings on faults. The solid tidal loading effect is evaluated by multiplying the tide-generating potential, calculated directly from the celestial coordinates of the Sun and Moon (Nakai 1979), by the eigenfunctions of spherical wavenumber n = 2 calculated for a tidal surface boundary condition. The ocean tidal loading effect is calculated by convolution of the Green's functions for surface vertical point loading problem, with the loading mass distribution given by spherical ocean tide models (NAO.99b, Matsumoto et al. 2000;and NAO.99L, Takanezawa et al. 2001) for each of 21 tidal constituents. We used the PREM Earth model (Dziewonski and Anderson 1981) with the topmost 3 km oceanic layer replaced by a solid layer (V p = 5 km/s, V s = 2.6 km/s, ρ = 2.6 g/cm 3 ) (Tsuruoka et al. 1995). TidalStrain.2 sums the temporal variations of the solid and oceanic tidal loading effects for six independent components of the strain tensor (Ozawa 1974), as estimated at the hypocenter of each event.
Then, it converts them to temporal variations of the normal stress (Δσ) on the assumed fault plane, based on the GCMT solution extracted in Sect. 2. We used a sampling interval of three minutes. See Hirose et al. (2019) for further details of this method. We defined dilatation and compression as positive and negative normal stress, respectively; consequently, positive Δσ values promote fault slip. We assigned phase angles of − 180° and + 180° to the minimum tidal stresses before and after an event, respectively, and 0° to the maximum tidal stress that occurred between these two minima. The phase angle at the earthquake occurrence time was estimated by linear interpolation over the time interval between − 180º and 0° or between 0° and + 180°.
The relation between tidal phase angle and tidal normal stress on the fault at the time of earthquake occurrence is illustrated in Fig. 1d. Tidal normal stress tended to be positive for phase angles of − 90° to + 90°, and negative otherwise. In the frequency distribution of tidal phase angles, the peak of the distribution was between − 20° and 0° (Fig. 1e), indicating that interplate earthquakes of M w ≥ 5.5 tend to occur at specific tidal phase angles.
We note that a minor bug in the calculation of ocean tidal loading effects in the TidalStrain software (used in Hirose et al. 2019) was fixed in TidalStrain.2 (see manual for details). Therefore, we recalculated the tidal phase angles of the normal stress for the 661 events from Hirose et al. (2019) using TidalStrain.2: 598 events were consistent, 50 differed by ± 1°, and 13 by ± 2°.

Schuster test
The possibility of earthquakes occurring at a particular tidal phase angle is commonly evaluated using the p-value, calculated by Schuster (1897) as: where N is the total number of earthquakes, ψ i is the phase angle of the ith earthquake, and D indicates the final distance from the origin after a two-dimensional random walk of N steps. Equation (1) is the complementary cumulative distribution function of the Rayleigh distribution and corresponds to the probability that the magnitude of the vector sum of a random set of earthquake phase angles will be greater than D (Hirose et al. 2019); we note that this approximation is sufficient only for N > 10 (Heaton 1975). Thus, the p-value (ranging from 0 to 100%) represents the significance level for rejecting the null hypothesis that earthquakes occur randomly with respect to the tidal phase angle, such that the confidence in rejecting the null hypothesis becomes greater as the p-value becomes smaller. In general, p-values of 5% (Emter 1997;Tanaka et al. 2002aTanaka et al. , 2006 or 10% ) are used to judge tidal correlations. For the selected earthquakes, the p-value is 3.93% (Fig. 1e), indicating that the earthquakes were, at least partially, tidally triggered, as pointed out in our previous study (Hirose et al. 2019).
In Eqs. (1) and (2), the p-value depends on the number of earthquakes, regardless of the tidal phase angle distribution (Tsuruoka and Ohtake 2002). Therefore, it is best to use a window of a fixed number of events rather than a fixed time interval when discussing the temporal variation of the p-value. Thus, we used fixed number windows of 30, 40, and 50 events with their ending time shifted by one day; for these number windows, the median periods are 1.75, 2.31, and 2.94 years, respectively. Figures 2 and  3 show temporal variations of the p-values. The plotting point of the p-value corresponds to the right-end of the number window. In particular, Fig. 3 shows variations of the p-value for the 1000 days preceding earthquakes of M w ≥ 7.0. We investigated whether these precursory decreases of p-values precede earthquakes of specific magnitude using Molchan's error diagram (Molchan 1997).

Molchan's error diagram
To evaluate the performance of our forecast model based on Schuster's p-value, we used Molchan's error diagram (Molchan 1997), which consists of two parameters: the alarmed fraction τ (= cumulative alarmed periods (duplicated period is excluded)/studied time interval = alarm rate AR/probability gain PG) and the failure-to-predict (miss) rate ν (= number of missed targets/known number of events = 1 − AR). This diagram allows a visual comparison of model performance compared to random guessing; an ideal forecast plots near the origin, indicating maximum success ( ν → 0 ) at minimum cost ( τ → 0 ). The alarm rate AR is the number of target earthquakes in alarmed space-time divided by the total number of target earthquakes, and the probability gain PG is the occurrence rate of target earthquakes in alarmed spacetime divided by the background occurrence rate. PG can also be used to evaluate the performance of a forecast model; when PG < 1, the forecast model is inferior to the uniform Poisson model. The smaller the forecast error e(τ , ν) = τ + ν , the better the model, and the combination of parameters that maximizes Peirce's skill score (SS p ; Peirce 1884; Talbi et al. 2013): is determined to be the optimal forecast model.

Confidence interval for a random guess
In the case where target earthquakes occur n times during the study period T, i.e., the average occurrence rate = n/T , the probability of obtaining k or more hits by chance over the entire duration of the alarmed periods ( τ T ) is described by the binomial distribution as: where n C i indicates the number of possible ways to select i objects from a set of n objects without repetition or order. By searching for values of τ such that α becomes smaller than 0.01 for a given value of ν = 1 − k/n , a 99% confidence interval on τ can be determined by randomly guessing.

Combination of parameters in the forecast model
Temporal variations of p-values were determined for fixed number windows (N) of 30, 40, or 50 events with their ending time shifted by one day (Fig. 2). As the threshold p-value in the forecast model, we used both the p-value itself and its rate of change. The threshold p-value ( p th ) is 1, 2, …, 100% or less, and the threshold rate of change of the p-value ( �p th = log 10 (p 2 /p 1 ) ) is − 0.01, − 0.02, …, − 1 or more negative, where p 1 is the p-value T f days before p 2 , and T f was set to 360, 720, or 1080 days. The alarm period ( T a ) was 30, 180, or 360 days after satisfying the threshold, and the target earthquakes were of M w 7.0, 6.5, 6.0, 5.5, or larger. We assumed the alarm area to be 15-35° S along the Tonga-Kermadec trench so that all earthquakes extracted in Sect. 2 were included (Fig. 1a).   Figure 4 shows the case for which SS P is the largest for each target magnitude when the forecast model is based on a threshold p-value. For target earthquakes of magnitude M w ≥ 7.0 (Fig. 4a), the best model performance (maximum SS P = 0.30, red star) was obtained for an alarm issued for T a = 360 days after the p-value based on the previous N = 40 events decreased below 14%. However, because this forecast is within the 99% confidence interval of a random guess, it cannot be considered a valid forecast; indeed, the maximum SS P value obtained for a random guess was 0.32 (blue diamond in Fig. 4a). In each plot, the blue line shows the boundary line equivalent to the maximum SS P value obtained by randomly guessing; even if a forecast deviates from the 99% confidence interval of a random guess (e.g., points near the bottom-right of Fig. 4a), it should be judged ineffective if it plots above the blue line.

Performance of the earthquake forecast model based on tidal information
For target earthquakes of magnitude M w ≥ 6.5 (Fig. 4b), the maximum SS P value obtained for the forecasts was smaller than that of a random guess; again, these forecasts cannot be considered valid. For target earthquakes of magnitude M w ≥ 6.0 or 5.5 (Fig. 4c, d, respectively), the maximum SS P values of the forecasts were only slightly larger than the maximum SS P values of random guesses. In addition, because the threshold p-values in these cases were around 30%, the forecasts were based on information not closely related to the tidal correlation; even though these forecasts out-performed random guesses, we cannot say that they are based on a correlation between background seismicity and tides. Figure 5 is the same as Fig. 4, but for forecasts based on the rate of change of the p-value. For target earthquakes of M w ≥ 7.0 (Fig. 5a), the maximum SS P values of the forecasts were always smaller than those of random guesses, and the forecasts cannot be considered valid. For target earthquakes of M w ≥ 6.5, 6.0, and 5.5 (Fig. 5b-d), the maximum SS P values of the forecasts were slightly larger than those of random guesses. The maximum SS P difference ( SS O max p − SS R max p , where SS Omax p and SS Rmax p are the SS P values of the forecast and random guessing, respectively) was obtained for the forecast that earthquakes of M w ≥ 6.5 will occur within T a = 30 days after p th decreases to values more negative than − 0.14 for T f = 1080 days (i.e., the p-value decreases to below ~ 72% (= 10 -0.14 ) of the value approximately 3 years earlier) for a number window of N = 30 events. However, even in this best case, the PG of the forecast model was approximately 1.7, far less than those of forecasts based on foreshocks (well over 100; Nakatani 2020; Hirose et al. 2021). In addition, because the alarm period occupied almost half of the entire studied time interval ( τ = 0.44 , orange periods in Fig. 2), the practicality of this forecast model is poor.
In summary, forecasts based on tidal correlation information may yield slightly better results than random guesses for certain parameter values, but, because even the best forecast only has a probability gain of about 1.7, these forecasts are not practical. These results suggest that characteristic precursory p-value variations are not Molchan's error diagrams as in Fig. 4, but for forecasts based on the rate of change of the p-value always associated with large earthquakes, indicating that precursory background seismicity is not always strongly correlated with Earth's tides. Note that we also performed the same investigation for the tidal shear stress and reached the same conclusion because a probability gain was 1.4 even in the best forecast.

Discussion and future work
Precursory variations of p-values were reported off Sumatra for the 2004 earthquake (M w 9.0; Tanaka 2010) and off Tohoku for the 2011 earthquake (M w 9.0; Tanaka 2012); these analyses respectively covered a large area of 1500 km × 500 km and a limited area of 200 km × 200 km near the mainshock epicenter. Even for the same megathrust earthquake, precursory p-value variations may be localized, and might not necessarily appear as a precursory phenomenon over a wide enough area. In fact, Wang and Shearer (2015) argued that the low p-value before the 2011 Tohoku earthquake was sensitive to the selection of spatiotemporal grids. The spatial extent over which background seismicity is highly correlated with tides thus appears to be a complex issue. Here, we investigated the correlation between background seismicity and tides along the Tonga-Kermadec trench, which extends for about 2000 km. In the future, it will be necessary to investigate this correlation over a narrower area.
In the cases of the Sumatra and Tohoku earthquakes, seismicity was correlated with tidal shear stress rather than tidal normal stress (Tanaka 2010(Tanaka , 2012. However, on the plate boundary along the Tonga-Kermadec trench, tidal normal stress controls the triggering of earthquakes and the apparent friction coefficient is relatively large (Tanaka et al. 2002b;Hirose et al. 2019); Hirose et al. (2019) includes a discussion of areas and tectonic backgrounds that are particularly sensitive to tidal normal stress. Precursory p-value variations may be sensitive to interplate friction. Indeed, the difference between the magnitudes of background seismicity (M w ≥ 5.0) and those of target earthquakes (M w 7-8) in this study is smaller than that off Sumatra and off Tohoku (M w ≥ 5.0 and M w 8-9). Accordingly, we might not be able to reliably determine precursory p-value variations. Further studies should investigate the use of smaller earthquakes after more seismic data have been accumulated.
Along the Tonga-Kermadec trench, tidal fluctuations with a cycle of about 12 h are predominant due to the M2 tidal constituent (Hirose et al. 2019). Therefore, we herein investigated a forecast model based on p-values due to the semi-diurnal tide mainly. However, Ide et al. (2016) investigated the relation between the spring tide (cycle ~ 15 days) and 11,397 events of M w ≥ 5.5 that occurred globally during 1976-2015; they found that 75% (= 9/12) of earthquakes of M w ≥ 8.2 occurred on days when the daily maximum tidal shear stress amplitude was in the top third of that during the preceding 15 days. This tendency is not obvious for earthquakes of M w < 8.2, and we have not verified their results for larger earthquakes because a M w ≥ 8.2 earthquake did not occur within this study area/period (Fig. 1).
According to Scholz and Small (1997), the subduction of large seamounts generally increases the normal stress across the subduction interface, thereby strengthening seismic coupling and increasing earthquake recurrence intervals. Thus, they suggested that a large M w 8-class earthquake might be likely to occur in the seismically quiescent area around 26° S, where the Louisville seamount chain is being subducted. Unfortunately, we cannot use p-value information to forecast M w 8-class earthquakes in that area due to the lack of background seismicity needed to reliably estimate the p-value.
The peak of the tidal phase angle distribution is between − 20° and 0° (Fig. 1e), and the frequency of background seismicity for those phase angles is about 1.5 times the average (= 8.09/5.55). If an alarm were automatically issued during periods when the tidal phase angle was between − 20° and 0°, the alarm fraction τ would be 0.06 (= 20°/360°). However, over the study period, no earthquakes of M w ≥ 7.0 occurred at such tidal phase angles (Fig. 1d). For target earthquakes of M w ≥ 6.5, 6.0, or 5.5, the case of M w ≥ 6.0 had the highest probability gain (1.6), and miss rate and SS P values of 0.91 and 0.03, respectively. The probability gains of < 2 obtained herein (best case, 1.7; Fig. 5b) are far smaller than that (> 100) of a previous forecast based on foreshocks (Nakatani 2020; Hirose et al. 2021). Therefore, it is necessary to develop a comprehensive forecast model (Aki 1981) that incorporates various precursors, such as foreshocks and seismic quiescence (Nakatani 2020).

Summary
We investigated the efficiency of an earthquake forecast model based on tidal correlation information using Molchan's error diagram for interplate earthquakes along the Tonga-Kermadec trench. For target earthquakes of magnitude M w ≥ 7.0, the forecasts were as ineffective as random guessing. For target earthquakes of M w ≥ 6.5, 6.0, or 5.5, the forecasts sometimes showed slightly better results than random guessing, but even the best forecast only had a probability gain of about 1.7. Therefore, the practicality of this forecast model alone is poor for seismicity in this region. These results suggest that characteristic changes in the p-value of the correlation between