Stress and pore fluid pressure control of seismicity rate changes following the 2016 Kumamoto earthquake, Japan

We quantitatively examined the influence of pore fluid pressure and coseismic stress changes on the seismicity rate changes that followed the 2016 Kumamoto earthquake, on the basis of two approaches. One is a numerical calculation of the classic stress metric of ∆CFS, and the other is an inversion analysis of pore fluid pressure fields with earthquake focal mechanism data. The former calculation demonstrated that seismicity rate changes were consistent with the expectation from ∆CFS in 65% of the target region, whereas they were not in the remaining 35% of the region. The latter analysis indicates that seismicity rates increased in the regions where pore fluid pressure before the Kumamoto earthquake sequence was remarkably enhanced above hydrostatic, regardless of values of ΔCFS. This suggests that the increase in pore fluid pressure is one of the important physical mechanisms triggering aftershock generation. We obtained evidence that pore fluid pressure increased around the southern part of the main rupture zone after the mainshock, examining temporal changes in types of focal mechanism data. The average increases in pore fluid pressure were estimated to be 17, 20, and 17 MPa at depths of 5, 10, and 15 km, respectively. These large increases in pore fluid pressure cannot be explained under the undrained condition. The spatial derivative of the pore fluid pressure field in the depth direction implies that fluid supply from greater depths may have controlled increases in seismicity rates that followed the large earthquake.


