Stationarity of aftershock activities of the 2016 Central Tottori Prefecture earthquake revealed by dense seismic observation

To clarify the relationship between earthquake occurrence and fluid, we analyzed data from a dense aftershock-observation network with 69 high-gain short-period seismographs installed immediately after the mainshock occurrence (October 21) in the aftershock area of the 2016 Central Tottori Prefecture earthquake. We determined the hypocenters and focal mechanisms of the aftershocks very precisely in the period from October 22 to December 15. We then investigated the temporal changes in the spatial distributions of hypocenters and T-axis azimuths of focal mechanisms. The distributions of aftershock hypocenters and T-axis azimuths are basically temporally stable, except those in limited portions in the shallow layer near the western edge of the aftershock area, where rapid decrease of aftershocks with T-axis azimuths of WSW to west was observed. If fluid rises from the lower crust due to fault rupture, the locations of aftershocks and focal mechanisms may change over time, especially in the deepest part of the aftershock region. However, the temporal change in these parameters was not apparent at depth. These observations suggest that the aftershock activity of the Central Tottori Prefecture earthquake was controlled mainly by stress concentration rather than strength reduction due to high fluid pressure.


Introduction
An earthquake occurs when shear stress acting on a fault approaches the strength of the fault. Therefore, to forecast earthquakes, it is important to clarify the spatiotemporal changes in the stress field near the fault as well as its strength. When the shear stress acting on a fault approaches the fault strength, the shear stress either increases or the fault strength decreases. In extreme cases, an earthquake may occur due to a decrease in the fault strength even if the shear stress is constant.
The Hubbert-Rubey law is a famous hypothesis concerning the mechanism by which fault strength is reduced. According to this hypothesis, effective normal stress decreases as pore-fluid pressure increases, causing the fault to become slippery (Hubbert and Rubey 1959). Evaluation of fluid-injection experiments conducted in the Rangely oilfield (Colorado, USA) has empirically and quantitatively confirmed the relationship between fluid pressure and earthquake occurrence (Raleigh et al. 1976). In fluid-injection experiments in Basel, Switzerland, the occurrence of induced earthquakes was explained by using estimated crustal stresses together with fluid Open Access *Correspondence: iio@rcep.dpri.kyoto-u.ac.jp 1 Research Center for Earthquake Prediction, Disaster Prevention Research Institute, Kyoto University, Gokasho, Uji, Japan Full list of author information is available at the end of the article pressure (e.g., Mukuhira et al. 2017). Spatiotemporal changes in fluid pressure in the source region of the induced earthquakes were also estimated only from the focal mechanisms of the induced earthquake (Terakawa et al. 2012). The relationship between fluid and natural earthquake swarm activities also has been investigated (e.g., Yukutake et al. 2011). Yoshida and Hasegawa (2018) found deep to shallow hypocenter migration and temporal changes in stress drop and fault strength to be related to the swarm earthquakes that started in the Tohoku inland region seven days after the 2011 Tohoku-oki earthquake. Furthermore, a numerical simulation using a rate-and state-dependent friction law and considering pore-fluid pressure change allowed  to successfully explain the observed data as being due to the decrease and subsequent recovery of fault strength that accompanied the rise and diffusion of high-pressure fluid from the lower crust.
These seismic activities consist of relatively small earthquakes; however, the relationship with fluids has also been discussed for the aftershocks of large earthquakes (Miller et al. 2004;Terakawa et al. 2010). Miller et al. (2004) showed that the aftershock activity of M5.7 and M6.0 can be explained by time-dependent changes in fluid pressure due to the diffusion of CO 2 fluid rising from depths under high pressure at depth.
In terms of the movement of fluids associated with fault movement, the relationship between a quartz vein and/ or gold deposit in a mine and the movement of a highangle reverse fault has been reported in detail (e.g., Sibson 1992;Cox 1995). In addition, there has been a report of a dramatic increase in groundwater discharge near the fault after a major earthquake (Sibson 1981). The 1930 Kern County earthquake (California, USA) resulted in an increased groundwater discharge a few hours after the earthquake; the latter remained high for at least 2 months (Sibson 1992). Sibson (1992) also studied the role of highpressure water in the occurrence of large earthquakes as well as aftershocks after a major earthquake.
As discussed above, various past studies have contributed to knowledge of the relationship between fluid activity and the occurrence of earthquakes. The timescale of these phenomena ranges from several days to several weeks, and it is thought that the movement and diffusion of fluid in the seismogenic region control it.
In recent years, new systems have been developed that allow the easy installation of a dense temporary observation network to record seismograms (Miura et al. 2010;Hayashida et al. 2019). Therefore it has become possible to build such a network immediately after the occurrence of a large earthquake. The day after the 2016 Central Tottori Prefecture earthquake (M JMA 6.6) occurred, a dense network installed, using the "0.1 Manten" and "Manten" systems (Miura et al. 2010;Iio et al. 2018;Hayashida et al. 2019), succeeded in capturing aftershock activity. We investigated the temporal change in aftershock activity, especially within several weeks immediately after the occurrence of the mainshock, to clarify the relationship between fluid and the occurrence of an earthquake.

