Changes in seismicity pattern due to the 2016 Kumamoto earthquake sequence and implications for improving the foreshock traffic-light system

Crustal deformation due to the 2016 earthquake sequence in Kumamoto, Japan, that culminated in a preceding earthquake of magnitude M6.5 and a subsequent M7.3 earthquake 28 hours later, caused stress perturbation on and around the causative Futagawa-Hinagu fault zone. Monitoring changes in seismicity pattern along this zone plays a role in understanding the process before and after major earthquakes. For this purpose, stress-dependent laws in statistical seismology can be used: the Gutenberg-Richter frequency-size law and the Omori-Utsu aftershock-decay law. We review the results obtained by using these laws in previous studies to show a zone of high stress near the eventual epicenters of the M6.5 and M7.3 earthquakes before the start of the Kumamoto sequence, and after it, showing a decreasing trend in stress along the Futagawa-Hinagu fault zone. Detailed analysis suggests aseismic slips along the causative fault zone. The aseismic preslip locally reduced stress just prior to the M7.3 earthquake near its epicenter. Recently, a system was proposed by Gulia and Wiemer (2019) that utilizes the Gutenberg-Richter frequency-size law to judge, immediately after a large earthquake, whether it was the mainshock or a foreshock to a future event. Based on the reviewed results and our new results, further research that takes into account the spatial variation of frequency-size distribution, allowing the exploration of the possibility of a local preslip of a future nearby earthquake, is needed to improve this system.


