Changes in Seismicity Pattern Due to the 2016 Kumamoto Earthquakes Identify a Highly Stressed Area on the Hinagu Fault Zone

A highly stressed area where eventual ruptures have often been observed to nucleate is characterized by low b values of earthquake frequency‐size distribution. Crustal deformation due to the occurrence of large earthquakes causes stress perturbation in nearby regions, so an investigation into spatiotemporal b values can play a crucial role in the distribution of postseismic hazards after the 2016 Kumamoto earthquake sequence along the Futagawa‐Hinagu fault zone, which culminated in the magnitude 7.3 mainshock. Together with an analysis of aftershock decay p value that can be used to infer stressing history, a highly stressed area with a characteristic dimension of 10 km at the southern end of the causative faults was found. Our observation is explained by postseismic deformation due to an afterslip on the causative faults and viscoelastic relaxation model. Similar to the Kumamoto mainshock rupture, which started at a low‐b‐value area, the observed highly stressed area shows a high likelihood of future earthquake ruptures.

Slow deformation caused by slow movements of the causative Futagawa-Hinagu FZ and flow in the Earth's interior has been observed in the aftermath of the Kumamoto earthquakes. Data from the Global Navigation Satellite Systems Earth Observation Network System detected postseismic crustal deformation of the Kumamoto earthquakes (Hiyama et al., 2016). To explain this deformation, several researchers considered models of joint afterslip and viscoelastic relaxation at depths driven by the Kumamoto coseismic ruptures (Fuwa & Ohzono, 2018;Moore et al., 2017;Pollitz et al., 2017). Viscoelastic modeling strongly favors a very weak lower crust associated with the Beppu-Shimabara graben and the locus of central Kyushu volcanism (Pollitz et al., 2017). The coseimic ruptures and postseismic deformation, which caused stress perturbation along the Futagawa-Hinagu FZ, play a crucial role in the distribution of seismic hazards (Moore et al., 2017). epicenters of the mainshock and two M6.5-class foreshocks. Since b values are known to be inversely dependent on differential stress (Scholz, 1968(Scholz, , 2015, the observed low b values in the epicenter area was interpreted as an indication of highly stressed area. Seismicity along the Futagawa-Hinagu FZ is decaying with time but remains active. If the stress dependence of the b value holds for the Kumamoto area, b values based on current seismicity could be used to identify regions of increased seismic hazard. Here we examined the spatiotemporal variation in b along the Futagawa-Hinagu FZ. To support b value analysis, we also investigated the p value of the Omori-Utsu (OU) aftershock decay law (Utsu, 1961) because p values can be used to infer stressing history (Dieterich, 1994).

Methods and Data
The GR law is given as log 10 N = A − bM, where N is the cumulative number of earthquakes with magnitude larger than or equal to M, A characterizes seismic activity or earthquake productivity of a region, and the constant b is used to describe the relative occurrence of large and small events (i.e., a high b value indicates a larger proportion of small earthquakes, and vice versa). Spatiotemporal changes in b are known to reflect a state of stress in the Earth's crust (Narteau et al., 2009;Smith, 1981) and to be influenced by asperities and frictional properties (Hirose et al., 2002;Schorlemmer et al., 2004;Tormann et al., 2013;Yabe, 2003) and by interface locking along subduction zones (e.g. Ghosh et al., 2008, Nanjo & Yoshida, 2018, Sobiesiak et al., 2007, 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;Mogi, 1962;Scholz, 1968). To estimate b values consistently over space and time, we employed the EMR (entire-magnitude range) technique (Woessner & 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 (Aki, 1965) in computing the b value to events with a magnitude above M c . Uncertainties in b were computed by bootstrapping (Schorlemmer et al., 2003). The inset of Figure 1b shows a good fit of the GR law to observations in the present cases. 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 a magnitude greater than or equal to a cutoff magnitude; and c, k, and p are constants. p = 1 is generally a good approximation (Omori, 1894). A good fit of the OU law to observations is seen in the present cases (supporting information Figures S1 and S2). The rate-and statedependent friction law (Dieterich, 1994) assumes p = 1 as a standard form. Variability in p is possible: p > 1 and p < 1 for special cases with rapidly and slowly decreasing stress, respectively. Similar to the GR case, we used the maximum-likelihood fit to determine the parameters for this law. Uncertainties in p were computed by bootstrapping. The codes of the OU-fitting method and the EMR method are available together with the seismicity analysis software package ZMAP (Wiemer, 2001).
Our data set is the earthquake catalog maintained by the Japan Meteorological Agency. We used earthquakes having M ≥ 0 since 2000 with depths shallower than 25 km within the study region.  used earthquakes until immediately before the Kumamoto earthquakes for their analysis, while we used them until 11 March 2019.
In addition to the Kumamoto earthquakes, two M5-class earthquakes occurred in the study region: one was the M5.0 event in 2000 near the epicenter of the mainshock (M7.3), and the other was the M5.1 event on 3 January 2019, further to the northwest of the Futagawa FZ. ERC (2019) announced that the latter earthquake was unrelated to the Kumamoto earthquakes. Both events are not considered as the Kumamoto earthquakes so that they were not dealt with in this study.

