We investigated the foreshock activity characteristics using the Japan Meteorological Agency Unified Earthquake Catalog for the last 20 years. Using the nearest-neighbor distance approach, we systematically and objectively classi- fied the earthquakes into clustered and background seismicity. We further categorized the clustered events into fore - shocks, mainshocks, and aftershocks and analyzed their statistical features such as the b-value of the frequency–mag- nitude distribution. We found that the b-values of the foreshocks are lower than those of the aftershocks. This b-value difference suggested that not only the stochastic cascade effect but also the stress changes/aseismic processes may contribute to the mainshock-triggering process. However, forecasting the mainshock based on b-value analysis may be difficult. In addition, the rate of foreshock occurrence in all clusters (with two or more events) was nearly constant (30–40%) over a wide magnitude range. The difference in the magnitude, time, and epicentral distance between the mainshock and largest foreshock followed a power law. We inferred that the distinctive characteristics of foreshocks can be better revealed using the improved catalog, which includes the micro-earthquake information. Keywords: Foreshocks, Aftershocks, b-value, Nearest-neighbor distance, JMA unified catalog accompanied by M6–7 foreshock activities (e.g., Kato Introduction et al. 2012, 2016; Marsan and Enescu 2012). For the Foreshocks are earthquakes that occur prior to the main- case of the 2016 Kumamoto sequence, the Japan Mete- shock, which is defined as the largest magnitude event orological Agency (JMA) issued aftershock probabili- in an earthquake sequence. The foreshock activity often ties following the April 14 M6.5 earthquake based on shows significant differences compared with the ordinary the Earthquake Research Committee (ERC) protocol seismic activity such as a decreased b-value (e.g., Suye- (ERC 1998, 2016). However, a larger earthquake, the hiro 1966; Enescu and Ito 2001; Nanjo et al. 2012; Tor- M7.3 Kumamoto Earthquake (the mainshock), occurred mann et al. 2015) and migration and acceleration prior to on April 16. This case illustrates that it is impossible to the mainshock (e.g., Kato et al. 2012; Marsan et al. 2014). determine whether an earthquake is a foreshock before While these foreshock characteristics may reflect the the occurrence of the mainshock. However, by statisti- accumulation of stress and/or the occurrence of slow slip cally analyzing past foreshock sequences, researchers within the seismogenic zone, Helmstetter et al. (2003) have attempted to probabilistically forecast the occur- showed that such features may also be explained by sim- rence of a larger event, the mainshock (Jones 1985; Aber- ply using the ETAS statistical model of seismicity (Ogata crombie and Mori 1996; Maeda 1996; Reasenberg 1999; 1988) and Gutenberg–Richter law (G–R law). Thus, they Ogata and Katsura 2012). In many of these studies, rela- may reflect stochastic rather than physical processes. tively large magnitude thresholds were used, for example, Two recent large earthquakes, the 2011 M9.0 earth- M ≥ 4 (Ogata and Katsura 2012), due to the incomplete quake off the Pacific coast of Tohoku (Tohoku-oki) detection of smaller events. Because the foreshock activ- and the 2016 M7.3 Kumamoto earthquakes, were ity has a wide range of magnitudes and diverse patterns, it is important to detect smaller events and poten- *Correspondence: email@example.com tially reveal foreshock patterns over a wide magnitude Meteorological Research Institute, Japan Meteorological Agency, Tsukuba, Japan Full list of author information is available at the end of the article © The Author(s) 2018. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 2 of 13 range (such as M0–6) to understand the foreshock Data and method characteristics. We used the JMA unified catalog and targeted the shal - The JMA published and routinely updates an earth - low (depth ≤ 30 km) inland seismicity in Japan, which quake catalog (JMA unified catalog) in cooperation with has a near-homogeneous completeness magnitude of the Ministry of Education, Culture, Sports, Science and approximately 1.0 (e.g., Nanjo et al. 2010). Note that we Technology (MEXT) using seismic waveforms of stations estimated the magnitude of completeness before deter- belonging to the JMA, National Research Institute for mining (and discussing) the b-values of foreshocks and Earth Science and Disaster Resilience (NIED), universi- aftershocks to avoid any possible bias, as explained in ties, and other institutes. The publication of the JMA uni - detail at the end of this section. The data span over the fied catalog started in October 1997, and the data have period from October 1, 1997, to July 25, 2017 (M ≥ 1; been archived for approximately 20 years. Figure 1 shows long-term period). We also investigated the period from the number of events in the catalog from October 1997 to April 1, 2016, to July 25, 2017 (M ≥ 0; short-term period) July 2017. The number of detected events increased since to confirm the effect of the improvement of the JMA uni - 2000 due to the deployment of the Hi-net NIED network fied catalog by automatic processing (Tamaribuchi et al. (Okada et al. 2004). However, M ≥ 1 earthquakes have 2016). This automatic processing method determines the been recorded constantly. Visual inspection of all events hypocenters of earthquakes that occur simultaneously by was impossible after the 2011 Tohoku-oki Earthquake searching for the optimal combination of P- and S-wave because a large number of earthquakes occurred in Japan, arrival times and maximum amplitudes using a Bayesian even in inland areas such as the northern part of the Iba- estimation technique. The system detects nearly twice raki Prefecture. Therefore, the M < 2 earthquakes have as many earthquakes as those listed in the conventional not been completely cataloged in the aftershocks area of JMA catalog. By effectively using this catalog, the seismic the 2011 Tohoku-oki Earthquake. In April 2016, the JMA activity can be investigated over a wide magnitude range. adopted a new automatic processing methodology, devel- First, we clustered the seismic activity to extract fore- oped by Tamaribuchi et al. (2016) to more accurately and shocks and aftershocks from numerous hypocenters. quickly determine small earthquakes; thus, the number Several automatic clustering methods have already been of M < 1 quakes significantly increased (Fig. 1b; see details proposed (e.g., Frohlich and Davis 1990; Ogata and Kat- in JMA, 2017). In this study, we report the characteris- sura 2014), which are suitable for the analysis of big data- tics of foreshocks and aftershocks over a wide magnitude sets. In this study, we used the nearest-neighbor distance range using a long-term interval (1997–2017). method (Baiesi and Paczuski 2004; Zaliapin and Ben- Zion 2013), which can clearly and objectively distinguish the clustered and background earthquakes. We used a M all a M ≥ 7 b M ≥ 1 M ≥ 2 C D A B Year Fig. 1 JMA unified catalog for 20 years starting in October 1997. a Epicentral distribution (October 1, 1997, to July 25, 2017, M ≥ 1). Gray dots indi- cate epicenters in the inland area (depth ≤ 30 km). Red stars indicate M ≥ 7 earthquakes (also shown as A, B, C, and D). b Frequency–time diagram by month. The blue, yellow, and red bars indicate the entire magnitude, M ≥ 1, and M ≥ 2, respectively Number of events Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 3 of 13 program by Kasahara (2016) for the analysis, which is We estimated the b-values of the clusters using the faster than the general implementation [the calculation maximum likelihood method (Aki 1965; Utsu 1965). To cost is reduced from O N to O N N ] and allows the calculate the standard error of the b-values, σ , we used method to be applied to large earthquake catalogs. The the equation of Shi and Bolt (1982). We determined the seismic activity was clustered by defining a new metric completeness magnitude, Mc, based on the “MAXC” (called “distance”) between events and determined using approach, proposed by Wiemer and Wyss (2000), which the epicentral distance, time difference between event considers the highest frequency of events in the noncu- pairs, and magnitude information. The new metric is mulative frequency–magnitude distribution. The magni - defined by the following equation: tude of completeness is defined as MAXC + 0.2 according to Woessner and Wiemer (2005). When b-values for dif- f −bm ferent distributions are compared, Mc is first calculated η = t r 10 , t > 0, ij ij ij ij (1) for each dataset, and the largest value is used as the mag- η =∞, t ≤ 0. ij ij nitude threshold. The parameter, η , is the distance between the events, ij i and j; t is the difference in the origin time (year); r is Results ij ij the epicentral distance (km); m is the magnitude of the Figure 2c, d shows the histogram of the nearest-neighbor parent event, i; d and b are the fractal dimensions of distance, η , between earthquakes based on short-term the epicentral distributions of earthquakes and b-value, data; the long-term data are also shown in Fig. 2e, f. The respectively. By dividing the aforementioned equation histogram of both datasets shows a bimodal distribution. −qbm into time and spatial components, that is, T = t 10 The threshold for the long-term data, log (η ) =−4.2 , ij ij 0 −(1−q)bm and R = r 10 , respectively, we can also is slightly larger than the threshold of Zaliapin and Ben- ij ij describe Eq. (1) as η = T R . q is a free parameter. Each Zion (2013) of log (η ) =−5.0 . This trend was also ij ij ij 0 parameter was set by following Zaliapin and Ben-Zion observed for the short-term data ( log (η ) =−4.5 ). W e (2013), that is, d = 1.6, b = 1.0, and q = 0.5. Based on this also experimented with other values for b and d in Eq. (1) f f algorithm, the distance, η , was calculated for each event, of the short-term data but did not show significant varia - ij j, and the shortest distance, η , between the parent event, tions in the results (Table 1). ij i, and the child event, j, was taken as the nearest-neigh- We extracted 58,509 clusters from the analysis of bor distance. 241,946 events for the short-term data. The clusters Second, based on Zaliapin and Ben-Zion (2013), we included 6103 families with multiple events in a cluster defined a threshold, η , that distinguishes whether two and 52,406 singles, which comprise one event. events are strongly or weakly linked. In other words, in An example of an extracted cluster for the 2016 M7.3 the case of η < η , the event pair is strongly linked (clus- Kumamoto Earthquake sequence is shown in Fig. 3. tered), while the event pair is weakly linked in the case of (Another example for the 2016 M6.6 central Tottori pre- η ≥ η . (The events may belong to background seismic - fecture earthquake is shown in Additional file 1: Fig. S1.) ity.) Weakly linked events are well described by a Poisson If the magnitude of the parent event (e.g., mainshock) is (random) process. The histogram of η was fitted assum - large, subsequent earthquakes that occur shortly after ing two mixed Weibull distributions for the strongly and the parent event may be included in the same cluster as weakly linked seismicity (Fig. 2a, b). The threshold, η , is their parent, as child events, due to the characteristics of defined as the intersection of the two fitted distributions. Eq. (1), although they occur at very large distances from By sequentially linking strongly related pairs, we can the parent event (e.g., in Hokkaido for the Kumamoto build clusters of events. For each cluster, we defined the mainshock). However, the events in the blue rectangle foreshocks, mainshock, and aftershocks in the same way in Fig. 3a dominate the sequence over the whole period, as in conventional studies. The mainshock is defined as that is, 98,682 events (97.2% of 101,489 events). In the the earthquake of largest magnitude in each cluster, while period from April 14 to 18 (same time span as in Fig. 3), the foreshocks and aftershocks are earthquakes occur- 9868 events (94.3% of 10,465 events) are observed in the ring before and after the mainshock, respectively. If there rectangle. Therefore, fewer distant events occurred out - are two or more earthquakes of the same largest magni- side the blue rectangle. Notably, 83,280 events (82.1% of tude, the one that occurred earlier is considered to be the 101,489 events) are hypocenters automatically processed mainshock. Events with η ≥ η that have no child events by the method of Tamaribuchi et al. (2016). As already are referred as singles. Clusters that have a mainshock, shown in Fig. 1, this automatic process substantially con- foreshocks, and/or aftershocks are referred as families; tributes to the improvement of the JMA unified catalog. that is, one “family” has two or more events in contrast to We reconfirmed the b-value difference between a “single,” which comprises one event. foreshocks and aftershocks for the 2016 Kumamoto Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 4 of 13 a Schematic view b Schematic view Weakly related Strongly Weakly (Background) related related (Clustered) (Background) Strongly related (Clustered) log(η) log(T) log(η0) False Positive (FP) False Negative (FN) c Short-term data d Short-term data e Long-term data f Long-term data Frequency log(R) Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 5 of 13 (See figure on previous page.) Fig. 2 Distance measure using the nearest-neighbor distance analysis. a, b Schematic view, c, d the short-term data, and e, f the long-term data. a, c, e log η histogram. The gray area shows the frequency of log η , the blue dotted line shows the optimal solution when assuming two Weibull distributions, and the red line shows the sum (theoretical value) of the mixed Weibull distribution. b, d, f Joint distribution with rescaled time T and rescaled distance R. The dashed line shows η Table 1 Parameter validity check Figure 4a shows the frequency–magnitude distribution of the dataset that includes both foreshocks and after- b d η FP (%) FN (%) f 0 shocks, while the foreshocks and aftershocks in Fig. 4b 1.0 1.6 − 4.5 1.9 3.3 were separately calculated. The b-value is 0.62 ± 0.02 for 0.8 1.6 − 3.3 3.6 3.9 foreshocks and 0.68 ± 0.01 for aftershocks in the case of 1.2 1.6 − 5.7 2.5 2.4 Mc = 2.0; the b-value of the foreshocks is slightly lower 1.0 1.4 − 4.7 2.1 3.9 than that of the aftershocks. We confirmed whether this 1.0 1.8 − 4.3 2.6 2.9 b-value difference was significant by using the ΔAIC test (Utsu 1999). The AIC (Akaike 1974) is often used as an The input parameters, b and d , take the values as shown above. The meaning of FP and FN is explained in Fig. 2a indicator of comparison between models. One model (model 1) assumed that all earthquakes (foreshocks and aftershocks) have the same b-value, while the other Earthquake sequence. Since the aftershock period con- model (model 2) assumed that foreshocks and after- tinues to the present, we decided to calculate the b-value shocks have their own (different) b-value. We calculated of aftershocks using the data from the mainshock until the AIC values of two models, and we obtained the dif- April 2017 (Fig. 4). However, the mainshock is not ference of their AIC; ΔAIC = AIC(model 1) − AIC(model included in the dataset of the foreshocks or aftershocks. 2). The b-value difference is considered to be statistically Apr. 16, 2016 01:25 M7.3 Oita Kumamoto Date (April 2016) Fig. 3 Space–time–magnitude distributions of the 2016 Kumamoto earthquake sequence. a Epicentral distribution, b spatiotemporal distribu- tion, and c magnitude versus time diagram. The light red open circles, red squares, and dark stars indicate foreshocks, aftershocks, and mainshock, respectively. The gray dots represent other clusters or singles Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 6 of 13 ab Foreshock Cum. Aftershock Non- cum. Mc cd Fig. 4 Frequency–magnitude distributions of the 2016 Kumamoto earthquake sequence. a, b Frequency magnitude distributions and c, d b-value versus magnitude: a, c for the entire sequence (foreshocks and aftershocks), mainshock excluded, and b, d for the foreshocks and aftershocks, separately. The red and blue colors indicate foreshocks (b ) and aftershocks (b ), respectively. The frequency–magnitude distributions show both the 1 2 cumulative (black) and noncumulative number (light gray) of events. The dashed black line indicates Mc estimated by the MAXC method significant if ΔAIC > 2. According to the ΔAIC test, this consideration. They adopted the time-dependent func - difference is significant (ΔAIC = 3.0; Fig. 4). However, the tion of the detection rate and estimated the b-value more b-value seems to depend on the lower magnitude thresh- robustly. We used the data starting 1 day after the April old (Mth) as shown in Fig. 4c, d. Immediately after a large 14 M6.5 earthquake for the foreshocks and that from earthquake, such as the Kumamoto Earthquake (M7.3), the April 16 M7.3 earthquake (mainshock) for the after- many small events are usually missing; thus, Mc = 2.0, shocks. In this case, the mainshock is included in the estimated by MAXC method, might underestimate the aftershock dataset. The resulting b-value is 0.65 ± 0.19 completeness magnitude shortly after the mainshock. for the foreshocks and 0.98 ± 0.22 for the aftershocks; Therefore, we also confirmed the b-value for the cases thus, the b-value of the foreshocks is lower than that of of Mc = 3.0 or 4.0. Even in these cases, the b-value of the the aftershocks as reported by Nanjo and Yoshida (2017) foreshocks is lower, that is, 0.74 ± 0.06 for foreshocks and ERC (2016). We also mention that the same b-value compared with 0.86 ± 0.03 for aftershocks (M ≥ 3.0) is obtained even if the mainshock is excluded from the and the b-value of 0.68 ± 0.13 for foreshocks compared aftershock dataset. with 0.79 ± 0.09 for the aftershocks (M ≥ 4.0). However, To extract detailed seismicity characteristics for the we cannot conclude that there is a statistically signifi - whole catalog, we classified all earthquakes of the long- cant difference according to the ΔAIC test (ΔAIC = 1.6 term data as families (mainshocks, foreshocks, and after- for M ≥ 3.0, and ΔAIC = -1.6 for M ≥ 4.0). To confirm shocks) and singles. Table 2 summarizes the classification the obtained b-values, we also used the method of Omi results. Among all the detected earthquakes, 12% were et al. (2013), which takes the temporal change of Mc into foreshocks and 62% were aftershocks. There are four b value Number of events Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 7 of 13 Table 2 Statistics of singles, mainshocks, aftershocks, and foreshocks using the nearest-neighbor distance analysis of long-term data Magnitude range Singles Families All Mainshocks Aftershocks Foreshocks All events: M ≥ 0 93,938 22% 17,700 4% 265,656 62% 52,332 12% 429,626 1 ≤ M < 2 87,014 24% 9990 3% 219,297 61% 41,847 12% 358,148 2 ≤ M < 3 6676 11% 5934 10% 39,884 65% 8960 15% 61,454 3 ≤ M < 4 245 3% 1480 17% 5674 65% 1315 15% 8714 4 ≤ M < 5 3 0% 248 22% 718 62% 184 16% 1153 5 ≤ M < 6 0 0% 35 27% 76 58% 20 15% 131 6 ≤ M < 7 0 0% 9 41% 7 32% 6 27% 22 7 ≤ M 0 0% 4 100% 0 0% 0 0% 4 sequences with a M7 mainshock: (A) the M7.3 West- of both aftershocks and foreshocks tend to be lower ern Tottori Earthquake on October 6, 2000; (B) M7.2 than the background b-value. In case of the large main- Iwate–Miyagi Nairiku Earthquake on June 14, 2008; (C) shocks (M ≥ 5.5), the lack of catalog data immediately main M7.0 Eastern Fukushima Earthquake on April 11, 2011; after the mainshock may lead to b-value underestima- and (D) M7.3 Kumamoto Earthquake on April 16, 2016 tion. Therefore, we also estimated the b-value by fixing (Fig. 1). The M7.3 Western Tottori Earthquake did not Mc = M − 3.5 for M ≥ 5.5. Under this condition, main main show a foreshock activity, while the M7.2 Iwate–Miyagi the tendency toward a low b-value of the foreshocks Nairiku Earthquake had a M1.3 foreshock at 8:11, only does not change, although the number of events (≥ Mc) 30 min before the mainshock. In the case of the M7.0 decreases (Fig. 5b and Table 3b). However, when cal- Eastern Fukushima Earthquake, the foreshock activity culating the b-value of aftershocks (≥ 200 events) for started with the M5.7 quake on March 11, shortly after sequences without foreshocks or with a small number the 2011 Tohoku-oki Earthquake. As already mentioned, of foreshocks (< 200 events), the variation is quite large two large foreshocks of M > 6.0 were recorded in the case (Fig. 5c–f ), and thus, it would be difficult to judge from of the 2016 Kumamoto Earthquake sequence, that is, the the b-values alone, in quasi-real time, if it is a foreshock April 14 M6.5 and April 15 M6.4 earthquakes. or an aftershock sequence. The anomalously high b-value We further extracted both foreshocks and aftershocks (Δb > 1.0, Group ID #6 in Table 3a) is associated with of families with 200 or more events and investigated the swarm-like activities around the Mt. Sharitake volcano in b-value difference between the foreshocks and after - the eastern part of Hokkaido. Notably, this b-value is out shocks. Aftershocks usually span over a considerably of range in Fig. 5a, b. longer period compared to the foreshock activity. To Earthquake forecasting applications would benefit focus on the b-value characteristics close to the time from examining how often the mainshock would occur of the mainshock occurrence, we chose a time win- (and its magnitude, timing, and location) when assuming dow of 200 events (for which we determined Mc and b). that ongoing clustered seismicity is a foreshock activity. Because b-values likely depend on the tectonic region Therefore, the percentage of families with a foreshock in which the seismic activity takes place, we compared activity was determined to verify how many mainshocks the b-values of the foreshocks and aftershocks with the were accompanied by foreshocks. Figure 6a shows that b-value of the background seismicity. The b-value of the the percentage of families (without singles) with a fore- background seismicity was calculated using hypocenters shock activity is approximately 30–40% in the range of other than the target cluster within a range of ± 0.4° from M1–5 for the mainshocks and approximately 70% in the the mainshock’s epicenter of the target cluster. The com - range of M5–7. When focusing on at least 100 samples, parison of the b-value differences for different magni - the percentage (rate) of foreshock occurrence appears tudes of the mainshock is shown in Fig. 5a and Table 3a. to be independent of the mainshock magnitude (M1–5). Most b-value differences are insignificant according The number of singles increases for smaller magnitudes. to the ΔAIC test; however, the overall tendency does Therefore, if singles are included in the calculations, the not depend on the mainshock magnitude (M ), and percentage of foreshock occurrence approaches zero in main the b-value of the foreshocks tends to be relatively low. the small magnitude range (see Fig. 6a, dashed red line). When the mainshock magnitude increases, the b-values We also observed that the frequency distribution of the Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 8 of 13 ab Mc=MAXC+0.2 Mc=MAXC+0.2 Mc=M –3.5 main Foreshock–Background Foreshock–Background cd Aftershock–Background Aftershock–Background Magnitude of mainshock Magnitude of mainshock ef Δb Δb Fig. 5 Comparison of b-value differences. a–d Comparison of the magnitude of the mainshock and b-value differences. The magnitude of the mainshock is shown on the horizontal axis; the vertical axis shows the b-value difference. a, b Foreshocks–aftershocks (with 200 foreshocks and 200 aftershocks), c, d foreshocks—background (with 200 foreshocks, red circle), and aftershocks–background (with 200 aftershocks, black circle). e, f His- togram of b-value differences. The frequency of the b-value differences is on the left axis, and the cumulative probability is shown on the right axis. a, c, e Mc was chosen as MAXC + 0.2 for the entire magnitude range, b, d, f Mc was chosen as MAXC + 0.2 for M < 5.5; M − 3.5 for M ≥ 5.5. main main main The b-value of the background was calculated using hypocenters other than the target cluster within the range of ± 0.4 degree of the mainshock of the target cluster magnitude difference between the mainshock and largest (M ≥ 1.0, Fig. 7a), and the time difference between the foreshock (Fig. 6b) follows a power law; that is, the main- largest foreshock and mainshock ( t ) attenuates as 1/�t shock magnitude tends to be close to the magnitude of between 0.1 and 10 days (Fig. 7c). However, approxi- the largest foreshock. The frequency of ΔM decreases by mately 1 week is required for the mainshock to occur in 1/10 per unit. 50% of the cases of M ≥ 2.5 foreshocks. The larger the Figure 7 shows the time difference and epicentral magnitude of the foreshock, the longer is the time dif- distance of the largest foreshock relative to its main- ference between the largest foreshock and mainshock. shock, respectively. Approximately 50% of the main- Because, generally, larger mainshocks have larger fore- shocks occurred within a day after the largest foreshock shocks (see above), a strong link between them at longer Δb Cumulative probability [%] Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 9 of 13 Table 3 Statistics of clusters with 200 or more foreshocks and aftershocks Group Mainshock Mc Foreshock Aftershock Background ΔAIC ID Lat (°) Lon (°) M N (≥ Mc) b σ N (≥ Mc) b σ N (≥ Mc) b σ b b b (a) Mc = MAXC + 0.2 for the entire magnitude range 1 31.84 130.29 4.2 1.2 131 0.95 0.09 140 1.15 0.10 4171 1.04 0.02 0.5 2 34.96 139.18 5.9 1.3 153 0.50 0.03 155 0.59 0.03 3234 0.86 0.01 0.3 3 36.33 137.63 5.6 1.7 92 0.75 0.08 87 0.90 0.09 2826 0.85 0.02 − 0.6 4 35.28 135.93 5.2 1.3 110 1.02 0.12 105 0.98 0.09 4055 0.89 0.01 − 1.9 5 42.50 140.83 4.9 1.8 65 1.12 0.13 71 1.16 0.11 37 0.80 0.11 − 2.0 6 43.74 144.70 4.8 1.6 42 2.53 0.25 56 1.06 0.13 360 0.93 0.06 14.9 7 37.30 136.84 5.3 2.7 71 0.75 0.08 23 0.83 0.12 19 0.68 0.12 − 1.9 8 34.96 139.13 5.1 1.4 78 1.11 0.10 89 0.87 0.09 5362 0.76 0.01 0.5 9 36.95 140.67 7.0 3.2 30 0.57 0.10 114 0.63 0.05 9 1.04 0.24 − 1.8 10 32.75 130.76 7.3 3.1 50 0.69 0.10 126 0.72 0.06 166 1.00 0.07 − 1.9 11 31.38 130.62 5.3 1.8 27 1.02 0.20 34 1.16 0.15 275 1.15 0.06 − 1.8 (b) Mc = M – 3.5 for M ≥ 5.5 main main 2 34.96 139.18 5.9 2.4 47 0.82 0.11 31 0.96 0.12 3234 0.86 0.01 − 1.5 3 36.33 137.63 5.6 2.1 45 0.74 0.11 40 1.00 0.15 2826 0.85 0.02 0.0 9 36.95 140.67 7.0 3.5 21 0.59 0.12 81 0.71 0.06 9 1.04 0.24 − 1.4 10 32.75 130.76 7.3 3.8 15 0.62 0.18 42 0.89 0.15 166 1.00 0.07 − 0.5 Lat., Lon., and M indicate the latitude, longitude, and magnitude of the mainshock, respectively. Mc indicates the completeness magnitude. The parameters, b and σ , are the b-value and its uncertainty, respectively. ΔAIC indicates the difference of AIC values (Utsu 1999). a Mc was chosen as MAXC + 0.2 for the entire magnitude range, b Mc was chosen as M − 3.5 for M ≥ 5.5 main main ab Fig. 6 Frequency–magnitude distribution and relationship between the largest foreshock and its mainshock for the long-term data. a Frequency– magnitude distribution and percentage of foreshock occurrence. The white solid and dashed histograms show the frequency of mainshocks versus their magnitude for all families (without singles) and clusters (with singles), respectively; the gray histogram shows the frequency of mainshocks that are accompanied by foreshocks (left axis). The red solid and dashed lines show the proportion of sequences that have foreshocks for all families and clusters, respectively (right axis). b Relationship between the magnitude of the largest foreshock and the corresponding mainshock magnitude for each earthquake cluster. The gray shading indicates M ≤ 1.0 (potentially incomplete data). The lower right inset shows the frequency distribution of the magnitude difference between the mainshock and largest foreshock. The dashed line has a slope of 1.0. The black dots, blue triangles, and red crosses indicate that the magnitude of the largest foreshock is above 1.0, 1.5, and 2.5, respectively Rate of Foreshock Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 10 of 13 a Δt b Δr c Δt d Δr Fig. 7 Cumulative probability and frequency density distribution of the time difference ( t ) and the epicentral distance ( r ) between the largest foreshock and mainshock for long-term data. a, c Time difference ( t ) and b, d epicentral distance ( r ). a, b Cumulative probability distribution. b, d Probability density distribution. The dashed lines have slopes of − 1.0, − 1.5, and − 2.0, as indicated in the figure. The black dots, blue triangles, and red crosses indicate that the magnitude of the largest foreshock is above 1.0, 1.5, and 2.5, respectively time differences is more commonly observed. The distri - Discussion bution of the distance between the largest foreshock and The b-value difference between foreshocks and after - mainshock ( r ) is shown in Fig. 7b. Approximately 90% shocks has been discussed in many studies. Recently, of the mainshocks occur within a range of 1 km from several studies pointed out that the b-value of foreshocks their largest foreshock. The frequency of r attenuates decreases before major subduction-zone earthquakes, as 1/�r (Fig. 7d). The distribution does not strongly such as the 2011 Tohoku-oki and 2014 Iquique earth- depend on the magnitude of the foreshock, in contrast to quakes, due to the stress accumulation before the meg- the time difference. Whether purely statistical in nature athrust events (e.g., Nanjo et al. 2012; Tormann et al. or reflecting some physical generation mechanisms, the 2015; Schurr et al. 2014). The b-value after the 2011 time difference and epicentral distance characteristics Tohoku-oki Earthquake increased relatively rapidly, may be useful for better hazard assessment of large earth- which may reflect the heterogeneous stress recovery pro - quakes (Fig. 7). cess after the megathrust event (Tormann et al. 2015, -1 Prob. Dens. [unit ] Cumulative Probability [%] Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 11 of 13 2016). As Tormann et al. (2016) suggested, the b-values The frequency distribution of the magnitude difference possibly react less to the absolute stress levels and more between the mainshock and largest foreshock follows a to the homogeneity of the stress field, which would have power law (Fig. 6b). Investigation of the physical nature been disturbed during the mainshock rupture and could of this power law in detail is beyond the scope of this likely recover more easily via aftershocks. Figure 5 shows study; however, our results suggest that the size of the not only a b-value difference between foreshocks and mainshock could be estimated from the foreshock activ- aftershocks for the large Kumamoto Earthquake (Fig. 4) ity in a probabilistic way (Jones 1985). As shown in Fig. 7, but also over a wide magnitude range of M3–7 for the the time difference and epicentral distance between entire analyzed catalog. If the b-value is related to dif- the mainshock and largest foreshock can also contrib- ferential stress, as discussed in laboratory and actual ute to the estimation of the probability of mainshock seismicity studies (Scholz 1968; Amitrano 2003; Schor- occurrence. lemmer et al. 2005; Scholz 2015), this would indicate the By further extending the completeness and homoge- stress accumulation before mainshocks. However, from neity of earthquake catalogs, probabilistic forecasting of the point of view of classifying foreshocks and after- earthquakes is expected to improve based on better sta- shocks based on their b-values, our results suggest that tistics and a wider, more complete magnitude range. The it is quite difficult to distinguish them. Further investiga - progress in the development of automatic processing tions are necessary to describe them probabilistically in techniques (e.g., Yoon et al. 2015; Tamaribuchi et al. 2016 more detail. for hypocenter determination and Omi et al. 2013 for If the “stress level” reflects the closeness to a critical b-value estimation) will become increasingly important. threshold of differential stress (Scholz 1968), the b-value change during foreshocks and aftershocks is likely due Conclusion to the spatial heterogeneity of the stress distribution. We analyzed the long-term and wide magnitude range Generally, when the increase in the stress at long wave- of foreshock activity using the JMA unified catalog from lengths is dominant, the rupture tends to expand (i.e., 1997 to 2017. By objectively clustering the inland, shal- larger earthquakes tend to occur more frequently), which low seismicity, we identified foreshocks and aftershocks results in lower b-values. Conversely, when the increase and analyzed their statistical characteristics. Our sys- in the stress at short wavelengths is dominant, the rup- tematic investigation supports the previous findings, ture does not spread smoothly, which results in higher mainly based on the analysis of individual sequences, b-values (Scholz 1968). A decrease in the b-values before that the b-value of foreshocks is slightly lower than that the mainshocks corresponds to the predominance of long of the aftershocks over a wide magnitude range of M3–7 wavelengths, relatively large stress levels, and is caused by (although the differences are rather subtle and, probably, the stress accumulation process before the earthquake. still difficult to use for prospective earthquake forecast - Short-wavelength (i.e., highly heterogeneous) stresses ing). This result suggests that factors, such as the stress become dominant after the mainshock, and the b-value interactions induced either seismically or aseismically, increases due to the stress release for the “wavelength” may also contribute to the mainshock-triggering/nuclea- that corresponds to the mainshock magnitude. Consider- tion process. We also analyzed the characteristics of the ing that the increase in stress at long wavelengths (e.g., magnitude and occurrence time differences, as well as those caused by the Earth tides; Ide et al. 2016) leads to epicentral distance between the mainshock and the larg- a decrease in the b-value and relatively high probability est foreshock and confirmed that they follow a power law. of the occurrence of a large earthquake, stress changes We expect that the characteristics of seismic activity will should be observed at various wavelengths to understand become clearer with further enhancement of the number the preparation process of huge earthquakes. of cataloged earthquakes due to progress in the process- Figure 6a shows that the percentage of clusters (with- ing for automatic hypocenter determination. out singles) with a foreshock activity is approximately Additional file 30–40% in the range of M1–4 for the mainshocks. This percentage (rate) of foreshock occurrence appears to be Additional file 1: Fig. S1. Space–time–magnitude distributions of the independent of the mainshock magnitude. Abercrombie Central Tottori earthquake (from October 21, 2016, to October 23, 2016). and Mori (1996) showed that the rate of foreshock occur- a Epicentral distribution, b spatiotemporal distribution, and c magnitude rence in California is approximately 40% in the range of versus time diagram. The light red open circles, red squares, and dark stars indicate foreshocks, aftershocks, and mainshock, respectively. The gray M5–6 for the mainshocks, without magnitude depend- dots represent other clusters or singles. ence. The result of this study suggests that the rate of foreshock occurrence might be relatively constant (30– 40%) for a wider magnitude range. Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 12 of 13 Abbreviations Earthquake Research Committee (1998) Aftershock probability evaluation AIC: Akaike information criteria; ETAS: epidemic type aftershock sequence; methods. http://www.jishi n.go.jp/main/yoshi n2-e/yoshi n2.htm. Accessed JMA: Japan Meteorological Agency; ERC: Earthquake Research Committee; 15 Sep 2017 MEXT: Ministry of Education, Culture, Sports, Science and Technology; NIED: Earthquake Research Committee (2016) Guidelines for the seismic forecast National Research Institute for Earth Science and Disaster Resilience; G–R law: information after big earthquakes. http://www.jishi n.go.jp/repor ts/resea Gutenberg–Richter law.rch_repor t/yosok u_info/. Accessed 15 Sep 2017 Enescu B, Ito K (2001) Some premonitory phenomena of the 1995 Hyogo- Authors’ contributions Ken Nanbu (Kobe) earthquake: seismicity, b-value and fractal dimen- KT performed the analysis and drafted the manuscript. YY participated in sion. Tectonophysics 338(3–4):297–314. https ://doi.org/10.1016/S0040 the design of the study. BE and SH contributed to interpretation of results -1951(01)00085 -3 and a part of the discussion. All the authors participated in discussions and Frohlich C, Davis SD (1990) Single-link cluster analysis as a method to evaluate contributed to revising an earlier draft of the manuscript. All authors read and spatial and temporal properties of earthquake catalogues. Geophys J Int approved the final manuscript. 100(1):19–32. https ://doi.org/10.1111/j.1365-246X.1990.tb045 64.x Helmstetter A, Sornette D, Grasso JR (2003) Mainshocks are aftershocks of con- Author details ditional foreshocks: how do foreshock statistical properties emerge from Meteorological Research Institute, Japan Meteorological Agency, Tsukuba, aftershock laws. J Geophys Res 108:2046. https ://doi.org/10.1029/2002J Japan. Faculty of Life and Environmental Sciences, University of Tsukuba, B0019 91 Tsukuba, Japan. Department of Geophysics, Graduate School of Science, Ide S, Yabe S, Tanaka Y (2016) Earthquake potential revealed by tidal influence Kyoto University, Kyoto, Japan. Department of Physical Science, College on earthquake size-frequency statistics. Nat Geosci 9:834–839. https :// of Science and Engineering, Ritsumeikan University, Kusatsu, Japan. doi.org/10.1038/ngeo2 796 Japan Meteorological Agency (2017) The seismological bulletin of Japan. Acknowledgementshttp://www.data.jma.go.jp/svd/eqev/data/bulle tin/catal og/notes _e.html. We thank the Editor Takuto Maeda and two anonymous reviewers for their Accessed 1 Mar 2017 helpful comments. The earthquake catalog is produced by the JMA in Jones LM (1985) Foreshocks and time-dependent earthquake hazard assess- cooperation with MEXT. We used a python package (AftFore: https ://githu ment in Southern California. Bull Seismol Soc Am 75(6):1669–1679 b.com/omita kahir o/AftFo re) coded by T. Omi to estimate the b-value. We used Kasahara A (2016) An implementation of Zaliapin and Ben-Zion (2013, JGR). generic mapping tools (GMT; Wessel et al. 2013) to draw the figures. Authors https ://githu b.com/kshra mt/trial _kshra mt. Accessed 6 Jan 2017 acknowledge the partial support from the Japan Society for the Promotion of Kato A, Obara K, Igarashi T, Tsuruoka H, Nakagawa S, Hirata N (2012) Propaga- Science (JSPS) KAKENHI Grant 16K05529 and 16H06477 (to Y. Yagi). tion of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki Earthquake. Science 335:705–708. https ://doi.org/10.1126/scien ce.12151 41 Competing interests Kato A, Fukuda JI, Nakagawa S, Obara K (2016) Foreshock migration preced- The authors declare that they have no competing interests. ing the 2016 Mw 7.0 Kumamoto earthquake, Japan. Geophys Res Lett 43(17):8945–8953. https ://doi.org/10.1002/2016G L0700 79 Availability of data and materials Maeda K (1996) The use of foreshocks in probabilistic prediction along the The data used in this article are available at the Data Management Center of Japan and Kuril trenches. Bull Seismol Soc Am 86(1A):242–254 NIED (http://www.hinet .bosai .go.jp/?LANG=en) and JMA (http://www.data. Marsan D, Enescu B (2012) Modeling the foreshock sequence prior to the jma.go.jp/svd/eqev/data/bulle tin/index _e.html). 2011, Mw9.0 Tohoku, Japan, earthquake. J Geophys Res 117(B6):B06316. https ://doi.org/10.1029/2011J B0090 39 Consent for publication Marsan D, Helmstetter A, Bouchon M, Dublanchet P (2014) Foreshock Not applicable. activity related to enhanced aftershock production. Geophys Res Lett 41(19):6652–6658. https ://doi.org/10.1002/2014G L0612 19 Ethics approval and consent to participate Nanjo KZ, Yoshida A (2017) Anomalous decrease in relatively large shocks Not applicable. and increase in the p and b values preceding the April 16, 2016, M7.3 earthquake in Kumamoto, Japan. Earth Planets Space 69:13. https ://doi. Fundingorg/10.1186/s4062 3-017-0598-2 Authors acknowledge the partial support from the Japan Society for the Pro- Nanjo KZ, Ishibe T, Tsuruoka H, Schorlemmer D, Ishigaki Y, Hirata N (2010) motion of Science (JSPS) KAKENHI Grant 16K05529 and 16H06477 (to Y. Yagi). Analysis of the completeness magnitude and seismic network cover- age of Japan. Bull Seismol Soc Am 100(6):3261–3268. https ://doi. org/10.1785/01201 00077 Publisher’s Note Nanjo KZ, Hirata N, Obara K, Kasahara K (2012) Decade-scale decrease in b Springer Nature remains neutral with regard to jurisdictional claims in pub- value prior to the M9-class 2011 Tohoku and 2004 Sumatra quakes. Geo- lished maps and institutional affiliations. phys Res Lett 39(20):L20304. https ://doi.org/10.1029/2012G L0529 97 Ogata Y (1988) Statistical models for earthquake occurrences and residual Received: 14 January 2018 Accepted: 18 May 2018 analysis for point processes. J Am Stat Assoc 83(401):9–27 Ogata Y, Katsura K (2012) Prospective foreshock forecast experiment during the last 17 years. Geophys J Int 191:1237–1244. https ://doi.org/10.1111/ j.1365-246X.2012.05645 .x Ogata Y, Katsura K (2014) Comparing foreshock characteristics and foreshock forecasting in observed and simulated earthquake catalogs. J Geophys References Res Solid Earth 119(11):8457–8477. https ://doi.org/10.1002/2014J B0112 Abercrombie RE, Mori J (1996) Characteristics of foreshock occurrence to large earthquakes in the western United States. Nature 381(6580):303–307 Okada Y, Kasahara K, Hori S, Obara K, Sekiguchi S, Fujiwara H, Yamamoto A Akaike H (1974) A new look at the statistical model identification. IEEE Trans (2004) Recent progress of seismic observation networks in Japan -Hi-net, Automat Contr 19(6):716–723. https ://doi.org/10.1109/TAC.1974.11007 05 F-net, K-NET and KiK-net–. Earth Planets Space 56:xv–xxviii. https ://doi. Aki K (1965) Maximum likelihood estimate of b in the formula log (N) = a − bM org/10.1186/bf033 53076 and its confidence limits. Bull Earthq Res Inst Tokyo Univ 43:237–239 Omi T, Ogata Y, Hirata Y, Aihara K (2013) Forecasting large aftershocks within Amitrano D (2003) Brittle–ductile transition and associated seismicity: one day after the main shock. Sci Rep 3:2218. https ://doi.org/10.1038/ experimental and numerical studies and relationship with the b value. J srep0 2218 Geophys Res 108(B1):2044. https ://doi.org/10.1029/2001J B0006 80 Reasenberg PA (1999) Foreshock occurrence before large earthquakes. J Geo- Baiesi M, Paczuski M (2004) Scale-free networks of earthquakes and after- phys Res 104(B3):4755–4768. https ://doi.org/10.1029/1998J B9000 89 shocks. Phys Rev E 69(6):066106 Tamaribuchi et al. Earth, Planets and Space (2018) 70:90 Page 13 of 13 Scholz CH (1968) The frequency-magnitude relation of microfracturing in rock Utsu T (1965) A method for determining the value of b in a formula log and its relation to earthquake. Bull Seismol Soc Am 58(1):399–415 n = a − bM showing the magnitude–frequency relation for earthquakes Scholz CH (2015) On the stress dependence of the earthquake b value. Geo- (in Japanese). Geophys Bull Hokkaido Univ 13:99–103 phys Res Lett 42(5):1399–1402. https ://doi.org/10.1002/2014G L0628 63 Utsu T (1999) Representation and analysis of the earthquake size distribu- Schorlemmer D, Wiemer S, Wyss M (2005) Variations in earthquake-size distri- tion: a historical review and some new approaches. Pure Appl Geophys bution across different stress regimes. Nature 437(7058):539–542. https :// 155(2–4):509–535 doi.org/10.1038/natur e0409 4 Wessel P, Smith WH, Scharroo R, Luis J, Wobbe F (2013) Generic mapping tools: Schurr B, Asch G, Hainzl S, Bedford J, Hoechner A, Palo M, Wang R, Moreno M, improved version released. EOS Trans Am Geophys Union 94(45):409– Bartsch M, Zhang Y, Oncken O, Tilmann F, Dahm T, Victor P, Barrientos S, 410. https ://doi.org/10.1002/2013E O4500 01 Vilotte JP (2014) Gradual unlocking of plate boundary controlled initia- Wiemer S, Wyss M (2000) Minimum magnitude of completeness in earthquake tion of the 2014 Iquique earthquake. Nature 512(7514):299–302. https :// catalogs: examples from Alaska, the western United States, and Japan. doi.org/10.1038/natur e1368 1 Bull Seismol Soc Am 90(4):859–869 Shi Y, Bolt BA (1982) The standard error of the magnitude-frequency b-value. Woessner J, Wiemer S (2005) Assessing the quality of earthquake catalogues: Bull Seismol Soc Am 72(5):1677–1687 estimating the magnitude of completeness and its uncertainty. Bull Suyehiro S (1966) Difference between aftershocks and foreshocks in the rela- Seismol Soc Am 95(2):684–698. https ://doi.org/10.1785/01200 40007 tionship of magnitude to frequency of occurrence for the great Chilean Yoon CE, O’Reilly O, Bergen KJ, Beroza GC (2015) Earthquake detec- earthquake of 1960. Bull Seismol Soc Am 56(1):185–200 tion through computationally efficient similarity search. Sci Adv Tamaribuchi K, Moriwaki K, Ueno H, Tsukada S (2016) Automatic hypocenter 1(11):e1501057. https ://doi.org/10.1126/sciad v.15010 57 determination for the Seismological Bulletin of Japan using Bayesian Zaliapin I, Ben-Zion Y (2013) Earthquake clusters in southern California I: identi- estimation (in Japanese with English abstract). Quart J Seis 79:1–13 fication and stability. J Geophys Res Solid Earth 118(6):2847–2864. https :// Tormann T, Enescu B, Woessner J, Wiemer S (2015) Randomness of megathrust doi.org/10.1002/jgrb.50179 earthquakes implied by rapid stress recovery after the Japan earthquake. Nat Geosci 8(2):152–158. https ://doi.org/10.1038/ngeo2 343 Tormann T, Wiemer S, Enescu B, Woessner J (2016) Normalized rupture poten- tial for small and large earthquakes along the Pacific Plate off Japan. Geo - phys Res Lett 43(14):7468–7477. https ://doi.org/10.1002/2016G L0693 09
Earth, Planets and Space – Springer Journals
Published: May 29, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.
Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.
All the latest content is available, no embargo periods.
“Hi guys, I cannot tell you how much I love this resource. Incredible. I really believe you've hit the nail on the head with this site in regards to solving the research-purchase issue.”Daniel C.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud
“I must say, @deepdyve is a fabulous solution to the independent researcher's problem of #access to #information.”@deepthiw
“My last article couldn't be possible without the platform @deepdyve that makes journal papers cheaper.”@JoseServera