Introduction
The 2016 Kumamoto earthquake sequence (Fig. 1) that caused serious damage to Kumamoto, in central Kyushu, Japan (Hashimoto et al., 2017), began with an earthquake of magnitude (M) 6.5 on April 14, 2016, at 21:26, on the northern section of the Hinagu fault zone (FZ).
It was followed by numerous events, including a M6.4 earthquake on April 15, at 00:03. These M6.5-class earthquakes, which were considered as foreshocks, involved a predominantly rightlateral strike slip. Eventually, on April 16, at 01:25, a M7.3 mainshock occurred on the Futagawa FZ that conjugated with the northern section of the Hinagu FZ, revealing a similar mechanism of fault motion as the M6.5-class foreshocks. The M7.3 mainshock started at the westernmost part of the Futagawa FZ, and its rupture propagated toward the east and terminated around Mt. Aso (Yagi et al., 2016). The Earthquake Research Committee (ERC) concluded that the Futagawa FZ and the northern section of the Hinagu FZ had been activated (ERC, 2016a). Over the first three years since the M7.3 mainshock (Fig. 1b), more than 27,000 aftershocks with M ≥ 1 occurred, including more than 100 events with M ≥ 4. The broader context of the Kumamoto earthquakes is that the Futagawa-Hinagu FZ is encompassed by the Beppu-Shimabara graben, a geological formation that runs across the middle of Kyushu, from Beppu Bay in the east to the Shimabara Peninsula in the west (Inset of Fig. 1). This is understood as being the result of crustal deformation caused by the rifting and spreading of the Okinawa Trough, which is viewed as a continuation of the Beppu-6 Shimabara graben (Tada 1985).
Monitoring crustal deformation and seismicity along the causative faults plays a role in understanding the process before and after major earthquakes. Numerical simulation and rock mechanics experiments suggested the occurrence of pre-and post-seismic slip transients. Some studies implied long-term precursory slip acceleration while other studies indicated a nucleation process that initially involves stable and slow rupture growth (Hori and Miyazaki, 2010;Yoshida and Kato, 2003;Dieterich, 1979;Lapusta et al., 2000). However, it is generally accepted that a precursory slip before a large earthquake is not always observed. This may be because this phenomenon does not always emerge in nature and/or due to difficulty in observing a small slip. On the other hand, a post-seismic slip was observed for many earthquakes. In the aftermath of the Kumamoto earthquakes, slow deformation caused by a post-seismic slip of the causative Futagawa-Hinagu FZ, together with flow in the Earth's interior, has been detected by data from the GEONET: GNSS Earth Observation Network System (Moore et al., 2017;Pollitz et al., 2017).
Recognizing earthquakes as foreshocks in real time would provide a valuable forecasting capability. However, decades of work have failed to establish a robust feature of individual foreshocks that distinguishes them from mainshocks or their aftershocks (Dascher-Cousineau et al., 2020). Foreshocks may be distinctive due to some precursory loading process or influence from the locked zone of the subsequent mainshock, neither of which will exist in the aftershock sequence 7 (Brodsky and Lay, 2014). A fundamental question regarding the Kumamoto sequence is whether or not any indication of an impending large earthquake had appeared before the M7.3 mainshock. The Japan Meteorological Agency (JMA) held a press conference after the M6.5 earthquake to warn of the possibility of large aftershocks with further damage. However, the JMA did not considered the occurrence of large earthquakes, according to the ERC (1998) protocol (see also ERC, 2016b), so no information was made public on the increased probability of M7 or larger earthquakes.
In a recent study, Gulia and Wiemer (2019) proposed a Foreshock Traffic-Light System (FTLS) that relies on abrupt changes in b-values of the Gutenberg-Richter (GR) law (Gutenberg and Richter, 1944), relative to background values. The GR law is given as log 10 N = a-bM, where a and b are constants, and N is the cumulative number of earthquakes with a magnitude larger than or equal to M. The approach utilizes high-resolution earthquake catalogs to monitor localized regions around the largest earthquakes and to distinguish foreshock sequences (reduced b-values) from aftershock sequences (increased b-values) (Gulia et al., 2020;Dascher-Cousineau et al., 2020). Evidence supporting the FTLS hypothesis emerged from investigating a timeseries of two sequences: the Norcia sequence in Italy and the Kumamoto sequence (Gulia and Wiemer, 2019). Their results of the Kumamoto case showed that b-values in the time interval between the M6.5 and M7.3 events first dropped from the background level and then increased over time, and that the b-values at the intermediate and ending times were below the background level, triggering a red traffic light (Table   8 1 of Gulia and Wiemer, 2019) and suggesting that seismicity in between was not considered as an aftershock sequence, rather as a foreshock one. Once the M7.3 mainshock occurred, b-values of the subsequent events increased strongly by 20-40%, triggering a green traffic light, suggesting that these events were aftershocks.
We first review our previous studies Nanjo and Yoshida, 2017) that used the b-value of the GR law and the p-value of the Omori-Utsu (OU) law (Omori, 1894;Utsu, 1961) to explore the changes in seismicity before and after the Kumamoto mainshock. The OU law is given as λ = k(c + t) -p , where t is the time since the occurrence of a mainshock, λ is the number of aftershocks per unit time at t with M greater than or equal to a cutoff magnitude, and c, k, and p are constants. The b-value is inversely proportional to differential stress (Scholz, 1968;2015), and the p-value can be used to infer stressing history (Dieterich, 1994). The main conclusion is that localand time-dependent variation in b-and p-values was linked with changes in stress state associated with aseismic slips and a coseismic one along the causative Futagawa-Hinagu FZ. Based on this review, we revisited the applicability of FTLS to the Kumamoto case and discuss the possibility of improving it.