Results
A map view ( Figure 1a) of b values based on seismicity over the period since 1 January 2017 reveals an area of low b values (b~0.7, light blue), which lies further to the southwest of the M6.5-class foreshocks (red stars). To show this current map, we did not consider the period from immediately after the mainshock to the end of 2016 to remove the effect of the main aftershock sequence that reveals strong temporal variability in b ( Figure S3). To create Figure 1a, we followed  and calculated b values for events falling within a cylindrical volume with a 5-km radius, centered at each node on a 0.02°× 0.02°grid, and plotted a b value at the corresponding node if at least 50 events in the cylinder yielded a good fit to the GR law. The characteristic dimension of this low-b-value area in approximately 0.1-0.2°is 10-20 km. We compared this map with the distribution of M ≥ 5 earthquakes that occurred during the period from the start of the Kumamoto sequence to the end of 2016 (stars, squares, and circles). This low-b-value area is located near the southern end of M ≥ 5 earthquakes.  absolute difference in b is larger than the sum of uncertainties of the b values (inset of Figure 1b). This significant change in b holds true (Δb = −0.12, logP b = −2.09) when compared with the b value based on seismicity since 2018 with (b, a, M c ) = (0.71 ± 0.03, 3.04, 0.5). We computed Δb and logP b for nine regions along the Futagawa-Hinagu FZ (Figures 2a, 2b, and 2d), where Region 6 is the same as that shown in Figures 1a and 1b. Δb ≥ 0 is observed for all regions (Figure 2a), except for Region 6, interpreted as an indication of overall stress release along this FZ with a locally stress-increasing region. Figure 2b shows that Δb values are significant, except for Regions 7 and 9. These regions are located farther south so that the influence of the coseismic ruptures on stress perturbation in the regions is interpreted to be less pronounced. We checked if the extent of the observed surface rapture (Shirahama et al., 2016) corresponds to that of the area with the significant b value change ( Figure 1) and found no clear correspondence.
As a preliminary analysis to modeling the OU law, the cumulative number of earthquakes was calculated as a function of time since 2000 along the Hinagu FZ (inset of Figure S1). Seismicity is still decaying, and the most recent rate of earthquakes has not yet dropped to normal levels. Figure S1 shows a good fit of the OU law to data. We next estimated p values for nine regions along the Futagawa-Hinagu FZ (Figure 2c). For Regions 1-6, p values are 1-1.4, although relatively high p values observed in the eastern regions (Regions 1 and 2) may be possible to relate to the presence of the volcanically active region near Mount Aso. For Regions 7-9, p values are below a typical value (p < 1) if unreliable estimates (open circles) are not taken into consideration. The difference in p between the northern and southern regions of the Hinagu FZ is statistically significant while considering uncertainties in p. The decay of seismic activity is faster in northern regions than in southern ones. The feature is not sensitive to considering large regions for estimation of p values ( Figure S2).