Introduction
Numerous aftershocks often follow large earthquakes in and around the main rupture zone. These phenomena suddenly increase the seismicity rate in some regions and decrease it in other regions after the mainshock. This compartmentalization has been usually explained using the Coulomb failure stress change (∆CFS) (e.g., Reasenberg and Simpson 1992;King et al. 1994;Kilb et al. 2002; Open Access *Correspondence: terakawa@seis.nagoya-u.ac.jp 1 Graduate School of Environmental Studies, Nagoya University, D2-2 (510) Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan Full list of author information is available at the end of the article Toda et al. 2011). The premise for ∆CFS is that the occurrence of earthquakes, in which aftershocks are not an exception, is governed by the Coulomb failure criterion: where τ s is the fault strength, σ n is the normal (positive in compression) stress, P f is the pore fluid pressure, and µ is the friction coefficient of the fault. The Coulomb failure stress (CFS) is defined by the difference between shear stress and fault strength of a receiver fault. The ΔCFS is its change caused by a preceding event: (1) τ s = µ(σ n − P f ), where Δτ and Δσ n are coseismic changes in shear and normal stresses, and ΔP f is that in pore fluid pressure. Theoretically, we expect that an increase in shear stress (the first term of the right side of Eq. (2) and/or a decrease in fault strength (the second term of the right side of Eq. (2) promote triggering aftershocks, whereas a decrease in shear stress and/or an increase in fault strength inhibit it.
The ΔCFS has been the most widely and usefully used for understanding the spatial pattern of aftershocks for four decades. Meanwhile, the applicability of the model has been examined and dissected (e.g., Hardebeck et al. 1998;Mallman and Zoback 2007). One of the essential causes of the inapplicability may be attributed to the fact that the evaluation of ΔCFS mainly depends on coseismic stress changes. In evaluating ΔCFS values, we need to know coseismic changes in stress and pore fluid pressure fields caused by the mainshock. Coseismic stress changes can be calculated from coseismic slip distributions using slip response function (e.g., Okada 1992;Fukahata and Matsuura 2005). In contrast estimating coseismic changes in pore fluid pressure is more difficult. On the assumption of isotropic homogeneous poroelastic material under the undrained condition, coseismic changes in pore fluid pressure are usually evaluated by where B (0 ≤ B ≤ 1) is Skempton's coefficient (e.g., Kilb et al. 2002). With this equation, the ΔCFS in Eq. (2) is rewritten by where µ ′ = µ(1 − B) is apparent friction coefficient of the fault. However, this treatment amounts to replacing the friction coefficient to a lower value, which results in underestimating the effects of decreases in fault strength on ΔCFS. Fluids associated with large earthquakes migrate from overpressurized fluid reservoirs in and around the source region (e.g., Nur and Booker 1972;Sibson 1990;Miller et al. 2004). These phenomena under the drained conditions drastically increase pore fluid pressure in some fault zones and reduce the fault strength, which may trigger numerous aftershocks (e.g., Hubbert and Rubey 1959;Terakawa et al. 2013Terakawa et al. , 2020a. Focal mechanism tomography (FMT) is an inversion technique that estimates 3-D pore fluid pressure fields from information on earthquake focal mechanism data relative to stress pattern or a deviatoric stress tensor normalized by maximum shear stress (Terakawa et al. , 2012Terakawa 2014). Analyses with this technique have revealed the existence of overpressurized fluid reservoirs in the source region of induced, (3) �P f ≡ B�σ n , (4) �CFS ≡ �τ − µ ′ �σ n , triggered, and natural earthquakes (e.g., Terakawa et al. , 2012Terakawa et al. , 2020bTerakawa 2014Terakawa , 2017.
In this study, we target the 2016 Kumamoto earthquake (M w 7.0) in the Kyushu island, Japan (Fig. 1), and quantitatively investigate the influence of coseismic changes in pore fluid pressure as well as coseismic stress changes on seismicity rate changes after the mainshock. This event occurred in the Beppu-Shimabara graben, along which active volcanoes exist. This area is also known for high geothermal and seismic activities, which makes it an appropriate site for this study. In Sect. 2, we first examine the seismicity rate changes for 3 years before and after the Kumamoto earthquake. In Sect. 3, we calculate the coseismic stress changes caused by the Kumamoto earthquake and evaluate the distribution of ΔCFS. In Sect. 4, we estimate the 3-D pore fluid pressure field before the Kumamoto earthquake using the technique of FMT. Finally, in Sect. 5, we examine relationships among the seismicity rate change, the pore fluid pressure field, and ΔCFS. We try to understand the roles of stress and fluids to seismicity rate changes after the large earthquake. Fig. 1 Tectonic setting of the Kyushu island, Japan. The light blue and yellow stars and focal spheres denote the hypocenters and focal mechanisms of the mainshock and largest foreshock of the 2016 Kumamoto earthquake. The red triangles denote active volcanoes. The black lines denote the iso-depth contours of the upper surface of the Philippine Sea plate (Hashimoto et al. 2004

Seismicity rate changes caused by the Kumamoto earthquake sequence
The 2016 Kumamoto earthquake sequence started with the occurrence of an M w 6.1 event, which was eventually the largest foreshock, on April 14, 2016, at the Hinagu fault (Fig. 1). This event was followed by the M w 7.0 mainshock along the Futagawa-Hinagu faults on April 16, 2016. The hypocentral depths of the foreshock and mainshock estimated by the Japan Meteorological Agency (JMA) were 11 and 12 km, respectively. This seismic sequence occurred in the Beppu-Shimabara graben with north-south extension, which runs east and west in the central part of the Kyushu island (Fig. 1). Controlled by the tectonic setting, the focal mechanism solutions for both the largest foreshock and the mainshock were strike-slip types with normal components of north-south tension, and the strikes, dip angles, and rake angles listed in the F-net seismic moment tensor catalog are 212°, 89°, and − 164° for the largest foreshock and 226°, 84°, and − 142° for the mainshock.
We examined the seismicity rate changes that followed the start of the Kumamoto earthquake sequence using the JMA unified hypocenter catalog. Concretely, we targeted the whole Kyushu island as the model region (longitude: 129.5°N-132.0°N, latitude: 31.0°E-34.0°E, depth: 0-20 km) and set 9273 evaluation points at 5-km intervals in the horizontal planes at depths of 5, 10, and 15 km. Then, we counted the numbers of events within 5 km of each evaluation point for the preseismic (April 14, 2013, to April 13, 2016 and postseismic (April 14, 2016, to April 13, 2019 periods, respectively, and evaluated the ratio of the number for the latter period (R a ) to that for the former period (R b ) as the seismicity rate change. When both R b and R a are zero, we set the ratio to be 1, which means no seismicity rate change. When one of them is zero, we replace zero to one, and calculate the ratio. This means that we represent seismicity rate changes as the number of events for the postseismic period with R b = 0 and R a ≠ 0, and the reciprocal number of events for the preseismic period with R b ≠ 0 and R a = 0.
The total numbers of earthquakes in the model region are 34,879 and 183,768 for the preseismic and postseismic periods, respectively. Figure 2 shows the seismicity rate changes. Seismicity rates increased along the Beppu-Shimabara graben, including the main rupture zone. On the western side of the southern part of the main rupture zone (the Hinagu fault), they also remarkably increased. In contrast, they decreased in the northern and southern Fig. 2 Seismicity rate changes that followed the Kumamoto earthquake sequence. Seismicity rate changes at depths of (a) 5 km, (b) 10 km, and (c) 15 km. The seismicity rate change is shown by the logarithm of the ratio of the number of events for the postseismic period (R a ) to that for the preseismic period (R b ). The stars denote the epicenter of the mainshock. The triangles denote active volcanoes. The black dashed lines show prefectural boundaries. The light green rectangles show the main rupture zone of the Kumamoto earthquake by Asano and Iwata (2016). The light green dashed lines show the study region , where seismicity rate changes are investigated in this study. Grey lines show active faults sides of the northern part of the main rupture zone (the Futagawa fault). The spatial characteristics in seismicity rate changes are independent of depths. This is consistent with the fact that the Kumamoto earthquake was primarily strike-slip faulting.

Coulomb failure stress changes caused by the Kumamoto earthquake sequence
To evaluate the Coulomb failure stress change (ΔCFS) caused by the Kumamoto earthquake sequence, we need to know the background stress pattern before the seismic sequence as well as coseismic stress changes. The information of the background stress pattern is used for determining the orientation of receiver faults and the positive direction of shear stress changes on them. In Sect. 3.1, we estimated the background stress field in the same model region in Sect. 2 by applying a dataset of earthquake focal mechanisms to the CMT data inversion method (Terakawa and Matsuura 2008). In Sect. 3.2, we calculated coseismic stress changes caused by the 2016 Kumamoto earthquake (the mainshock) and the largest foreshock to evaluate ΔCFS.

Background stress pattern before the Kumamoto earthquake sequence
The CMT data inversion method (Terakawa and Matsuura 2008) is an inversion method to statistically derive each of the six components of a deviatoric stress tensor, based on the Bayesian statistical inference and Akaike's Bayesian Information criterion (ABIC) (Akaike 1977(Akaike , 1980. In this method, centroid moment tensor (CMT) data are related to the stress pattern in and around the source, using the definition of the moment tensor of earthquakes (e.g., Matsuura et al. 2019). We put two catalogs together not to duplicate the same events, and made a dataset for the analysis. One dataset consists of 2494 earthquake focal mechanism data (M ≥ 1) during the period from January 24, 1996, to April 7, 2016, obtained by Matsumoto et al. (2018). The other one consists of 3,115 focal mechanism data (M ≥ 3) in the Kyushu district (longitude: 127°N-135°N, latitude: 27°E-37°E) during the period from January 1, 1997, to April 13, 2016, listed in the F-net seismic moment tensor catalog. The number of data in the original dataset is 5306 (Fig. 3). We used the CMT data converted from these focal mechanism data using the well-known relationship between moment magnitude and seismic moment (Hanks and Kanamori 1979).
To discretize the problem, we distributed 8029 (31 × 37 × 7) tricubic B-splines (basis functions) with equally spaced 10-and 5-km local supports (grid intervals) in horizontal and vertical directions, respectively, to represent the stress field in the model region in Sect. 2. The original dataset has 2516 events within the model region, and they comprise the final dataset. We determined the best estimates of the model parameters (expansion coefficients of B-splines) from the dataset using the inversion formula of Yabuki and Matsuura (1992). Note that only the relative values of the six stress components were meaningful; the absolute values were not.
The background stress pattern in the Kyushu island is strike-slip to normal faulting types with north-south tension, which is consistent with the Beppu-Shimabara graben structure (e.g., Terakawa and Matsuura 2010;Matsumoto et al. 2015Matsumoto et al. , 2018 (Fig. 4). The consistency of focal mechanisms of the Kumamoto earthquake and the largest foreshock with the background stress pattern also indicates that these events played a role that released tectonic stress in the wide range of the Kyushu island.
We also evaluated the uncertainty for the stress pattern using the average inner tensor product (ITP) between the best stress tensor and a possible stress tensor (Terakawa 2017). We define a metric for uncertainties as Fig. 3 Earthquake focal mechanism data for the preseismic period of the Kumamoto earthquake sequence. Focal mechanism data are represented with lower hemisphere projections of focal spheres. The color scales denote hypocentral depths. The star, the black dashed lines, and the grey lines are the same as in Fig. 2. The black and red dashed lines show the model region for the background stress field and pore fluid pressure field and the study region for detailed investigation, respectively Nakagomi et al. Earth, Planets and Space (2021) where smaller values of E mean that the stress pattern is more reliable. Uncertainties in the source region of the Kumamoto earthquake sequence were smaller than 0.3 (Fig. 4). The values of 0.1, 0.4, and 2.0 roughly correspond to the rotation angles of the maximum compressive principal stress axis of 13°, 30°, and 90°, respectively.

Distribution of ΔCFS caused by the Kumamoto earthquake sequence
To evaluate ΔCFS, we calculated a coseismic stress change field using the analytical slip response function by Fukahata and Matsuura (2005) with the source rupture model by Asano and Iwata (2016). In this calculation, we considered both the effects of the Kumamoto earthquake (the mainshock) and the largest foreshock (Fig. 5). The rupture model for the mainshock consisted of the northern and southern segments with lengths of 28 and 14 km, respectively. The fault length for the largest foreshock was set to be 14 km. Tables 1 and 2 summarize the elastic parameters for the medium and the detailed information on the two rupture models for the mainshock and the largest foreshock, respectively. To estimate the coseismic stress change in as much detail as possible, we obtained smooth distributions of coseismic slip on the fault segments (Additional file 1: Figure S1) by applying an inversion formula based on the Bayesian statistical inference with ABIC to the original model of Asano and Iwata (2016).
Here, we focused on coseismic stress change fields in the region (longitude: 130.2°N-131.6°N, latitude: 32.0°E-33.4°E, depth: 5-15 km) with small uncertainties of the background stress pattern (Fig. 4), because we cannot evaluate ΔCFS without reliable background stress patterns. Hereafter, we call the smaller region for detailed analyses as the study region. We did not evaluate the coseismic stress change field within 1 km from the source because of technical difficulty in handling stress near the source.
In the northern and southern sides of the main rupture zone, the pattern of coseismic stress change was primarily characterized by reverse and/or strike-slip faulting with north-south compression. In contrast, it was characterized by normal faulting and/or strike-slip faulting with east-west tension in the east and west sides of the source region. The maximum and minimum compressive principal stress axes of the coseismic stress change field were opposite to those of the background tectonic stress field, indicating that the Kumamoto earthquake sequence released some of the accumulated tectonic stress in a wide region in and around the source. The deviatoric stress magnitude (the maximum shear stress) of coseismic stress changes decreased with distance from the source faults. The maximum values of the coseismic stress changes were 4.5, 7.3, and 4.7 MPa at depths of 5, 10, and 15 km, respectively.
Next, we calculated ΔCFS in Eq. (2), taking maximum shear stress planes of the background tectonic stress field (Fig. 4) as receiver faults. We took directions of the resolved shear traction of the background stress field as positive in shear stress, indicating that Δτ increases/ decreases when shear stress is accumulated/released by the two preceding events. In this calculation, we evaluated the values of ΔCFS on both maximum shear stress planes at each evaluation point and plotted the larger one on the map. Here, we ignored the effects of changes in pore fluid pressure. Figure 6 shows the distribution of ΔCFS due to the largest foreshock and the mainshock of the Kumamoto earthquake. In regions of the eastern and western sides of the main rupture zone, values of ΔCFS were positive. The values of ΔCFS were especially positive along the Beppu-Shimabara graben except for the main rupture zone. Conversely, in the northern and southern sides of the main rupture zone, they were negative. We expect that seismicity rates in regions with positive and negative   values of ΔCFS increased and decreased for the postseismic period, respectively. The characteristic distribution of ΔCFS is roughly consistent with the seismicity rate changes (Fig. 2). Further investigations are shown in Sect. 5.
Pore fluid pressure field before the Kumamoto earthquake sequence FMT (Terakawa et al. 2010, 2012Terakawa 2014) is an inversion technique to estimate 3-D pore fluid pressure fields from earthquake focal mechanism data as a continuous function defined in the model region, based on the same Bayesian inversion scheme with ABIC (Akaike 1977(Akaike , 1980Yabuki and Matsuura 1992) as that used in the CMT data inversion method. The primary assumption for FMT is the Wallace-Bott hypothesis that seismic slip occurs in the direction of a resolved shear traction acting on a preexisting fault (Wallace 1951;Bott 1959) governed by Coulomb's failure criterion. We, therefore, can expect diverse earthquake focal mechanisms even in a uniform stress field. Two end-member models can explain the variation in the focal mechanism (e.g., Zoback 1992; Rivera and Kanamori 2002). One model attributes the variation to differences in frictional coefficients of rocks, assuming that the pore fluid pressure field is uniform everywhere. The other model assumes variation in a pore fluid pressure field with a constant friction coefficient (e.g., Byerlee 1978;Zoback and Townend 2001). In FMT, we estimate pore fluid pressure fields based on the latter model, examining a type of focal mechanism relative to the stress pattern in the 3-D Mohr diagram with a constant friction coefficient. The detailed procedure is summarized in Terakawa et al. (2012). We targeted the same model region as that for estimating the background stress field in Sect. 3.1 and discretized the pore fluid pressure field by a superposition of the same tricubic B-splines. Given the background tectonic stress pattern (Fig. 4), a rock density of 2700 kg/ m 3 (Table 1), and a fixed friction coefficient (0.6), we determined the 3-D pore fluid pressure field before the Kumamoto earthquake sequence from the final dataset in Sect. 3.1 (Fig. 3). To select the true fault plane from two nodal planes, we evaluated the misfit angle between the observed slip vector and that expected from the stress pattern, which is attributed to the errors for focal mechanism data and the background stress pattern. Considering these uncertainties, we used only focal mechanism data with a misfit angle of < 25° for the FMT analysis. The number of data for the FMT analysis was 2351. In cases where the misfit angles for both the nodal planes were less than the threshold, we determined the one at which the seismic slip could be triggered under lower pore fluid pressure. The average misfit angle of the focal mechanism data to the stress pattern was 9.9°.
To compare the degree of overpressure regardless of depths, we define the overpressure coefficient C by The stars, the light green rectangles, the triangles, the black dahsed lines, and the grey lines are the same as in Fig. 2 where P H and P L are the hydrostatic and lithostatic pressures, respectively. Because ambiguity in the background stress pattern causes the modeling error (Terakawa 2017;Terakawa et al. 2020b), we focused on the results of pore fluid pressure fields in the study region. Figure 7 shows the distribution of the overpressure coefficient at depths of 5, 10, and 15 km. Figure 8 shows the estimation errors (the 68% confidence regions) for the estimates of overpressure coefficients. We found that pore fluid pressure was highly enhanced above hydrostatic near the hypocenters of the largest foreshock and the mainshock.
Overpressurized fluids were distributed in the wider regions at shallower depths. This tendency is consistent with the results of seismic tomography in the whole Kyushu island, where higher Poisson's ratios were distributed in the wider regions at a depth of 4 km than that at a depth of 13 km (Zhao et al. 2018). Figures 1 and 7 show that the Hinagu fault, the northermost part of which corresponds to the southern segment of the main rupture zone, plays the role of a structural boundary that separates overpressurized fluid reservoirs. Regions with overpressurized fluids were located on the eastern side of the northern part of the Hinagu fault at depths of 5 and 15 km and on the western side of the central part of Hinagu fault at depths of 5-15 km. The former corresponds to a localized conductive region detected by magnetotelluric (MT) and telluric data survey, and the latter may do a relatively conductive region beneath which a conductive zone exists at depths of ~ 20 km (Aizawa et al. 2017). This comparison may not necessarily indicate a clear coincidence between pore fluid pressure fields estimated in this study and the resistivity structure obtained by the MT survey. In general, we cannot estimate pore fluid pressure in highly conductive regions, where seismicity is inactive (Fig. 8). Therefore, overpressurized fluid reservoirs estimated from focal mechanism data may correspond to a relatively conductive zone detected by the MT survey. Aizawa et al. (2017) also detected highly conductive zones in the Beppu-Shimabara graben, where active volcanoes exist. In this conductive zone, seismicity may be lacking, because stable sliding is more dominant than dynamic rupture of earthquakes under extremely high pore fluid pressure. We cannot evaluate pore fluid pressure in the corresponding regions because of a lack of data (Fig. 8).
We also investigated the distribution of pore fluid pressure on the mainshock faults (Additional file 1: Figure S2). This figure indicates that pore fluid pressure was enhanced above hydrostatic below the hypocenters of the mainshock and the largest foreshock. Comparing the results with the coseismic slip distribution Fig. 7 3-D pore fluid pressure field before the Kumamoto earthquake sequence. Pore fluid pressure fields at depths of (a) 5 km, (b) 10 km, and (c) 15 km. The excess pore fluid pressure above hydrostatic is shown with overpressure coefficient C (see text). The red lines show contour lines for C = 0.13. The stars, the black dashed lines, and triangles are the same as in Fig. 2. The white lines denote active faults. The grey rectangles show the main rupture zone of the mainshock (Additional file 1: Figure S1), we can see that coseismic slip became larger in regions with low pore fluid pressure or high strength. These results are consistent with previous studies of seismic tomography in and around the source region of the Kumamoto earthquake, where the mainshock occurred in a high velocity and low Poisson's ratio and V p /V s (e.g., Shito et al. 2017;Zhao et al. 2018). This suggests that large shear stress accumulated in the regions with high strength was released by the Kumamoto earthquake sequence.

Physical quantities controlling seismicity rate change: stress and pore fluid pressure
We compared the seismicity rate change (Fig. 2) with the distribution of ΔCFS at each evaluation point (Fig. 6). Figure 9 shows the relationship between ΔCFS and the seismicity rate change. Seismicity rates tended to increase and decrease in the regions with ΔCFS > 0 and ΔCFS < 0, respectively. Seismicity rates also tend to increase more as values of ΔCFS are larger. These tendencies did not depend on depth.
We divided evaluation points in the study region into four groups: (1) the regions where seismicity rates increased, as expected from positive ΔCFS values (Group 1), (2) those where seismicity rates increased despite negative ΔCFS values (Group 2), (3) those where seismicity rates decreased, as expected from negative ΔCFS (Group 3) values, and (4) those where seismicity rates decreased despite positive ΔCFS values (Group 4). Additional file 1: Figure S3 shows the distribution of evaluation points (Groups 1-4) with remarkable seismicity rate changes on ΔCFS. The original number of the evaluation points was 2463, but we did not calculate coseismic stress changes at 16 points because of difficulty in handling coseismic stress changes near the source region. The proportions of the four groups to the targeted points (2447) were 50%, 28%, 15%, and 7%, respectively. The results indicate that seismicity rate changes at 65% of the target region were consistent with the coseismic stress changes caused by the largest foreshock and the mainshock of the Kumamoto earthquake.
However, in 35% of the target region, seismicity rate changes were inconsistent with the results of ΔCFS. Many events of Group 4 primarily occurred in regions where seismicity rate changes were not remarkable and/or the effects of ΔCFS were small relative to the stress drop (Fig. 9), indicating that these results were affected by uncertainties in evaluating seismicity rate changes and ΔCFS and/ or some local conditions. In contrast, the occurrence of events of Group 2 indicates that seismicity rates remarkably increased in regions where seismicity should have been inhibited by the coseismic stress change ( Fig. 9 and S3). This suggests that some physical mechanisms other than coseismic stress changes control aftershock generation. The numbers of events that occurred within 3 and 5 km from the evaluation points of Group 2 for the postseismic  Fig. 7. The red cycles show data used in the FMT analysis, whose hypocenters are within 2.5 km from each depth. The white stars, white lines, triangles, the black dashed lines, and grey rectangles are the same as in Fig. 7 period were 20,253 and 57,915, which were 14% and 39% of all the 149,250 events for the period, respectively.
Next, we examined the effects of both coseismic stress changes and pore fluid pressure before the Kumamoto earthquake sequence on seismicity rate changes (Fig. 10). Besides the positive correlation of seismicity rate changes to ΔCFS, this figure shows that seismicity rates increased in regions where pore fluid pressure before the Kumamoto earthquake sequence was high (with C > 0.13), regardless of values of ΔCFS (Figs. 7, 10). Additional file 1: Figure S4 shows the distribution of evaluation points with remarkable seismicity rate increases (Groups 1 and 2) on the excess pore fluid pressure field. This figure shows positive correlation between seismicity rate increases and pore fluid pressure. These results indicate that the pore fluid pressure before the large earthquakes is also one of the essential factors that control the seismicity rate change.
If pore fluid pressure in a fault zone was enhanced by large earthquakes, this decreases effective normal stress, which results in fault strength reduction. Therefore, this process can trigger events (e.g., Hubbert and Rubey 1959). On the basis of the idea that negative effects from coseismic stress changes (ΔCFS) were compensated by increases in pore fluid pressure, we estimated the necessity minimum of increases in pore fluid pressure P min f that explains unexpected increases in seismicity rates (at the evaluation points of Group 2) by where normal stress is positive for compression. Additional file 1: Figure S5 shows the distribution of P min at evaluation points of Group 2. This figure shows that non-negligible increases in pore fluid pressure played an important role in triggering seismicity near the main rupture zone where large shear stress was released by the mainshock and the large foreshock (Fig. 5). The maximum values of P min f were estimated to be 5.7, 9.2, and 5.5 MPa near the main rupture zone at depths of 5, 10, and 15 km, respectively. These values of P min f correspond to overpressure coefficients ΔC of 0.069, 0.055, and 0.022, respectively.
Dynamic triggering can be another physical mechanism that can explain seismicity rate changes in Group 2 (e.g., Hill et al. 1993;Hill and Prejean 2015;Kilb et al. 2002;Miyazawa 2016;Uchide et al. 2016). The direct and indirect effects of the process are related to changes in stress and pore fluid pressure, respectively. It would be interesting to consider the effects on seismicity rate changes in the future.

Discussion
The type of focal mechanism data relative to the stress pattern is a measure for pore fluid pressure in the source region of events (e.g., . Given the background stress pattern before the Kumamoto earthquake (Fig. 4), for simplicity, we tried to quantitatively evaluate temporal changes in pore fluid pressures from focal mechanism data. For the analysis, we need to bear in mind possibilities that the background stress pattern was altered by the large earthquakes (e.g., Yoshida et al. 2016;Mitsuoka et al. 2020). Therefore, in this investigation, we targeted the outside regions of the main rupture zone where deviatoric stress magnitudes of coseismic stress change fields were less than 1 MPa, which are small enough not to alter the background stress pattern (e.g., Terakawa and Hauksson 2018).
The dataset of focal mechanisms for the postseismic period consists of 3999 data during the period from April 16, 2016, to September 9, 2017 (Matsumoto et al. 2015;Mitsuoka et al. 2020;Shito et al. 2020) and 784 data listed in the F-net seismic moment tensor catalog during the period from April 14, 2016, to April 30, 2019. We combined these data not to duplicate the same events in the same way as making the original dataset for the analysis in Sect. 3. The number of data in the final dataset was 4635 (Additional file 1: Figure S6).
We targeted 11 regions (five, four, and two evaluation points at depths of 5, 10, and 15 km, respectively) of Group 2 with more than seven focal mechanism data for both the preseismic and postseismic periods (Fig. 11a) and compared pore fluid pressure triggering events for the preseismic period with those for the postseismic period (Fig. 11b). In the figure, we evaluated the average values of overpressure coefficients C at the target regions with standard errors (68%) for both the periods. In all the target regions, the level of pore fluid pressure (C) triggering seismicity increased after the Kumamoto earthquake. We found significant temporal increases in pore fluid pressure in eight regions (four, three, and one regions at depths of 5, 10, and 15 km, respectively). The average increases in overpressure coefficient C (and pore fluid pressure) in the eight regions were 0.04-0.47 (3.1-39 MPa), 0.03-0.18 (5.8-31 MPa), and 0.07 (17 MPa) at depths of 5, 10, and 15 km, respectively. These values overwhelm those of the necessity minimum of increases in pore fluid pressure that were 0.16-1.09 MPa, 0.07-0.58 MPa, and 1.36 MPa at depths of 5, 10, and 15 km, respectively (Figs. 11a; Additional file 1: S5). These estimates support that the aftershocks in the regions were triggered by decreases in fault strength because of increases in pore fluid pressure caused by the Kumamoto earthquake. We could not target regions of Group 2 near the main rupture zone because of the condition of coseismic stress magnitude for the analysis. However, the eight estimates of temporal increases in pore fluid pressure are at least comparable to typical stress drop, and tend to be larger than the maximum values of P min f near the source (Additional file 1: Figure S5). This implies that increases in pore fluid pressure can be significant to enhance seismicity even near the main rupture zone.
We also examined temporal changes in pore fluid pressures in the regions of Group 1, where seismicity rates increased, as expected from positive ΔCFS, in the same way as above (Additional file 1: Figure S7). The number of the target regions was 62 (17, 32, and 13 evaluation points at depths of 5, 10, and 15 km, respectively). In 39 regions (11, 20, and 8 evaluation points at depths of 5, 10, and 15 km, respectively), we recognized significant temporal increases in pore fluid pressure. The ratio of the number of regions with significant increases in pore fluid pressure is 62.9%, which is smaller than that of Group 2 (72.7%). However, the amounts of increases in pore fluid pressure are the same as those in the target regions of Group 2. These results suggest that a part of aftershocks in Group 1 may have been controlled not only by coseismic stress changes but also by increases in pore fluid pressure. Additional file 1: Figure S4 indicates that seismicity rates significantly increased especially around the northern part of the Hinagu fault, which may have been controlled by overpressurized fluids. However, it would be difficult to separate the effects of pore fluid pressure and those of coseismic stress changes on seismicity rate changes.
One plausible mechanism for temporal increases in pore fluid pressure is enhancement by coseismic stress changes under the undrained condition. In this case, increases in pore fluid pressure are product of Skempton's coefficient B (0 ≤ B ≤ 1) and changes in normal stress ∆σ n (Eq. 3). Figure 12 shows the upper limit of increases in pore fluid pressure with B = 1 under the undrained condition. This figure indicates that these values are too small to explain the observed temporal changes in pore fluid pressure with focal mechanism data (Fig. 11).
Another plausible mechanism is fault-valve behavior (e.g., Sibson 1990), or fluid migration from overpressurized fluid reservoirs in the deep. We investigated the distribution of the partial derivative of the excess pore fluid pressure field in the depth direction (Fig. 13). Positive values of this quantity mean that fluids can migrate upwards. The eight regions of Group 2 with significant temporal increases in pore fluid pressure (Fig. 11) are located in regions with positive values of the partial derivatives. More interestingly, these regions are concentrated around the overpressurized fluid reservoirs where regions with positive partial derivatives continue at depths of 5-15 km (Figs. 7,13). This implies that fluid supply from greater depths may play an important role in triggering aftershocks in the regions. Origins of overpressureized fluids in the Kyushu district are of great interest. Beneath the Kyushu island the Philippine Sea plate is subducting from the east at the Nankai trough and the Ryukyu trench. The Beppu-Shimabara graben in the central Kyushu is also located between the Median Tectonic Fig. 11 Temporal changes in pore fluid pressure that followed the Kumamoto earthquake (Group 2). a The 11 target regions. b Temporal changes in overpressure coefficients. The circles, diamonds, and triangles in a and b denote depths of the target regions at 5, 10, and 15 km, respectively. The horizontal and vertical coordinates of symbols in b indicate averages of overpressure coefficients of pore fluid pressure that trigger events for the preseismic and postseismic period, respectively. The bars show standard errors (68%). The color scales of each symbol in (a) show the difference of the average of C for the preseismic period from that for the postseismic period. White color means that we did not find any significant temporal increases in C there. The grey rectangles are the same as in Fig. 8 Line and the Okinawa trough (Fig. 1). It will be interesting how this complex tectonic setting is related to the existence of overpressurized fluids in the Kyushu island.
In Sect. 3.2, we calculated ΔCFS on the maximum shear plane, taking directions of the resolved shear traction of the background stress field as positive in shear stress changes. In this calculation, we do not consider Fig. 12 Theoretical estimates of increases in pore fluid pressure due to the Kumamoto earthquake at depths of (a) 5 km, (b) 10 km, and (c) 15 km. The color scales of symbols show overpressure coefficients of upper limits of theoretical increases in pore fluid pressure under the undrained condition (with B = 1). The red circles denote the evaluation points, where we observed significant increases in C after the mainshock. The grey rectangles are the same as in Fig. 8   Fig. 13 Partial derivatives of pore fluid pressure field in the depth direction. The partial derivatives at depths of a 5 km, b 10 km, and c 15 km. Values are truncated outside the region from − 1 to 1 MPa/km. The squares and circles denote evaluation points of Groups 1 and 2 where we observed temporal increases in pore fluid pressure triggering events. The stars, grey lines, and rectangles are the same as in Fig. 2 that the Kumamoto earthquake and the largest foreshock reversed the direction of resolved shear traction on the background stress field. Thus, we implicitly assume that the magnitudes of background stress field are substantially higher than those of coseismic stress changes. Therefore, it may be difficult to examine influences of stress and pore fluid pressure on seismicity rate changes in the main rupture zone where large coseismic stress changes would have altered the background stress pattern (e.g., Yoshida et al. 2016;Mitsuoka et al. 2020). However, our investigation of seismicity rate changes is rational outside of the source region (e.g., Terakawa and Haukkson 2018).
Finally, earthquake generation is controlled by the absolute stress field, which is the sum of background stress fields and coseismic stress change fields. The magnitude of the background stress field is usually unknown, which makes it difficult to examine the effect on seismicity rate changes (e.g., Saito et al. 2018;Matsuura et al. 2019). Recently, Noda et al. (2020) demonstrated the effects of absolute stress fields on aftershock generation in and around the source region of the Kumamoto earthquake, focusing on changes in shear strain energy due to large earthquakes. Terakawa et al. (2020a) proposed an energetics-based stress metric (ΔEFS) that rationally evaluates earthquake generation by considering the effects of changes in shear strain energy, volumetric strain energy, and pore fluid pressure. The stress metric can generalize ΔCFS, which is defined on a specific fault plane, into a failure stress change described with scalar metric defined in three-dimensional stress space. The essential difference of ΔEFS from ΔCFS is that the former depends on background stress fields as well as coseismic stress changes. These recent studies show that considering background stress fields is essentially important in evaluating seismicity rate changes.
Furthermore, Terakawa et al. (2020a) quantitatively showed that increases in pore fluid pressure triggered 20% of aftershocks for 1 year after the 1992 Landers earthquake. In Sect. 5, we showed that 14-39% of aftershocks of the Kumamoto earthquakes were triggered by some factors other than coseismic stress changes. To understand the roles of pore fluid pressure more quantitatively, we need to take the absolute stress field into account in evaluating seismicity rate changes. This further investigation will provide new insight on earthquake generation.

Conclusions
We investigated the influences of stress and pore fluid pressure on seismicity rate changes that followed the 2016 Kumamoto earthquake sequence, based on the numerical calculation of the coseismic stress change by the mainshock and the largest foreshock, and the inversion analysis of pore fluid pressure. We quantitatively demonstrated that pore fluid pressures before the Kumamoto earthquake, as well as coseismic stress changes, controlled aftershock generation. Seismicity rate changes in two-thirds of the whole study region were consistent with the results of ΔCFS. This indicates that the coseismic stress change can primarily control aftershock activity. Seismicity rates increased in regions where pore fluid pressure before the mainshock was high, regardless of the values of ΔCFS. On the basis of temporal changes in types of focal mechanisms relative to the background stress pattern, we obtained the evidence that pore fluid pressure nonnegligibly increased after the Kumamoto earthquake around the southern part of the main rupture zone. Besides coseismic stress changes, pore fluid pressure before the large earthquake is one of the essential factors that control the seismicity rate change.
Additional file 1: Figure S1 Coseismic slip distribution of the mainshock and the largest foreshock of the Kumamoto earthquake (Asano and Iwata, 2016). (a) The northern segment of the mainshock, (b) the southern segment of the mainshock, and (c) the largest foreshock's segment. The stars in (b) and (c) denote the hypocenters of the mainshock and the foreshock. Figure S2 Excess pore fluid pressure on (a) the northern segment, and (b) the southern segment of the mainshock. (c) and (d) show the estimation errors (the 68 % confidence region). The black circles in (c) and (d) denote data used in the FMT analysis, whose hypocenters are within 10 km from the planes. The star is the hypocenter of the mainshock. Figure S3 Distribution of seismicity rate changes at evaluation points on ΔCFS. Seismicity rate changes at depths of (a) 5 km, (b) 10 km, and (c) 15 km. The blue and red pluses denote evaluation points of Groups 1 and 2 where seismicity rates for the postseismic period are more than ten times larger than those for the preseismic period. The red and blue minuses denote evaluation points of Groups 3 and 4 where seismicity rates for the postseismic period are more than ten times smaller than those for the preseismic period. The color scales in the background show values of ΔCFS (which are truncated outside the region from −1.0 to 1.0). The stars and the grey lines are the same as in Figure 2. The black rectangles show the main rupture zone. Figure S4 Distribution of seismicity rate increases on pore fluid pressure fields. Seismicity rate increases at depths of (a) 5 km, (b) 10 km, and (c) 15 km. The pluses and circles denote evaluation points of Groups 1 and 2 where seismicity rates for the postseismic period are more than ten times larger than those for the preseismic period. The stars, the red lines, the white lines, and the grey rectangles are the same as those in Figure 7. Values indicate overpressure coefficients ΔC for P min f . Figure S6 Earthquake focal mechanism data for the postseismic period. Representation of focal mechanism data, the star, the black and red dashed lines, and grey lines are the same as in Figure 3. Figure S7 Temporal changes in pore fluid pressure that followed the Kumamoto earthquake (Group 1). (a) The 62 target regions. (b) Temporal changes in overpressure coefficients. The representation of these figures are the same as in Figure 11.