A review of studies on b-and p-values associated with the Kumamoto earthquakes
Understanding the stress state is of key importance for earthquake forecasting and hazard assessment. However, it is rather difficult to measure stress directly underground. Statistical parameters of seismicity might provide a useful way to monitor stress evolution on faults. Therefore, monitoring changes in seismicity pattern along a FZ plays an important role in understanding the process before and after major earthquakes. A set of three papers (Nanjo et a., , 2019Nanjo and Yoshida, 201) conducted a systematic study on the seismicity associated with the entire process of the Kumamoto sequence.
Using earthquake data since 2000 from the JMA catalog,  showed that seismicity before the start of the Kumamoto sequence revealed a zone of low b-values near the eventual epicenters of the M7.3 mainshock and two M6.5-class foreshocks (Fig. 2a). The precursory duration to the M7.3 mainshock was found to be unknown based on observation of an approximately constant b-value over the past 16 years . A decrease in b tracking stress buildup, as expected from a laboratory experiment (Scholz, 1968;Lei, 2003), might not be observable for the decade-scale observation of active faults, because it was too short to observe such a decrease in b for the Kumamoto case. The mechanism of stress buildup within a FZ is uncertain, but one hypothesis is that a steady slip of the deeper continuation of faults that does not produce earthquakes, but still involves motions across the fault, forces the upper crust around the faults to deform, and thus concentrate stress. The stress field remains difficult to measure directly (Narteau et al., 2009), but it would be helpful if some indirect evidence exists, such as stress inversion by a focal mechanism, or cases of stress buildup. A specific case includes thermally activated creep on major strike-slip faults such as the San Andreas Fault in California (Turcotte et al., 1980;Turcotte and Schubert, 2002).
These faults extend deep into the lithosphere, where they are likely to behave plastically, displaying a steady-state creep in deep fault zones.
Assuming that the M6.5 foreshock was a mainshock, Nanjo and Yoshida (2017) showed that the OU law with p > 1 was correlated with seismic activity for the entire area during the period of 1.16 days between the M6.5 and M7.3 earthquakes. The same analysis was conducted to obtain p > 1 for the northern area (above 32.72N) and p < 1 for the southern area (below 32.72N), indicating rapid stress relaxation in the northern area, rather than in the southern one. The b-values significantly increased over time for the northern area ( Fig. 2b), supported by the Utsu test (Utsu, 1999;Schorlemmer et al., 2004). An increase in b over time was observed for the southern area, but this increase was insignificant. The authors suggested that stress around the Futagawa FZ had begun to decrease preceding the M7.3 mainshock because they found a large p-value (p > 1) and a time-dependent increase in b for the northern area. Nanjo and Yoshida (2017) also pointed out that earthquakes occurred near the hypocenter of the M7.3 mainshock in addition to most seismicity around the northern section of the Hinagu FZ.
These events occurred later (0.98-0 days preceding the M7.3 mainshock). One explanation, based on Onaka (1993), may be that they were produced by a nucleation process of the M7.3 mainshock. This process is considered to be promoted by the gradual propagation of an aseismic slip (afterslip) toward the M7.3 hypocenter on and around the Hinagu FZ, interpreted from a detailed observation of seismicity migration starting immediately after the occurrence of the M6.5 earthquake and terminating at early times after this earthquake . The finding of the occurrence of earthquakes close to the M7.3 hypocenter at later times after the M6.5 earthquake implied the growth of a quasi-static preslip before the M7.3 mainshock. This preslip might have relaxed stress around the Futagawa FZ that hosted the forthcoming M7.3 mainshock. This stress relaxation provided feedback, causing a significantly larger p-value than p = 1 and a significant increase in the b-value ( Fig. 2b) in the northern area, including a part of the Futagawa FZ. This relaxation was also inferred from the result of a stress inversion analysis performed by the National Research Institute for Earth Science and Disaster Resilience (2016) that suggested the reduction in normal stress at the Futagawa FZ, i.e., a decrease in compressive stress on the fault. It was considered that afterslip of the M6.5 earthquake was not the main cause of the observed change in seismicity, but might set up the conditions that allowed this change in seismicity. This is because a detailed observation  showed that the afterslip occurred only on and around the Hinagu FZ while the major slip was limited to an early period after the M6.5 event.

