A New Method for Imaging Seismic Quiescence and Its Application to the Mw = 8.3 Kurile Islands Earthquake on 15 November 2006

In the present paper, a new method referred to as the Poisson probability map (PMAP) method is presented for identifying and visualizing seismic quiescence. With the PMAP, the P-value is defined as the probability that consecutive earthquakes occur according to a homogeneous Poisson process: the smaller the P-value, the less frequently the longer time interval is observed, i.e. the more significant the seismic quiescence. The PMAP method was applied to the sequence which preceded the Kurile Islands earthquake that occurred on 15 November 2006 [Mw = 8.3 and the centroid = (154.33 °E, 46.71 °N)]. The seismic quiescence is identified by a small P-value of 9.0 × 10–5 that was found to start in 1990.1, which lasted for 15.4 years and ended in 2005.5 within a circular area centered at (153.8 °E, 47.1 °N) and with a radius of 26 km. This seismic quiescence has not previously been recognized using any other method.


Introduction
The seismic quiescence hypothesis is an empirical model in which the number of small earthquakes decreases in and around the assumed focal area of a great earthquake when its occurrence time is approaching (e.g., Inouye 1965;Utsu 1968;Mogi 1969;Kelleher and Savino 1975;Ohtake et al. 1977;Ogata 1992). Discussion continues as to whether this hypothesis is correct, and therefore additional case studies are required. Various methods have been proposed to investigate seismic quiescence, including the RTL method (Sobolev and Tyupkin 1997;Huang 2004;Gentilia et al. 2019) and the ZMAP method (Wiemer and Wyss 1994). The ZMAP method is a statistical technique for systematically searching for spatial and temporal changes in seismicity. For a more specific explanation of the ZMAP method, see Wiemer and Wyss (1994) and Katsumata (2011). Since the ZMAP computer program has a rich graphical user interface and is very easy to operate, it has been commonly used to investigate the seismic quiescence that occurs before large earthquakes (e.g., Wyss et al. 1999;Wu and Chiao 2006;Katsumata 2011). In the ZMAP method, a study area is divided into a grid, N earthquakes are selected around each node, and the Z-value is calculated using the following equation: where R bg and R w are the occurrence rates for background earthquakes and those occurring in a time window T w , respectively, S bg and n bg are the variance and the number of bins in the background period, respectively, and S w and n w are the variance and the number of bins in time window T w , respectively. A simple example of the calculation of the Z-value is described in supplementary material 1. Note that large positive Z-values correspond to a decrease in the seismicity rate. A great earthquake occurred on 15 November 2006 in the central Kurile Islands region. The centroid moment tensor (CMT) solution for this earthquake was determined by the Global CMT project (Dziewonski et al. 1981;Ekström et al. 2012; available at https://www.globalcmt.org/CMTsearch. html). The centroid was located at 46. 71°N and 154.33°E at a depth of 13.5 km. The seismic moment was 3.51 9 10 21 Nm (M w = 8.3), and the fault plane was (strike, dip, slip) = (215°, 15°, 92°). This event was a low-angle megathrust interplate earthquake on the upper surface of the Pacific Plate subducting beneath the North American or the Okhotsk Sea Plate (e.g., Seno et al. 1996;Takahashi et al. 1999). This great earthquake filled the southwestern half of a seismic gap (Fedotov 1968;Kelleher and McCann 1976;McCann et al. 1979;Nishenko 1991). This seismic gap is located along the Kurile Trench between the great 1963 Kuril Islands (M w = 8.5) earthquake (Kanamori 1970;Beck and Ruff 1987) and the great 1952 Kamchatka (M w = 9.0) earthquake (Kanamori 1976;Johnson and Satake 1999).
Although Katsumata (2017) investigated the seismicity in and around the Japan subduction zone from 1975 to 2012 in order to search for the longterm seismic quiescence preceding great earthquakes using the ZMAP method, no seismic quiescence was detected before the 2006 Kurile Islands earthquake. The reasons for this oversight appeared to be due to the following issues with the ZMAP method: (1) When no earthquake is observed in a time window T w , i.e. R w = S w = 0, the Z-value is approximately the same whether T w is long or short. (2) The length of the time window T w is arbitrary and must be empirically assumed. (3) The number N of earthquakes around each node is also arbitrary and is the same for all nodes. In order to address these issues, a new method is presented herein. The new method calculates the probability that no earthquake will occur. The longer the period when the earthquake does not occur, the smaller the probability; therefore, it becomes a measure of seismic quiescence. The proposed method is applied to the M w = 8.3 Kurile Islands earthquake on 15 November 2006, and the suitability of the proposed method is discussed.

Data
An earthquake catalog compiled by the International Seismological Center (ISC) (International Seismological Center 2013) was used to investigate the seismicity in and around the focal area of the 2006 Kurile Islands earthquake. The study area was defined by the circle of investigation (CI) used in the M8 algorithm (Keilis-Borok and Kossobokov 1990). The diameter of the CI was calculated using D = exp (M w -5.6) ? 1.0 = 15.88°for a meridian length of 1763 km (Fig. 1). The center of the CI was located at the centroid of the main shock, and earthquakes within a radius of 882 km were downloaded from the Reviewed ISC Bulletin (ftp://isc-mirror.iris. washington.edu/pub/ffb/catalogue/). Although the ISC has finished rebuilding the bulletin from 1964 to 1979 (Storchak et al. 2017), old data stored in the ''prerebuild'' directory (ftp://isc-mirror.iris. washington.edu/pub/prerebuild/ffb/catalogue/) were used in the present study to maintain the temporal homogeneity of the earthquake catalog. Earthquakes satisfying the following three conditions were selected: (1) the earthquake occurred between 1 January 1964 and 14 November 2006, (2) the body wave magnitude was m b C M C , and (3) the depth of the hypocenter was B 60 km. Here, M C is the minimum magnitude of completeness, and the number of earthquakes versus magnitude was examined for the catalog in order to estimate M C . In the present study, M C is defined as the magnitude for which at least 90% of the frequency-magnitude distribution can be modeled using a power law (Wiemer and Wyss 2000). The result indicated that M C ranged from 4.2 to 4.7, depending on the period. Figure 2 shows a temporal change in M C . In general, M C fluctuates in time due to slight changes in observational conditions, such as installation of new seismic stations, closing of old seismic stations, temporary interruption in observation, or changes in the combination of seismic stations used for magnitude calculation. Taking the fluctuation into account, M C was set to be 5.0 in the present study. Consequently, the number of earthquakes with m b C M C was 2239 within the CI ( Fig. 1 and supplementary material 2).
Seismicity typically comprises two kinds of activity: background activity and cluster activity. Aftershocks and earthquake swarms are typical of cluster activity. In the present study, the temporal change in background activity was investigated in order to search for seismic quiescence before the 2006 Kurile Islands earthquake, and thus the cluster activity was removed from the ISC catalog using a stochastic declustering method (Zhuang et al. 2005), which is based on an epidemic-type aftershock sequence (ETAS) model (Ogata 1988(Ogata , 2004. The ETAS parameters in model 5 of Zhuang et al. (2005) were determined as follows: A = 0.330, a = 1.58, c = 0.0177, p = 1.217, D 2 = 0.00827, q = 1.96, and c = 1.21. After the declustering process, 1096 earthquakes were identified as background seismicity and are the basis of the analyses in the present study ( Fig. 1 and supplementary material 3). Although M C fluctuates slightly in time, a cumulative number curve in Fig. 2 shows that the background seismicity rate is almost constant for earthquakes with M = 5.0 and larger.

Method
In order to search for a seismic quiescence, a rectangular area was defined to include the CI, i.e. 40-54°N and 143-163°E. This area was divided into a grid with an interval of 0.1°in latitude by 0.1°in longitude. A circle centered on a node (x, y) was drawn with a radius of R [km], and R was increased Earthquakes considered in the present study. a and c Earthquake distributions before and after the declustering process, respectively. Each earthquake list is provided in supplementary material 2 and 3. Earthquakes were selected from the ISC catalog and occurred from 1 January 1964 to 14 November 2006, with m b C 5.0 and depth B 60 km. The large circle is the circle of investigation (CI), and the solid line along the Kurile Islands indicates the plate boundary (Bird 2003). b and d Space-time plots before and after the declustering process, respectively gradually so that the number of earthquakes was N within the circle. This circle is defined as the resolution circle, and the radius R is defined as the spatial resolution (Wiemer and Wyss 1994).
In order to resolve the issues with the ZMAP method, another parameter was introduced in the present study. First, we assume that earthquake occurrence follows a stationary Poisson process in time (e.g., Utsu 1969;Gardner and Knopoff 1974), and the average rate of occurrence l is defined as: where T is the observation time in years, and N is the number of earthquakes that occurred in the resolution circle during T. In this case, an average of k earthquakes occur during dt years: Assuming that an average of k earthquakes occur, the probability that n earthquakes occur during dt years is: Therefore, the probability that no earthquake occurs during dt years, i.e. dt is the interval between earthquake occurrence, is obtained by substituting 0 for n in Eq. (4): In the present study, P(0) in Eq. (5) is used to measure the seismic quiescence quantitatively. As mentioned above, there are N earthquakes within the resolution circle. Suppose that their origin times are t 1 , t 2 , …, t i , t i?1 , …, t N in order of occurrence. Here, dt is defined at time t as follows: In this case, P(0) is obtained as follows: where P N (0) is a function of N ranging from 5 to 40 in the present study. Finally, the P-value is defined at node (x, y) at time t as: where x ranges from 143°E to 163°E with an interval of 0.1°, y ranges from 40°N to 54°N with an interval of 0.1°, and t ranges from 1964.8 to 2006.8 with an interval of 0.1 years. If the spatial resolution R is larger than 50 km for all N values at a node, then the P-value was not calculated at the node. The smaller the P-value, the more significant the seismic quiescence. A simple example of the P-value calculation is as follows. Suppose that this is an earthquake catalog from 1964.0 to 2007.0, i.e. 43.0 years, and a time slice is produced at time t = 2000.0. First, N = 15 earthquakes are selected around a spatial node, and they are arranged from oldest earthquakes in chronological order (Fig. 3a). Let the occurrence times of these 15 earthquakes be t 1 , t 2 , …, t 15 . The value of i that satisfies the condition t i \ t \ t i?1 is i = 12, i.e. t 12 = 1990.137 and t 13 = 2005.581. In this example, dt = t t 12 = 2000.0 -1990.137 = 9.863 years. The average rate of occurrence is l = 15/43.0 = 0.3488 events/year, k = ldt = 0.3488 9 9.863 = 3.440, and thus P(0) = exp(-k) = exp(-3.440) = 0.0321. Second, N = 28 earthquakes are selected around the node (Fig. 3b). The value of i that satisfies the condition As a result, the P(0) becomes the minimum when N = 28. At this spatial node, the minimum value, i.e. P(0) = 0.00162, is defined as the P-value at the time t = 2000.0. In actual analysis, N increases from 5 to 40 by 1 for finding the smallest P(0) at each node. Moreover, the time t increases from 1964.8 to 2006.8 by 0.1 years. Using the above method, we were able to create a P-value spatial map every 0.1 year.

Results
As in the ZMAP method, the temporal change in the P-value was visualized. One frame of time was produced and displayed every 0.1 year, and the frames were used to visualize small P-values ( Fig. 4  and supplementary material 4). Note that the smaller the P-value, the more significant the seismic quiescence. These visual time slices are called a Poisson probability map (PMAP) in the present study. By investigating all frames of time using PMAP, a node with a very small P-value was found in the focal area of the 2006 Kurile Islands earthquake. The parameters are as follows: x = 153.8°E, y = 47.1°N, t = 2005.5, N = 26, R = 26 km, dt = 15.4 years, and P = 9.0 9 10 -5 or log 10 P = -4.05. The P-value at this node is the smallest in and around the focal area of the 2006 earthquake between t = 1964.8 and 2006.8. The second smallest P-value was obtained for the following parameters: x = 153.5°E, y = 47.0°N,   Wyss and Habermann (1979) pointed out the seismic quiescence corresponding to the second smallest P-value at t = 1976.5. Interestingly, the smallest and second smallest P-values were located very close to each other in space. However, they were separated by 29 years. This strongly suggests that not all seismic quiescence leads to a subsequent great earthquake. The seismic quiescence at t = 2005.5 was investigated in detail (Fig. 5a). The node is located 59 km northwest of the centroid of the 2006 main shock. In the present study, the seismic quiescence area that suffered from a significant decrease in seismicity is defined as a circle centered at (x, y) with a radius of R. The seismic quiescence area is much smaller than the focal area with a slip larger than 1 m during the main shock and is located close to the area with a slip larger than 10 m. This can be considered as an asperity with a large slip. If this is the case, the seismic quiescence area is strongly associated with the asperity that was ruptured by the subsequent main shock. The temporal features of the seismic quiescence found above are as follows (Fig. 5b). Twentythree earthquakes occurred from 1964.0 to 1990.1, and thus the average occurrence rate was 23/26.1 = 0.88 events/year. Although earthquakes occurred rather regularly during this period, since 1990.1, earthquakes have completely ceased to occur. The seismic quiescence started at t = 1990.1, and the duration was dt = 15.4 years (from 1990.1 to 2005.5  .1, and after 1990.1, monotonically decreased to below -4 (Fig. 5c). The log 10 P-value appears to fluctuate between -2 and 0 when earthquakes occur randomly without any long-term seismic quiescence (Fig. 5d). Compared to usual fluctuations of the P-value, the value of log 10 P = -4.05 is two orders of magnitude smaller.

Discussion
In the present study, a long-term seismic quiescence was found to be located very close to an asperity that was ruptured by the subsequent main shock of the 2006 Kurile Islands earthquake. As a measure of the rarity of the seismic quiescence, the P-value was introduced. The smaller the P-value, the rarer the occurrence of seismic quiescence. In this case, P = 9.0 9 10 -5 or log 10 P = -4.05 is the smallest value in and around the focal area. Although Katsumata (2017) failed to detect seismic quiescence before the 2006 Kurile Islands earthquake using the Year (d) Figure 5 a The cross indicates the spatial node at the center of the seismic quiescence area. Red circles are epicenters selected around the node and are located within the seismic quiescence area. The black star indicates the epicenter of the main shock. The coseismic slip distribution model for the 2006 Kurile Islands earthquake is shown as contours of 1, 4, 7, 10, and 13 m (Lay et al. 2009). The plate boundary is shown by solid lines (Bird 2003). b Space-time plot (red circles) and cumulative number curve (black line) of the epicenters in (a). The time period for the seismic quiescence is shown in gray. c Change in P-value as a function of time at the node shown as the cross in (a). The time period of seismic quiescence is shown in gray. d Temporal change in P-value without seismic quiescence. Red and blue lines indicate temporal changes in Pvalues at nodes A and B in (a), respectively ZMAP method, the seismic quiescence 16.7 years before the 2006 earthquake was successfully detected using the PMAP method in the present study. The next problem to consider is how frequently such small P-values are observed. The P-value was obtained at 3657 spatial nodes within the study area and at 421 temporal nodes between 1964.8 and 2006.8. Thus, the total number of P-values calculated is 3657 9 421 = 1,539,597. Note that the P-values are not independent from each other. Since the distance between neighboring nodes are 0.1°, their resolution circles are partly overlapped, and the same earthquakes are sampled. The P-value ranged from 1.6 9 10 -5 to 1.0, i.e. -4.8 B log 10 P B 0, in the entire study area (Fig. 6a). The number of P-values less than or equal to 1.0 9 10 -4 is 170, which is 0.01% of the total number. The spatial distributions of the nodes are shown for log 10 P B -2.0, -3.0, and -4.0 in Fig. 6b, c, and d, respectively. These nodes appear to be clustered rather than scattered uniformly along the Kurile Trench. When log 10 P B -4.0, the 170 nodes are clearly divided into four groups. The Closed stars indicate earthquakes that occurred within the CI with M w C 7.5 with a centroid at a depth B 100 km listed on the Global CMT catalog from 1976 to 2006. The P-value was calculated at nodes in the gray area. The plate boundary is shown by solid lines (Bird 2003).
c and d Distribution of nodes with log 10 P B -3.0 and log 10 P B -4.0, respectively first group is located near the 2006 earthquake, as described above. The second group is located around a node (44.2°N, 149.4°E) 14 km east of the centroid of the 1978 M w 7.6 and M w 7.5 earthquakes (Fig. 7a,  b). The third group is located at a node (  Red circles indicate earthquakes that occurred within the seismic quiescence area, and the center of the area is shown by a cross. The plate boundary is shown by solid lines (Bird 2003). b, d, and f Space-time plot of epicenters within the seismic quiescence area in red circles, and the cumulative number curve is shown by black lines. Vertical red lines indicate the time of occurrence of the main shock M w 7.6 earthquake (Fig. 7c, d). The fourth group is located around a node (43.5°N, 146.7°E) 76 km west of the centroid of the 1994 M w 8.3 earthquake (Fig. 7e, f). The seismic quiescence usually ends before the main shock; however, in Fig. 7f, it ends after the main shock. This is a natural result of removing aftershocks by a declustering process. Six main shocks within the CI are listed in the Global CMT catalog between 1976 and 2006 with M w = 7.5 or larger and a centroid depth = 100 km or shallower. Good correlation is observed between the seismic quiescence and the six large earthquakes in the study area. However, there is still insufficient evidence of a causal relationship. If causality is not proven, it cannot be concluded that seismic quiescence is a precursor to large earthquakes. We estimated the probability that the seismic quiescence is observed by chance through conducting a numerical simulation. The numerical simulation is based on a method developed by Zhuang et al. (2004). We produced a simulated earthquake catalog including the background and clustered seismicity in the Kurile Islands region based on the ETAS model with an assumption of the ETAS parameters obtained in this study. An R code used to produce the simulated catalog is presented in supplementary material 5. We applied the declustering process to the simulated catalog and calculated the P-value in the same way as when analyzing a real earthquake catalog. The smallest P-value was searched in the focal area, i.e. 45.5-48.0°N and 152-156°E and in the time period between t = 2005.0 and 2006.8. The procedure described above was repeated 500 times and made a histogram of the smallest P-value (Fig. 8). In the analysis of the actual earthquake catalog, the smallest P-value at t = 2005.5 was log 10 P = -4.05. In the numerical simulation, the smallest P-value was less than log 10 P = -4.05 in 20 out of 500 times, i.e. 4%. Therefore, we concluded that the probability of this seismic quiescence being observed by chance is quite low.

Conclusions
Katsumata (2017) pointed out that no long-term seismic quiescence was observed before the M w = 8.3 Kurile Islands earthquake on 15 November 2006 using the ZMAP method. In the present study, a new method, PMAP, was developed in order to improve on the performance of the ZMAP method, and it was applied to the 2006 earthquake. The seismic quiescence was found to have started at the beginning of 1990, i.e. 16.7 years before the occurrence of the 2006 earthquake, and the seismic quiescence area was located very close to the large asperity ruptured by the 2006 main shock. Thus, seismic quiescence that was not detected by the ZMAP method was detected by the PMAP method, indicating that PMAP is more effective, at least for the 2006 earthquake. supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.