Data
We installed a dense aftershock-observation network with 69 high-gain short-period seismographs in the aftershock region of the 2016 Central Tottori Prefecture earthquake that occurred in the San'in seismic zone (Kawanishi et al. 2009). The purpose was to clarify the absolute value of shear stress around the fault and the fault strength, in addition to any temporal change in aftershock activities.
The distribution of the aftershock-observation stations is shown in Fig. 1. From October 22 to November 8 2016, 10 "Manten" three-component seismographs, 10 Geospace GSX three-component seismographs, and 49 "0.1 Manten" single-component seismographs were installed in and around the aftershock area. To explain the detail of the installation, 20 three-component seismographs and 25 single-component seismographs were installed during October 22-24, and the remaining 24 singlecomponent seismographs were installed during November 7-8. However, data were not recorded at many of the first 25 stations, as there was initially a problem with writing the data to the SD card of the single-component seismograph, After October 29, all the stations initially installed were operating normally. As a result, 26 temporary observation stations were operating until October 28; this increased to 45 up to November 8 and 69 from November 9. The "0.1 Manten" single-component and GSX three-component seismographs were operated until December 12, and March 2017, respectively. Only the "Manten" three-component seismographs are currently still operating. Three-component seismographs are basically fixed with adhesives on concrete structures such as small dams and walls along the road to allow them to be installed as quickly as possible. In addition to the concrete structure, single-component seismographs were installed in soil such as that in a shoulder along a road. The "Manten" system recorded waveform data continuously at 250 Hz and the others at 100 Hz.
Temporary observation-station data were merged with neighboring permanent-station data, and then event detection and automatic processing were performed using automatic processing software developed by Horiuchi et al. (2009). About 40,000 earthquakes were automatically processed during the period from October 22 to December 15. Earthquakes with 30 or more automatically processed P-arrival times were manually checked and used in the following analysis.