Method and Data
To estimate b-values consistently over space and time (Figures 3 and 4), we employed the EMR (entire-magnitude range) technique (Woessner and Wiemer, 2005), which also simultaneously calculates the completeness magnitude M c , above which all events can be detected by a seismic network. EMR applies the maximum-likelihood method to compute the b-value to events with a magnitude above M c . Uncertainties in b were computed by bootstrapping (Schorlemmer et al., 2003).
If logP b , the logarithm of the probability that the b-values are not different, is equal to or smaller than -1.3 (logP b ≤ -1.3), then the change in b is significant (Utsu, 1999;Schorlemmer et al. 2004).
Spatiotemporal changes in b are known to reflect a state of stress in the Earth's crust (Narteau et al., 2009;Scholz, 2015). Patches with small b-values on active faults have been observed to coincide with locations of subsequent large earthquakes (Nanjo, 2020; 14 Tormann et al., 2013;Wang et al., 2021). Investigations using b-values have been applied to locate asperities and to estimate frictional properties on the plate interface along subduction zones (e.g. Ghosh et al., 2008;Nanjo et al., 2012;Nanjo and Yoshida, 2018;Sobiesiak et al., 2007;Schurr et al., 2014;Tormann et al., 2015). The results of laboratory experiments indicate a systematic decrease in the b-value approaching the time of the entire fracture (Lei, 2003;Scholz, 1968). Foreshocks have been known to show small b-values (Gulia and Wiemer, 2019;Gulia et al., 2020;Dascher-Cousineau et al., 2020).
Our dataset was the earthquake catalog maintained by the JMA. We used about 10 6 earthquakes having M ≥ 0 during the period from January 1, 2000 to March 11, 2019 with depths shallower than 25 km within the study region (Fig. 1). This period is the same as that considered for a previous paper , and subsets of that period in two other papers Nanjo and Yoshida, 2017). A b-value analysis is critically dependent on a robust estimate of completeness of the processed earthquake data. In particular, underestimates of M c lead to systematic underestimates of b-values. We always paid attention to M c when assessing M c at each node and time window. As shown in other studies using the JMA catalog Schorlemmer et al., 2018), typical values for M c are presently 0.5~1 in the Kyushu district that includes the Kumamoto region. Small earthquakes in the early stage after a large event are, in many cases, missing from the earthquake catalogs, as they are ''masked'' by the coda of the large event and overlapped with each other on the seismograms. According to Nanjo and Yoshida (2017) who showed time-dependent decrease of M c during the period between the M6.5 and M7.3 quakes, we examined a M c -time pattern (Fig. 4).

Results and Discussion
A review of previous studies Nanjo and Yoshida, 2017) indicates a systematic investigation into the seismicity before, during, and after the Kumamoto earthquakes that started in April of 2016, by analyzing the b-and p-values (section 2). A brief summary is as follows: (1) the focal areas of the M  6 earthquakes were more stressed than others before the start of the Kumamoto sequence; thereafter, stress generally decreased along the causative faults, as inferred from the spatial distribution of b-and p-values. (2) By using a more detailed analysis, the two statistical parameters suggested aseismic preslip and afterslip along the causative faults. The preslip locally reduced stress just prior to the largest earthquake near its epicenter.
Using the idea of the FTLS, analysis of b-values for seismicity during the time interval between the M6.5 and M7.3 earthquakes (Fig. 2b), in contrast to b-values for seismicity before the M6.5 earthquake (Fig. 2a), implied a yellow FTLS status for the southern area that included a part of the Hinagu FZ and a green FTLS status for the northern area that included a part of the Futagawa FZ hosting the eventual M7.3 earthquake. Here, we examined the inconsistency with the result obtained by Gulia and Wiemer (2019), who showed that a red traffic light was triggered before the M7.3 mainshock for the purpose of improving the FTLS. quake are not different, again a reproduction from Nanjo and Yoshida (2017). A general trend of the time-dependent increase in b in the interval between the M6.5 and M7.3 events, which had been shown by Gulia and Wiemer (2019), was also found in Fig. 2b. However, if we considered the southern and northern areas separately, our analysis of b-values confirmed a yellow FTLS status for the former area and a green FTLS status for the latter area. For the northern area, which includes a part of the Futagawa FZ on which the impending M7.3 mainshock occurred, b-values showed an increase over time, to values exceeding 1.0 (bottom panel of Fig. 2b), significantly larger than the background levels (b = 0.81~0.83  0.05~0.07) found by Gulia and Wiemer (2019) and those (b = 0.6~0.8) in the area of the M6.5 and M7.3 epicenters (Fig. 2a). On the other hand, for the southern area, which includes a part of the northern section of the Hinagu FZ on which the M6.5 foreshock had already occurred, b-values were similar to the background level.
The earthquake samples for b-value estimation in a grid search (Fig. 2)  To eliminate doubts about this bias, we conducted a b-value timeseries analysis (Fig. 3) based on events in two regions that do not overlap with each other: one region (blue) includes only a part of the Hinagu FZ, and the other (red) includes only a part of the Futagawa FZ. In the latter region, an area close to the M7.3 mainshock epicenter was not included, because it was difficult to distinguish between events from the Hinagu FZ and events from the Futagawa FZ. Fig. 4 shows the time dependence of M c , confirming Nanjo and Yoshida (2017), where M c was simultaneously calculated when applying the EMR method for obtaining the b-value. In creating Figs. 3 and 4, we used a moving window approach, whereby the window covered 80 events, and plotted b and M c at the end of the moving window that they represent. For both regions (Fig. 3) Our results for the northern area (indicative of a green FTLS status) suggested a scenario of a reduction in stress induced by growth of aseismic preslip on the Futagawa FZ, as discussed in section 2 (also see Nanjo and Yoshida, 2017), subsequently leading to a much larger chance of a stronger M7.3 event. The results by Gulia and Wiemer (2019)