Discussion and Conclusion
The M7.3 mainshock started at the westernmost part of the Futagawa FZ, and its rupture propagated toward the east (Yagi et al., 2016). M5-class aftershocks which occurred along this FZ lay in and around the eventual rupture area of the mainshock (Figure 1). On the other hand, two large foreshocks (M6.5 and 6.4) and M5-class foreshocks and aftershocks along the northern section of the Hinagu FZ indicate that this section has been activating (ERC, 2016).
Our analysis based on seismicity after 1 January 2017 pinpointed an area of low b values along the Futagawa-Hinagu FZ (Figure 1). This is the only area where the b value has significantly decreased, compared to the b value of the pre-Kumamoto seismicity (Figures 1 and 2). This significance was confirmed by an Utsu test and bootstrapping errors. Furthermore, the low-b-value area for seismicity in and after 2017 is located at the southernmost area of the activated northern section of the Hinagu FZ, beyond which no M5-class earthquake occurred (Figure 1a).
Seismicity along the Futagawa-Hinagu FZ was well modeled by the OU law. Large p values (p = 1-1.4) were obtained when we used data along the Futagawa FZ and the northern section of the Hinagu FZ, while they were small (p < 1) in the southern section (Figures 2c and S2). Several researchers demonstrated postseismic crustal deformation of the Kumamoto earthquakes detected by Global Navigation Satellite Systems Earth Observation Network System data (Hiyama et al., 2016) and considered models of joint afterslip and viscoelastic relaxation in the Earth's interior, driven by the Kumamoto coseismic ruptures (Fuwa & Ohzono, 2018;Moore et al., 2017;Pollitz et al., 2017). Pollitz et al. (2017) assumed afterslip planes representing the Futagawa FZ and the northern section of the Hinagu FZ (Figure 1). The south end of afterslip planes representing the latter is near the south end of M5-class aftershocks. They found that a shallow afterslip dominates the postseismic relaxation in the near field, while viscoelastic relaxation of the lower crust and mantle dominates at a greater distance. The Futagawa FZ and the northern section of the Hinagu FZ are included in the afterslip-dominant regime causing rapidly decreasing stress, and the southern section of the Hinagu FZ is in the viscoelastic-flow-dominant regime causing slowly decreasing stress.
Although the overall trend in stress decreases along the Futagawa-Hinagu FZ, the afterslip on the northern section of the Hinagu FZ, in contrast to the lack of an afterslip on the southern section of the same FZ, may not lead to a decrease in local stress at the boundary between the northern and southern sections but to rather a constant or even an increase in local stress there. If we consider a model for seismicity rate resulting from stressing history (Figure 8 of Dieterich, 1994), this explains our result (Figure 2c): High p values (p > 1) in the Futagawa FZ and the northern section of the Hinagu FZ are due to rapidly decreasing stress, low p values (p < 1) in the southern section of the Hinagu FZ are due to slowly decreasing stress, and p~1 near the boundary between the northern and southern sections is due to constant or increasing stress. This is supported by the b value analysis for the main aftershock sequence ( Figure S3) implying that along-fault regions generally exhibited stress decreasing (b value increase), except for a region around the southern end of the northern section of the Hinagu FZ (Region 6) where relatively time-independent b values are shown.
These results explain our observation (Figures 1 and 2) that b values became small near the south end of the northern section after the main aftershock sequence, compared with the b values of the pre-Kumamoto sequence. We further conducted detailed analysis of b values by dividing seismicity from 1 January 2017 to 11 March 2019 into two sub-data sets: one before 1 January 2018 and the other after it. Then, the same mapping procedure as that for Figure 1 was applied to each sub-data set to create the maps in Figures 3a and 3b. Low b values (blue to light blue) broadly distributed before 1 January 2018 (Figure 3b) became narrowly distributed after this date (Figure 3a). The characteristic dimension of unusually low b values in an approximately 0.1°× 0.1°area is ≈10 km. For earthquakes in a cylindrical volume with a radius r = 5 km, shown in Figures 3a and 3b, b values that were large (b > 0.8) before the start of the Kumamoto sequence showed a decrease over time, to values of b~0.6 (Figure 3c), with a large fluctuation due to active aftershocks after the start of the sequence. The same feature is seen for r = 3 km (Figure 3c), interpreting low b values as an indication of a highly stressed area in the south end of the same section. Figure 4 shows a map figure to summarize all of the observations and to spatially visualize our stress model along the Futgawa-Higanu FZ.
A simple check was conducted to assess whether the decrease in b in Figures 1-3 is correlated with the change in faulting type . This was achieved by using the focal mechanism catalog of F-net by National Research Institute for Earth Science and Disaster Resilience (Okada et al., 2004). The predominant type of faulting for the Futagawa-Hinagu FZ was strike-slip ( Figure S4), consistent mechanism of fault motion with the mainshock and M6.5-class foreshocks, although some events involved normal faulting. No correlation between faulting type and time (predominantly, strike-slip faulting at any time) was found in the circle of r = 5 km, shown in Figures 3a and 3b.
Overall consistency between our seismicity analysis and the results of geodesy indicates that changes in seismicity pattern in and around the Futagawa-Hinagu FZ reflect postseismic deformation following the Kumamoto earthquakes.
As observed in our previous study on seismicity prior to the Kumamoto earthquakes , the location of stress concentration (low b value) was regarded as a candidate of the initiation of dynamic ruptures for the M6.5-class foreshocks and mainshock. This was consistent with other observations in previous studies that correlated patches of low b value, with the large earthquakes in California, Tohoku-oki, Sumatra, and others (Nanjo et al., 2012;Nuannin et al., 2005;Tormann et al., 2013Tormann et al., , 2015Wiemer & Wyss, 2002). The combined evidence in those studies indicates that b values can be used to predict a highly stressed area, into which nucleation points of eventual ruptures are included (gray segment in Figure 4). The northern section of the Hinagu FZ relaxed stress due to an afterslip, so a future rupture would propagate along the remaining southern section (outlined line in Figure 4).
Calculation of Coulomb stress imparted by the M7.3 mainshock and M6.5 foreshock to the surrounding major FZs does not provide information to explain aftershocks on and very near the Futagawa FZ and the northern section of the Hinagu FZ that ruptured during these earthquakes (Stein & Toda, 2016), so we think that it cannot be used to interpret the observed changes in b and p based on such aftershocks. However, the Coulomb stress was imparted by the two earthquakes to the southern section of the Hinagu FZ associated with remote aftershocks, bringing this section closer to failure. This supports our implication for the possibility of a future earthquake rupture on the same section.
Based on the migration of foreshocks that occurred during the period between the M6.5 foreshock (the start of the Kumamoto earthquake) and the mainshock, an aseismic slip on the northern section of the Hinagu FZ was known to propagate northward to the nucleation point of the mainshock that eventually occurred on the Futagawa FZ (Kato et al., 2016;Nanjo & Yoshida, 2017). It is worthwhile monitoring seismicity migration representing an aseismic slip propagating toward the southern section of Hinagu FZ as well as spatial and temporal changes in b values along the same FZ.