Methods
Only the manually picked data for the stations with epicenter distances within 50 km were used for the hypocenter determination. The hypocenter determination program we used is hypomh_ps (Hirata and Matsu'ura 1987; Kawanishi et al. 2009), in which P-and S-wave velocity structures can be specified independently. The one-dimensional velocity structure used was that estimated in the aftershock area of the 2000 Western Tottori Prefecture earthquake (Shibutani et al. 2005), which was adjacent to the current study area. We first determined hypocenters without station corrections. Then, we calculated the averages of O-C values for P-and S-wave arrival times at each station for earthquakes that had ≥ 30 manually picked P-wave arrival time data and used them as station correction information. While calculating the averages of the O-C values, the values > 0.5 s for P-waves and > 1.5 s for S-waves were omitted. Most errors in the hypocenter calculation are a few tens of meters in the horizontal direction and about twice that in the vertical direction.
Next, the focal mechanisms for aftershocks were determined using an algorithm presented by Maeda (1992) and described in detail by Iio et al. (2017). They involved a fault-plane solution with the least number of inconsistent polarities being extracted via a grid search with a null-axis spacing of ~ 8° in the dip direction. The two nodal planes were then estimated for each null axis by Fourier analysis (Aoki 1986).
In this study, we used only the focal mechanisms for earthquakes with a magnitude greater than 0.0, P-wave polarities of 15 or more, 10 or fewer fault-plane solutions, and focal mechanisms with an uncertainty of < 60° (estimated from the Kagan angle (Kagan 1991) between the most different two solutions). The number of focal mechanisms that survived the above screening process was 10,066, about half of which had only one or two faultplane solutions, and the average uncertainty of the focal mechanisms is 13.2°.
In this study, we selected one optimal fault-plane solution from multiple solutions for each event and analyzed the P-axis and T-axis of the solution. If all the best solutions for an event had no inconsistent polarities, we selected the solution with the nodal planes that were the furthest from the polarity data on the focal sphere; however, if the best focal mechanisms had inconsistent polarities, we selected a solution with a nodal plane closest to the inconsistent polarities.

Results
The hypocenter distribution of aftershocks for which focal mechanisms were determined is shown in Fig. 1. Numerous aftershocks occur even in layers as shallow as 1-5 km. The distribution is about 2 km shallower than that estimated using only permanent-station data (e.g., Kubo et al. 2017). The rectangles in the figure show the mainshock fault estimated from Interferometric Synthetic Aperture Radar (InSAR) and Global Navigation Satellite System (GNSS) data assuming uniform slip (GSI 2017); this is the simplest assumption and can be used as a rough reference model. The aftershock distribution is basically consistent with this model except in deeper regions where such a geodetic fault model does not have good resolution. Figure 2 is an example of the distributions of the P-and T-axis orientations of the optimal focal mechanisms. In Fig. 2a, P-and T-axes of aftershocks located at a depth of 4 ± 0.5 km were indicated by bars projected onto a horizontal plane. We found that there were many aftershocks near both horizontal ends of the fault. Both P-and T-axis azimuths are different between on the right and on the left of the mainshock fault around those ends, suggesting that the magnitude of the differential stress is small compared to the coseismic stress change. The plunge of each axis is indicated by the color of the bar. We also observed many normal-fault type focal mechanisms with large P-axis plunges located in the region with Y = 2-4 km and − 2-0 km. The average azimuth of T-axes that plunge < 30° is N199° E; this is consistent with the stress field in the seismic zone (Kawanishi et al. 2009). Figure 2b shows the projection of P-axes on the vertical cross section along the estimated fault. Only the earthquakes located within 0.5 km of the cross section are shown. It shows that normal-fault type events are conspicuous in some regions, particularly, near the edge of the fault in the S20° E side from the shallow to deeper regions. These normal-fault type events also suggest that the magnitude of the differential stress is small compared with the coseismic stress change.
Next, we examined temporal changes in the hypocenter and focal mechanism distributions. Such temporal changes were reported for the several days to several weeks after the mainshock in various regions. Therefore, we investigated the difference between the two periods of October 22 to 29, and October 30 to December 15, by dividing the entire period into two periods with equal logarithmic lengths. Furthermore, we divided the analysis area into two parts, east and west, at X = − 4.25 near the fault. Distributions of T-axis azimuths are projected on the vertical cross section along the fault, for aftershocks that occurred in the eastern and western sides of the fault over the two periods, in Fig. 3a-d, respectively. T-axis azimuths from 90° (East) through 180° (South) to 270° (West) are represented by solid colors that change from dark blue to green to brown. T-axis azimuths close to the average (N199° E) are represented by green to yellow and extreme ones away from it by dark blue or brown. It is seen that T-axis azimuth color patterns in the eastern and western sides are horizontally symmetrically inverted as a result of the coseismic stress change. In the analysis of temporal change of the focal mechanisms, we set the lower limit of magnitude to 1.0 since it was difficult to precisely determine the focal mechanisms for small earthquakes in October, owing to the small number of observation stations that were in operation.
It seems that the hypocenter distributions of the aftershocks are basically similar for both periods on each side, despite the fact that the exact locations of the individual aftershock hypocenters differ. Furthermore, in terms of the temporal changes in focal mechanisms, the spatial Iio et al. Earth, Planets and Space (2020) 72:42 Fig. 2 Orientation of P-and T-axes of optimal focal mechanisms. a P-and T-axes of aftershocks located at 4 ± 0.5 km projected onto the horizontal plane. b P-axes of aftershocks located within 0.5 km from the mainshock fault projected onto the vertical cross section along the fault. The plunge of each axis is indicated by axis color. Rectangles show the mainshock fault estimated from InSAR and GNSS data assuming uniform slip (GSI 2017) distributions of T-axis azimuths are also basically similar between the two periods, as seen by the distributions of the colored circles, which have a similar pattern in the two periods. If fluid rises from the lower crust due to the mainshock fault rupture, aftershock hypocenters and/ or focal mechanisms may change over time, particularly in the deepest part of the aftershock region, where the greatest temporal changes in fluid pressure are expected to occur. This is because pore pressure is thought to be lithostatic in the lower crust, but hydrostatic in the upper crust. In this case, if fluid rises from the lower crust due to fault rupture, the change in fluid pressure is at its maximum at the bottom of the upper crust where it was probably hydrostatic before the mainshock. However, the temporal changes in the distribution of aftershock hypocenters and T-axis azimuths are not apparent at depth. Only slight decreases of aftershocks with large T-axis azimuths can be seen in the shallow parts (depths of about 1-4 km) on the western side as shown in Fig. 3d. Thus, T-axis azimuths at depths of 0-5 km are projected on the horizontal plane in Fig. 3e, f. It seems that the decreases of the aftershocks with the large T-axis azimuths occurred around the points X = − 5 and Y = 7, and X = − 5 and Y = − 1, near the western edge of the aftershock distribution.
To quantitatively examine the temporal changes in focal mechanisms, we conducted statistical analyses using the Kolmogorov-Smirnov test (KS-test). The KS-test is a non-parametric method that can test whether two datasets have different probability distributions by comparing cumulative frequency distributions. The cumulative frequency distributions of the T-axis azimuths are plotted for the two periods for 18 subregions in Fig. 4. The aftershock area was divided into 9 subregions 5 km by c Those for aftershocks in the former period in the western side of the fault. d Those for aftershocks in the latter period in the western side of the fault. e Those projected on the horizontal plane for aftershocks in the former period at depths of 0-5 km. f Those for aftershocks in the latter period at depths of 0-5 km. T-axis azimuths from 90° (East) through 180° (South) to 270° (West) are represented by solid colors that change from dark blue to green to brown. Rectangles show the mainshock fault estimated from InSAR and GNSS data assuming uniform slip (GSI 2017) Iio et al. Earth, Planets and Space (2020) Fig. 4a-i, and aʹ-iʹ, respectively. The distributions for the P-axis azimuth, P-axis plunge, and T-axis plunge are shown in Additional file 1: Figure S1a-c. For almost all the subregions, the null hypothesis that the two distributions are identical could not be rejected with 5% significance level, while in the three subregions indicated by arrows, the null hypothesis could be rejected. It is likely that the result of Y = 5-10 and Z = − 5-0 shown in Fig. 4c is reliable due to the large amount of data. For P-axis azimuths shown in Additional file 1: Figure S1a, the data of not only Y = 5-10 and Z = − 5-0, but also Y = − 5-0 and Z = − 5-0 show a different distribution, although the number of data in the latter case is not sufficient. These results are consistent with those shown in Fig. 3, as the decrease of the aftershocks with large T-axis azimuths occurred in the shallow layer in the western side.
Thus in the analysis of the aftershocks of the Central Tottori Prefecture earthquake, we found that the distributions of the hypocenters and T-axis azimuths are basically similar for the two periods, however, in the limited portions in the shallow layer near the western end of the aftershock area, a decrease of the aftershocks with large T-axis azimuths were seen in the later period.
Since in the above we have examined the temporal change of aftershock activity for the two periods, we next continuously examined the temporal change of the T-axis azimuths. Figure 5 shows the change in the T-axis azimuths over time. The vertical axis is set as the distance Y in the direction along the fault in Fig. 5a, b, and as the depth Z in Fig. 5c, d. The time axis is logarithm of the number of days that elapsed after the mainshock. The two periods used in the above are divided at the dotted line at the center of the horizontal axis. As shown in Figs. 3 and 4, in Fig. 5a-d, the results in the east and west sides of the fault are plotted separately in the upper and lower panels, respectively. Figure 5e shows the temporal change in the shallow region including the subregion in which the temporal change is detected as shown in Fig. 4c. The direction X is set perpendicular to the fault as is the vertical axis in Fig. 5e. A basically stable pattern can be seen in all panels. In detail, however, in each panel, the aftershock area expanded in the first 2 days (up to about  Fig. 4c is displayed in e with the distance X in the direction perpendicular to the fault on the vertical axis. The arrow approximately indicates the location of the estimated fault 0.3 on the logarithmic axis), and the T-axis azimuths in the expanded area tended to be large (reddish), although there are few data. These features can be observed near the northern end of the aftershock region, in the shallow layer, and in both sides in the direction away from the fault, in Fig. 5a-e, respectively. However, in Fig. 5b, c, it appears that the expanded activities disappeared with time and returned to their original state in the second period. Furthermore, as seen near Y = 7 in Fig. 5a and near Z = − 12 in Fig. 5c, the activity may continue narrowly instead of expanding. The temporal change shown in Fig. 4c corresponds to the fact that the reddish points decreases in the second period around the shallow layers of Z = − 3-0 in Fig. 5c and around X = − 5 in Fig. 5e. Temporal changes in the three-dimensional aftershock distributions are shown in Additional file 1: Figure S2.
It is inferred from the spatial pattern of the stress change due to the mainshock slip that aftershocks with large T-axis azimuths in the regions where the temporal changes were detected were not explained by increases in Coulomb failure stress (CFS) on their fault planes due to the mainshock slip. In the first place, it was seen that the decreases of aftershocks with large T-axis azimuths occurred near the western edge of the aftershock region, where the magnitude of the coseismic stress change was relatively small. Changes in seismic activities that may not be attributed to stress changes are generally interpreted as suggesting high pore fluid pressures, however as noted above, seismic activity shows an expansion only at the very beginning in the limited portions, then begins to shrink there, which is a contrast to what happens in seismic swarms, where there are many cases when the hypocenters migrate. Even in the regions where it can be seen that aftershock activities occurred regularly, the activity appeared to continue in a rather narrowing sense. As the aftershock activity basically looks steady in the logarithmic time axis, it is presumed that the obtained results do not indicate that the aftershock is mainly controlled by fluid diffusion. However, it is currently unknown what happened in the region where the temporal change is detected.

Discussion and conclusion
The relationship between earthquake occurrences and fluids is an important problem in seismology and has been attracting attention in recent years. In the inland region of the Tohoku district, the relationship between earthquake occurrence and fluid migration from the lower crust after the 2011 Tohoku-oki earthquake was demonstrated empirically via the temporal change in the focal mechanisms and stress drops (Yoshida and Hasegawa 2018). However, in this study, we only detected temporal changes in hypocenters and T-axis azimuths only in the limited portions in a shallow section near the edge of the aftershock area; clear temporal changes were not seen in other areas. It seems that aftershock activity did not spread over time, except for the very beginning, and that the activity may continue in a narrow way than expanding. Our findings suggest that the aftershock activity of the Central Tottori Prefecture earthquake was not mainly controlled by high fluid pressure.
Focusing the tiny temporal changes at the beginning to explore different possibilities, we notice the possibility that the fluid that rose from depth immediately after the mainshock (about 2 days) stayed longer in the aftershock area for 2 months or more. On the contrary, it might also be possible that high-pressure fluid migrated rapidly after the first 2 days and disappeared in about a few weeks in the aftershock area, resulting in the decrease of the aftershocks with large T-axis azimuths in the shallow layer on the western side. If the data are reviewed from that perspective, the tendency might be seen in Fig. 5 that the reddish points decrease with time but green to yellow ones increase, indicating that aftershocks with extreme T-axis azimuths away from the average (N199° E) decrease with time while normal ones closer to it increase.
However, it has been pointed out that the focal mechanisms of the aftershocks immediately after the 2000 Western Tottori Prefecture earthquake were very similar to those that occurred about 20 years later (Yukutake and Iio 2017;Iio et al. 2019). This suggests that the aftershock activity may not be controlled by high-pressure fluid migration but by stress in the San'in seismic zone. To fully examine the changes, it is necessary to analyze aftershock data more extensively in a longer period, as they are always controlled by the mainshock slip and the effect of fluid could be masked by the effect of the slip.