(indicative of a red FTLS status)
suggested an alternative scenario in which stress on the Futagawa FZ increased, leading to a subsequently much larger chance of a stronger M7.3 event. The discrepancy between our study and the Gulia and Wiemer (2019) study is that the former study investigated spatial variation of b-values, allowing us to find clues of preslip on the nearby Futagawa FZ, which were not found by the latter study. The Kumamoto case might not be a good application of the FTLS.
The picture changes markedly after the M7.3 mainshock: the b-values in the Kumamoto source area that largely fluctuated during the first several months (Fig. 3c), increased significantly, compared with the background (also see section 2 and Nanjo et al., 2019). This is consistent with the findings by Gulia and Wiemer (2019), who suggested that the chance of a subsequent large event was substantially small, although the Kumamoto aftershock sequence included many small events.
No subsequent large event has occurred so far.
The FTLS system is quite new, and some cases, except for the Kumamoto one, would provide support for it. We agree with Gulia et al. (2020) that it is too early to use the current version of the FTLS routinely for making decisions about civil protection and public communications. More extensive tests of sensitivity and robustness are needed, as pointed out by those authors. Based on our studies, it would be very worthwhile to improve the FTLS in ways that would allow spatial variability of b-values to be taken into consideration, and if an indication of preslip appeared, this would trigger a red traffic light. Also, the problem may be a lack of another type of sensitivity, namely only 3 levels of sensitivity, red, yellow and green. Perhaps a wider range of sensitivities would be more useful, but then the name would likely have to change to something other than the FTLS, since traffic lights do not have more than 3 colors.

Conclusion
This paper reviews previous studies on spatial and temporal patterns of seismicity in the Kumamoto region, in Japan, where some large earthquakes took place in 2016, including one with M7.3 2019;Nanjo and Yoshida, 2017). Based on stress-dependent laws of statistical seismology, the stress evolution during this earthquake sequence could be illustrated.

20
Based on the outcomes, we provide additional insight into the FTLS (Gulia and Wiemer, 2019).
Using the b-value of the GR frequency-magnitude law and the p-value of the OU aftershock-decay law, it was shown that a zone of low b-values (indicative of high stress) lay near the eventual epicenters of the M6.5 and M7.3 earthquakes before the start of the Kumamoto sequence, and that after it, the overall trend in stress decreased along the Futagawa-Hinagu FZ.
Detailed analysis suggested aseismic preslips and afterslips along the causative FZ. Stress reduction just prior to the M7.3 mainshock near its epicenter was suggested to be locally induced by preslips.
Based on the review, we noted that the observed b-values were indicative of a green FTLS status before the M7.3 event, in contrast to a red FTLS status shown by Gulia and Wiemer (2019), who did not study spatial variability of b-values to help detect aseismic preslips that occurred locally on the Futagawa FZ. In our view, taking this b-value variability into account, together with using other seismological and geological information such as p-values and the spatial configuration of active faults, is important to improve the FTLS. It would also be of interest for future work to improve FTLS, taking into account the triggering of a red traffic light when signals to suggest ongoing preslips are detected.

Declaration of competing interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The JMA earthquake catalog used in this study was obtained from http://www.data.jma.go.jp/svd/eqev/data/bulletin/index_e.html. We obtained the afterslip planes shown in Fig. 2 from those shown in Fig. 1 of Pollitz et al. (2017). b-values were calculated for events falling in a cylindrical volume with radius r = 5 km, centered at