Nanjo et al. Earth, Planets and Space (2016) 68:187
DOI 10.1186/s40623-016-0558-2

Open Access

FULL PAPER

Seismicity prior to the 2016 Kumamoto
earthquakes
Kazuyoshi Z. Nanjo1 , Jun Izutsu2 , Yoshiaki Orihara3, Nobuhiro Furuse4, Shoho Togo3, Hidetoshi Nitta3,
Tomohiro Okada3, Rika Tanaka3, Masashi Kamogawa3 and Toshiyasu Nagao4*

Abstract
We report precursory seismic patterns prior to the 2016 Kumamoto earthquakes, as measured by four different methods based on changes in seismicity that can be used for earthquake forecasting: the b-value method, two methods ofseismic quiescence evaluation, and an analysis of seismicity density in space and time. The spatial extent of precursory patterns differs from one method to the other and ranges from local scales (typically, asperity size) to regional
scales (e.g., 2° × 3° around the source zone). The earthquakes were preceded by periods of pronounced anomalies,
which lasted in yearly scales (1.5 years), or longer (>3 years). We demonstrate that a combination of multiple methods
detected different signals prior to the Kumamoto earthquakes. This indicates great potential to reduce the hazard at
possible future sites of earthquakes relative to long-term seismic hazard assessment. We also found that the seismic
quiescence in a regional-scale area, detected by using the two methods of seismic quiescence evaluation, was a common precursor to the 2016 Kumamoto earthquakes and 2015 Off Satsuma Peninsula earthquake. The result allows
us to interpret both events as the onset that occurred at a section along the tectonic line from the Okinawa Troughthrough the Beppu–Shimabara graben.
Keywords: Crustal deformation, Earthquake, Physical property of rocks, Time-series analysis, Probabilistic forecasting,Computational seismology, Statistical seismology, 2016 Kumamoto earthquakes, Critical phenomenon
Introduction
The 2016 Kumamoto earthquakes, including a magnitude M7.3 mainshock that struck on April 16, 2016, as
well as active foreshocks and aftershocks, fulfilled a section of the Futagawa–Hinagu fault zones. The zones areencompassed by the Beppu–Shimabara graben, a geological formation that runs across the middle of Kyushu,
from Beppu Bay in the east to the Shimabara Peninsula in
the west. This is understood as being the result of crustal deformation caused by the rifting and spreading of the
Okinawa Trough, which is viewed as a continuation of
the Beppu–Shimabara graben (Tada 1985). The November 14, 2015, M7.1 earthquake occurred Off Satsuma
Peninsula in the Okinawa Trough. This event and the*Correspondence: nagao@scc.u‑tokai.ac.jp
4Institute of Oceanic Research and Development, Tokai University, 3‑20‑1,
Orido, Shimizu, Shizuoka 424‑0902, Japan
Full list of author information is available at the end of the article

Kumamoto earthquakes are considered to be the onsetthat occurred under the same tectonics (Ishibashi 2016).
A wide variety of approaches have been applied toearthquake prediction and forecasting. Most proposed
prediction and forecasting methods rely on the concept of a diagnostic precursor, i.e., some kind of signal
observable before earthquakes that indicates with highprobability the location, time, and magnitude of animpending event. Precursor methods include changes
in strain rates, seismic wave speeds, and electrical conductivity; variations of radon concentrations in groundwater, soil, and air; fluctuations in groundwater levels;
electromagnetic variations near and above the Earth’s
surface; and seismicity patterns. Despite this, the searchfor diagnostic precursors has not yet produced a successful short-term prediction scheme (e.g., Keilis-Borok
2002; Scholz 2002; Kanamori 2003; International Commission on Earthquake Forecasting for Civil Protection2011).

© The Author(s) 2016. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/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.

Nanjo et al. Earth, Planets and Space (2016) 68:187

Despite such a notable lack of success, there has been aresurgence of research on earthquake predictability motivated by better monitoring networks and data on past
events, new knowledge of the physics of earthquake ruptures, and a more comprehensive understanding of stressevolution and transfer. The California and international
working groups, RELM (Field 2007) and CSEP (Jordan
2006), have been supporting research on earthquake predictability, conducting scientific experiments under rigorous, controlled conditions and evaluating them using
accepted criteria specified in advance. These groups point
to future directions of model development of earthquake
prediction and forecasting as well as its testing. However,making full use of currently available resources for and a
new knowledge and understanding of seismology, there
needs predictability research that is not directly associated with the RELM and CSEP, but rather that addresses
more general questions about the improvement of different hypotheses on earthquake generation and different
concepts related to diagnostic precursor.
Different methods exist to measure, map, and evaluate possible episodes before mainshocks. Examples
include the b-value method (e.g., Schorlemmer and Wiemer 2005; Nanjo et al. 2012), RTL-/RTM-algorithms(e.g., Sobolev and Tyupkin 1997; Nagao et al. 2011), theZ-value method (e.g., Wiemer and Wyss 1994; Wyss et al.
2004; Katsumata 2015), and seismicity density analysis (Lippiello et al. 2012). To our knowledge, there is no
standardized approach that encourages researchers torely exclusively on a single method. Instead, multiple
methods must be used to find evidence of precursory episodes. One disadvantage is that the results obtained bydifferent methods may not be easily compared, although
one may gain more confidence in impending earthquakes when using different methods. Also, additional
insight may be gained because of intrinsic differences inthe statistical characterization of seismic patterns. Earlyattempts were made by Wyss et al. (2004) who used the
RTL-algorithm and Z-value methods and by Enescu andIto (2001) who used the b-value and Z-value methods,
and a fractal dimension approach.
This paper reports the first results of recognizing seismic patterns as possible precursory episodes to the 2016
Kumamoto earthquakes using existing four different
methods: the b-value method, two kinds of seismic quiescence evaluation methods (RTM-algorithm and Z-value),and seismicity density analysis. These methods are basedon different assumptions, selection of sampling volumes,
algorithms, and definitions of anomalies; thus, they are
representative of a wide variety of methods that can beused to detect precursory anomalies. Before presenting
further details of our study, we present a brief overview
of the methods that were used.

Page 2 of 10

Methods
The b‑value method

The b-value is the slope of the Gutenberg–Richter (GR)frequency–magnitude distribution of earthquakes(Gutenberg and Richter 1944), log10N = a − bM, where
N is the cumulative number of earthquakes of M or
greater, a is the earthquake productivity of a volume, andb is used to describe the relative occurrence of large andsmall events (i.e., a high b-value indicates a larger proportion of small earthquakes, and vice versa). Spatial andtemporal changes in b are known to reflect a structural
heterogeneous structure (e.g., Mogi 1962) such as strong
coupling along subduction zones and a magma chamber (e.g., Tormann et al. 2015). In the laboratory and theEarth’s crust, the b-value is also known to be inversely
dependent on differential stresses (Scholz 1968, 2015). In
this context, measurements of spatial temporal changes
in b could act as a stress meter to help image asperities,
the highly stressed patches in faults where future ruptures are likely to occur (e.g., Schorlemmer and Wiemer
2005; Nanjo et al. 2012; Tormann et al. 2015).
The RTM‑algorithm

This is a type of weighted method that assesses the time,
space, and the size of earthquakes, and is a modified RTLalgorithm (e.g., Sobolev and Tyupkin 1997, 1999; Huang
and Nagao 2002; Huang 2004, 2006; Sobolev 2011). The
RTM-algorithm measures the level of seismic activity inmoving windows by counting the number of earthquakesthat are weighted by their size and inversely weighted bytheir distances in time and space from the point of observation. A detailed description is provided by Nagao et al.
(2011). A decrease in the RTM value implies a decrease
in seismicity compared to the background rate aroundthe investigated location (a seismic quiescence). A recovery stage from quiescence to the background level canbe broadly considered as foreshock activation. The RTM
method evaluates both seismic quiescence and the following stage of activation. In addition, the location of
extreme, if it is an anomaly, can be found by performing
RTM-calculations with the centers of the sampling circles
at the nodes of a grid.
The Z‑value method

The Z-mapping approach measures the difference inseismicity rate, within moving time windows, to thebackground rate by a standard deviation, Z (Wiemer
and Wyss 1994; Wyss et al. 2004; Katsumata 2015). The
purpose is to detect possible periods of anomalously low
seismicity just before the mainshock near its epicenterand to evaluate the statistical significance of such quiescence compared with all other parameters that may haveoccurred at random times and locations. The standard

Nanjo et al. Earth, Planets and Space (2016) 68:187

deviation Z is defined as Z = (R1 − R2)/(S1/n1 + S2/n2)1/2,
where R1 and R2 are the mean rates, S1 and S2 are thevariances, and n1 and n2 are the number of events in thefirst and second periods to be compared. The larger theZ-value, the more significant the observed difference.
Seismicity density analysis

An increase in the number of smaller magnitude events
before large earthquakes is often observed. The linear
density probability of earthquakes occurring before andafter events defined as mainshocks, ρ(Δr), as a function of
the distance from the mainshock hypocenter, Δr, displays
a symmetric behavior (Lippiello et al. 2012). This behavior indicates that the size of the area fractured during themainshock is associated with the spatial organizationof seismicity before the mainshock, because the spatialorganization after the mainshock reflects the mainshock
rupture area (Lippiello et al. 2012). The analysis based on
Lippiello et al. (2012), hereinafter referred to as the seismicity density analysis, can be implemented to forecast
the size of the mainshock. The objective of this paper isto detect such symmetric behavior for the Kumamoto
earthquakes and to evaluate the possible size of the preparatory areas to the earthquakes.

Data
Our dataset is the earthquake catalog maintained by theJapan Meteorological Agency (JMA). We mainly usedearthquakes since 2000, the start of modernization of theseismic networks feeding data to JMA (Obara et al. 2005).
Since aftershocks add noise to both the Z-value method
and the RTM-algorithm, we used a declustered catalog
as a basic input for these methods. We eliminated aftershocks using the JMA standard declustering program(for details, see Appendix 1 and Fig. A1 in Additional
file 1), which is classified as a link algorithm (Frolich andDavis 1990).
A reliable estimate of completeness magnitude Mc,
above which all earthquakes are considered to be detectedby a seismic network, is vital for seismicity-related studies. We paid attention to catalog completeness and performed a completeness check as a preprocessing step ofindividual analyses while referring to a comprehensive
analysis of Mc in Japan (Nanjo et al. 2010).
Results
The b‑value method

A map view (Fig. 1a) of b-values based on seismicity overa period from 2000 to April 14, 2016, 21:25 (immediately before the Kumamoto earthquakes) reveals a zoneof low b-values near the eventual epicenters of the mainshock and two M6.5-class foreshocks. The characteristic
dimension of unusually low b-values in an approximately

Page 3 of 10

0.3° × 0.3° area in Fig. 1a is ≈30 km. The M5 earthquake
on June 8, 2000, occurred near the mainshock hypocenterso that we excluded seismicity related to this event andcreated maps of b-values (Fig. A2 in Additional file 1).
These maps show that a zone of low b-values is insensitive to this exclusion, indicating that the 2000 M5 earthquake is not the main cause of the low spatially mappedb-values in Fig. 1. Comparison with the mainshock rupture zone (National Research Institute for Earth Scienceand Disaster Resilience 2016) shows the influence ofstructural heterogeneity on spatial distribution in b-values in such a way that rupture propagation of the mainshock terminated an area near Mt. Aso with b ~ 1, higher
than b-values within the rupture zone (b = 0.6–0.8).
The high-temperature area around the magma chamberof Mt. Aso may have contributed to termination of the
rupture (Yagi et al. 2016) and high b-values (e.g., Roberts et al. 2015). We compared the map in Fig. 1a withthe sequence of foreshocks, mainshock, and aftershocks
with M5 or greater. Areas of low b-values contained 84%of the sequence (Fig. A3 in Additional file 1). Two of thethree aftershock regions, which lay further to the northeast, were observed around the high b-value zones near
Mt. Aso and Mt. Yufu, as observed in previous studies
about b-value characterization of volcanic systems (e.g.,
Roberts et al. 2015).
For earthquakes in the cylindrical volume with a radius
R = 10 km, centered at the epicenter of the mainshock(Fig. 1c), b-values after the fluctuation due to aftershocks of the 2000 M5 earthquake were stable (see also
Fig. A4 in Additional file 1). The precursory duration of
the b-value is unknown, because the time variation in b
shows no clear change point during the analyzed period.
A similar time-series analysis was performed for seismicity in the Futagawa, Hinagu, Aso-Kuju, and Beppu fault
zones (see Figs. A5–A9 in Additional file 1). The b-values show no systematic increase or decrease for all fault
zones. A summary of the results is provided in Table 1.
RTM‑algorithm

As described in “Data” section, we first eliminated aftershocks. Next, we conducted a completeness analysis of
the JMA catalog for the study region (Fig. 2) and obtained
typically Mc = 2 so that we primarily set a lower threshold at Mmin = 2. Then, through an extensive parameter
survey, we found anomalous seismic quiescence patterns
as explained next.
A temporal variation of RTM at an almost equidistantpoint from the 2015 Off Satsuma Peninsula earthquake
and the 2016 Kumamoto mainshock showed a significant decrease prior to the occurrence of both (Fig. 2a, seealso Table A1 in Additional file 1). A deviation from thebackground level started at 2014.8 years for mid-OctoberNanjo et al. Earth, Planets and Space (2016) 68:187

Page 4 of 10

Fig. 1 a Map of b-values obtained from events from January 1, 2000, to April 14, 2016, 21:25 (immediately before the Kumamoto earthquakes).
Overlapped on the map are the hypocenter (yellow star) and rupture zone (rectangle) of the M7.3 mainshock; M6.5-class foreshocks (red stars); M ≥ 5foreshocks (squares); and M ≥ 5 aftershocks (circles). Triangles: volcanos; red line: faults. We calculated b-values (Woessner and Wiemer 2005) for
events with depths 25 km or shallower falling in a cylindrical volume with radius R = 5 km, centered at each node on a 0.02º × 0.02º grid (we did
not use a fixed number of earthquakes) and plotted a b-value at the corresponding node only if at least 50 events in the cylinder yielded a good fit
to the GR law. b GR fitting is shown for points A, B, and C in a. c Plot of b-values as a function of time, as obtained from seismicity in the cylindricalvolume with R = 10 km, centered at the mainshock. We used a moving window approach, whereby the window covered 500 events. Uncertainty
estimates are according to Shi and Bolt (1982). Gray squares and open circles show the b-values obtained from foreshocks and aftershocks, respectively. Horizontal line regional average b = 0.93; vertical line moment of the mainshock

Table 1 Characteristic precursory phenomena
MethodSpatial extent

b-value method

Area of 0.3° × 0.3° containing a part of the source zone of the Unknown
Kumamoto mainshock

RTM-algorithm

Area of 2° × 3° between the Kumamoto and Off Satsuma
Peninsula mainshocks

1.5 years for the Kumamoto mainshock (1 year for the Off
Satsuma Peninsula mainshock)Z-value method

Area with R = 100 km between the Kumamoto and Off
Satsuma Peninsula mainshocks

1.5 years for the Kumamoto mainshock (1 year for the Off
Satsuma Peninsula mainshock)Seismicity density analysis

Area with Δr = 8–40 km containing a part of the source zone 3 years or longer for the Kumamoto mainshockof the Kumamoto mainshock

2014, and the strongest deviation in 2015.8 years for midOctober 2015 was about −20. During the critical period,
the R, T, and M functions attained values of about −2.5
to −3.0. During a recovery stage from the quiescenceto the background level, the 2015 Off Satsuma Peninsula earthquake and the 2016 Kumamoto mainshock
occurred. This property is similar to that documented byother studies that used the RTL-/RTM-algorithm (e.g.,

Duration

Sobolev and Tyupkin 1997, 1999; Huang and Nagao 2002;
Huang 2004, 2006; Nagao et al. 2011; Sobolev 2011).
Since the RTM-algorithm is statistical and nonlinear, wewere unable to identify which of earthquakes contributed
to the recovery stage; this topic is beyond the scope of
our work in this paper.
A map view of RTM values on October 1, 2015, reveals
a highly significant change between the 2015 Off Satsuma

Nanjo et al. Earth, Planets and Space (2016) 68:187

Page 5 of 10

Fig. 2 a Temporal variations of the RTM (black), R (blue), T (green), and M (red) at the point indicated by a cross in b. Red arrows show the momentsof the 2015 Off Satsuma Peninsula earthquake and the mainshock of the 2016 Kumamoto earthquakes; both epicenters are indicated by stars in b.
b Areas of seismic quiescence before the two earthquakes (stars). The figure, generated by a grid with 0.1° × 0.1° spacing, shows the most quiescent period on October 1, 2015, indicated by a black arrow in a (45 days before the 2015 Off Satsuma Peninsula earthquake). The cross is the pointof time variation shown in a. Parameters used were R0 = 100 km, Rmax = 200 km, T0 = 1 year, Tmax = 2 years, and Mmin = 2 (Table A1 in Additional
file 1). For parameter definition, see Nagao et al. (2011)

Peninsula earthquake and the 2016 M7.3 Kumamoto
mainshock indicated by stars in Fig. 2b, in which theparameter set was the same as that for Fig. 2a. The characteristic dimension of seismic quiescence in an approximately 2° × 3° area is 200–300 km. We created a mapof RTM values using the same parameter set as that for
Fig. 2b except that a higher threshold (Mmin = 2.4) wasapplied. The position of the anomalous zone is not very
sensitive (see panel E of Fig. A10 in Additional file 1). We
also used three different sets of parameters with Mmin = 2
and created RTM maps to support a more robust anomalous pattern (see panels A–C of Fig. A10 and Table A1 in
Additional file 1).
To strengthen the related conclusion, it would be betterto conduct a reliability analysis of the revealed anomalies in detail, e.g., following the approach taken by Huang(2006). However, the main purpose of this paper was toprovide the first results on how to recognize possible precursory episodes to the 2016 Kumamoto earthquakes, as
described in “Introduction.” Thus, we presented essential parts of the reliability analysis such as completenesscheck, quality test of a declustered catalog, and somelevels of a parameter survey. For a full justification of
the present conclusion, further detailed investigations
involving the robustness of the temporal and spatial pattern should be conducted.

Our results do not indicate that the observed seismic
quiescence precedes either of the mainshocks. The 2015
Off Satsuma Peninsula earthquake and the 2016 Kumamoto earthquakes are considered to be the onset that
occurred at a section along the tectonic line from the
Okinawa Trough through the Beppu–Shimabara graben (Ishibashi 2016). Our results, combined with these
tectonics, indicate that a series of earthquakes along thetectonic line was preceded by seismic quiescence. These
results are summarized in Table 1.
The Z‑value method

Similar to the RTM-algorithm, the processing steps
were performed to ensure data quality. An extensive
parameter survey was then conducted to choose themost suitable parameter set. The cumulative number ofearthquakes with Mmin = 2.4 as a function of time for acircle with sampling radius R = 100 km, centered at apoint about 60 km away to the northeast from the epicenter of the 2015 Off Satsuma Peninsula mainshock,showed an anomaly of a few earthquakes during about
1 year before this mainshock (Fig. 3) and about 1.5 years
before the Kumamoto mainshock. With a look-ahead
time window Tw = 1 year, sampling radius R = 100 km,
and Mmin = 2.4, comparison of the seismicity rate during
the last year before the Off Satsuma Peninsula mainshockNanjo et al. Earth, Planets and Space (2016) 68:187

Page 6 of 10

Fig. 3 Cumulative number of earthquakes (blue, scale on left) with M ≥ 2.5 as a function of time, up to the occurrence of the Off Satsuma Peninsula
earthquake in 2015.8 (star). Curve in red shows the Z-values (scale on right), resulting from every position of the time window Tw = 1 year. The
sample is from a circle with R = 100 km, centered at the black dot in the inset. Also included in the inset for reference are stars showing the 2015 Off
Satsuma Peninsula and the 2016 Kumamoto mainshocks and gray dots showing the centers of the collecting volumes with R = 100 km to create Fig.
A11 in Additional file 1

with the background seismicity rate resulted in Z = 2.9.
This representative result shows a typical precursory quiescence pattern to a large earthquake. At other points,
typical precursory pattern showing fluctuation in Z at
low levels followed by a rapid increase to high values isdegraded or unobservable, as shown in Fig. A11 in Additional file 1 (see a further parameter survey in Appendix
2 in Additional file 1 that supports our results). Figure 3
also shows that during the fluctuation in Z, Z is below
zero at most times, indicating preceding seismic activation to quiescence. Following the RTM-algorithm, weused a lower threshold at Mmin = 2.0. Quiescence patterns were observed, but this characteristic was degraded(Fig. A12 in Additional file 1). In summary, in the areawith R = 100 km centered at a point about 60 km away
to the northeast from the epicenter of the Off Satsuma
Peninsula earthquake, the Z-value analysis detected atypical precursory quiescence, which lasted about 1 year
(Table 1).
Seismicity density analysis

Before assuming the approach based on Lippiello et al.
(2012) for Japan, a complete check of the JMA catalog
since 1998 was performed as a preprocessing step. We
decided to use inland earthquakes only and set a lower
threshold at Mmin = 1. The previous Z-value and RTM

approaches were applied to data including both inlandand offshore earthquakes (Figs. 2, 3). Mc of the JMA catalog was generally lower in inland regions than in offshore
regions (Nanjo et al. 2010), which explains why Mmin = 1
for this seismicity density analysis is different from
Mmin = 2 and 2.4 for the Z-value and RTM approaches.
As shown in the inset in Fig. 4 for mainshock magnitudes 4 ≤ m < 5, ρ(Δr) before and after mainshocks wasvery similar in the whole spatial range for all magnitude
ranges 3 ≤ m < 6, consistent with Lippiello et al. (2012).
The decay of seismicity density with distance is well modeled by an inverse power law ρ(Δr) ~ Δr−η at Δr ≫ 0,
where η is a constant of 1–2 (Richards-Dinger et al.
2010). We assumed a typical η = 1.35 (Felzer and Brodsky 2006).
We next compared the size of spatial organization
before and after mainshocks with the size of asperity (Fig. 4). As a representative of the former size, we
detected the deviation point from a scaling relation
ρ(Δr) ~ Δr−η (the inset in Fig. 4). This deviation point isdefined as the characteristic distance Δrc, below which
the scaling is no longer valid, due to large variance ofseismicity density or low seismicity density at very shortdistances to the mainshock hypocenter, above a location
uncertainty of ~100 m in the JMA catalog. We created aplot of Δrc versus m for 3 ≤ m < 6 in Fig. 4. The latter sizeNanjo et al. Earth, Planets and Space (2016) 68:187

Fig. 4 Plot of the characteristic distance Δrc as a function of mainshock magnitude m for seismicity prior (circles) and subsequent (diamonds) to mainshocks in Japan. The Kumamoto earthquakes (see also
Fig. 5) are shown as a set of data of prior and subsequent seismicity as
follows: blue line segment from m6.5 to 7.3 with a blue circle, indicating seismicity during a 3-year period prior to the M6.5 shock; blue line
segment from m6.5 to 7.3 with a blue diamond, indicating seismicity
subsequent to the M7.3 shock. As a reference, data for seismicity inthe period subsequent to the M6.5 shock until the M7.3 shock areincluded (red diamond). Thin solid line is drawn by using extrapolation from the scaling (thick solid line) of the characteristic asperity
radius la with m of 6.6 to 9.2, based on Somerville et al. (2015). Thindashed lines are drawn by using extrapolation from the one-standarddeviation limits (thick dashed line). The inset shows a plot of ρ(Δr) as afunction of Δr from the mainshock with m = 4 ~ 5 for seismicity prior(circles) and subsequent (diamonds) to mainshocks. Data are fittedwith ρ(Δr) ~ Δr−η with η = 1.35. The deviation point, below which an
inverse power law is no longer valid, is marked by Δrc

was based on the scaling between asperity area Sa andm (Somerville et al. 2015). Assuming a circular asperity
Sa = πl2a where la is the characteristic asperity radius, we
included the la − m relation in Fig. 4. Data (gray points)
show a positive correlation with large scatter, and theΔrc − m correlation appears to be similar to the la − mrelation. To support this idea, data from the Δrc − m correlation need to be increased.
We then considered the Kumamoto case in Fig. 5.
When using the same approach as Lippiello et al. (2012),
the M7.3 shock and a preceding M6.5 shock were close
to each other in space and time of the seismicity density
analysis for m = 6–7. Therefore, we used seismicity subsequent to the M7.3 quake and seismicity prior to theM6.5 quake. Although Fig. 5 has a slightly large scatter,data (blue stars) of seismicity subsequent to the M7.3
quake are approximated by ρ(Δr) ~ Δr−η with η = 1.35
for distances Δr ≥ 8 km (Δrc = 8 km). As a reference,
we applied the same procedure to seismicity subsequent
to the M6.5 shock until M7.3 shock (red triangles) andPage 7 of 10

Fig. 5 Plot of ρ(Δr) as a function of Δr for the Kumamoto earthquakes. Blue stars: seismicity subsequent to the M7.3 shock. As a reference, data on seismicity subsequent to the M6.5 shock until the M7.3
shock are marked by red upward-pointing triangles. Also includedin this figure is an inverse power law with an exponent η = 1.35.
The characteristic distance Δrc, below which an inverse power law is
no longer valid, is marked for respective data. The inset shows a plot
for seismicity prior to the M6.5 shock for different periods: 18 years(crosses), 12 years (gray downward-pointing triangles), 6 years (yellow
diamonds), and 3 years (blue circles). Data for the 3-year period are
fitted with an inverse power law with η = 1.35. The deviation point,
below which an inverse power law is no longer valid, is marked byΔrc. Gray shading indicates distances larger than 55 km, beyond whichdata are referred to as background density levels

ρ(Δr) ~ Δr−η with η = 1.35 being a reasonable approximation for Δr ≥ 4 km (Δrc = 4 km). Analysis of seismicity 3,
6, 12, and 18 years before the M6.5 shock (inset of Fig. 5)
reveals that the inverse power law is an approximation
for data at distances around Δr ≥ 40 km (Δrc = 40 km),with density levels above the background implied fromdistant seismicity (gray shading). An inverse power law
with η = 1.35 was clearly observed for the 3-year period,
although the exponent η showed a gradual increase overthe past 18 years. Data of Δrc for the Kumamoto case areincluded in Fig. 4. Although this figure still has a slightly
large scatter, it is reasonable that a positive correlation
between Δrc and m is comparable with the la − m relation. The similar dependences on m indicate that the sizeof the spatial organization of seismicity is governed by
the size of asperity (see Table 1 for a summary).

Discussion and conclusions
Although long-term probabilistic seismic hazard assessment can be made in Japan (e.g., National Seismic Hazard Maps for Japan) and other seismically active regions,
it is generally accepted that immediate local precursory

Nanjo et al. Earth, Planets and Space (2016) 68:187

phenomena are not seen ubiquitously. In order to test
the hypothesis of a local precursor for a fault system,
the Parkfield Earthquake Predictability Experiment was
initiated in 1985. The expected Parkfield earthquake
occurred beneath the heavily instrumented region on
September 28, 2004. No local precursory changes wereobserved (Lindh 2005). In the absence of local precursory signals, the next question is whether broader
anomalies develop and in particular whether there isanomalous seismicity.
Before attempting to identify precursory phenomena,a fundamental question is whether there are different
forecast methods applicable to seismicity in a region ofinterest, in particular in Japan. This is the question that
we addressed in this paper. We decided to use the b-value
method, the RTM-algorithm, seismicity density analysis,
and the Z-value method and applied them to seismicity
before the Kumamoto earthquakes. Our study is the firstreport on this theme after completing the modernization of seismic monitoring in Japan. Before this, an earlyattempt was made by Enescu and Ito (2001), who used
the b-value method, the Z-value method, and fractal
dimensions to discover premonitory quiescence followedby activation of seismicity about 2 years before the occurrence of the 1995 M7.2 Kobe earthquake. These authors
conducted their study before the completion of the modernization of seismic monitoring in Japan and used the
local network data of Kyoto University, besides the JMA
catalog data, to ensure a relatively low and stable magnitude of completeness during the whole analyzed time
period.
The combination of multiple methods may enhance
the reliability of the revealed anomalies when comparedto the reliability of a single method. We performed somecomprehensive comparisons of the results from eachmethod to support this claim. The properties of precursory episodes are summarized in Table 1. Estimates ofthe duration and spatial extent contain uncertainties,
which depend on the approaches taken in the analysis.
The durations were almost identical in both the RTMalgorithm and the Z-value method. In these methods,
the resulting maximum expression of the anomalies did
not coincide with the epicenters of the Kumamoto mainshock, nor of the Off Satsuma Peninsula mainshock.
Furthermore, both methods did not identify the same
locations for the anomalous maximum. Wyss et al. (2004)
indicated a reason for the difference between methods in
pinpointing the location of an anomaly. The weighting ofresults based on the size of the earthquakes, which is only
done in the RTM method but not in the Z-value method,results in some differences in the estimated significancein most samples. Therefore, the maximum expressionof the anomaly was not observed at exactly the same

Page 8 of 10

locations. The reason for the maximum expression of
anomaly not pointing to the epicenters may be due to thefact that the onset of both earthquakes may have been
at a segment along the tectonic line from the Okinawa
Trough (Ishibashi 2016). The large area of the anomaly
may reflect the nature of the process leading up to thephenomenon in the regional tectonic condition. Nonetheless, the two methods based on different assumptions,different algorithms, and different definitions of anomaly,
arrived at very similar results. This strongly suggests thatthe observed anomalies are real, and can be determinedwith considerable reliability.
In contrast to these two methods, the b-value method
and seismicity density analysis seem to have narrower
spatial extents (Table 1). The close match found in spatial extent between the different approaches implies thatasperities (highly stressed patches) may act as an indicatorof the preparation process to an impending earthquake. A
duration of 3 years or longer from the seismicity density
analysis (Table 1) is longer than that revealed by both the
RTM-algorithm and the Z-value method. The precursory
duration from application of the b-value method to analyzed data for the past 16 years is unknown. The timing of
the earthquake remains uncertain from low b-values andspatial organization prior to the Kumamoto earthquakes.
This is consistent with another observation of active faults
in which the heterogeneous pattern of b-values at Parkfield was, to a high degree, stationary for the past 35 years
and the 2004 M6 earthquake eventually occurred at azone of low b-values (Schorlemmer and Wiemer 2005). A
decrease in b tracking stress buildup, as expected from alaboratory experiment (Scholz 1968), may not be observable for the decade-scale observation of active faults inJapan and California, because it is too short to observe
such a decrease in b. Instead, the stationary nature of
b-values, as observed in Fig. 1, is reasonable. The mechanism of stress buildup within the fault zones is uncertain,
but one hypothesis is that a steady slip of the deeper continuation of faults that does not produce earthquakes, but
still involves motions across the fault, forces the upper
crust around the faults to deform and thus concentratestress. This is still difficult to measure directly.
Overall, our findings indicate that the methods weadopted do not allow the Kumamoto earthquakes to bepredicted exactly. The spatial extent of precursory patterns that were detected differed from one method toanother and ranged from local scales (equivalent to anasperity dimension) to regional scales (typically, 100–
200 km). The Kumamoto earthquakes were preceded byperiods of pronounced anomalies, which lasted yearly
scales (1.5 years) or longer (>3 years). Given the widelydifferent scales of anomalies in time and space, a combination of multiple methods was able to detect different

Nanjo et al. Earth, Planets and Space (2016) 68:187

signals prior to the Kumamoto earthquakes. This strongly
suggests a great potential to reduce the hazard at possible future sites of earthquakes relative to long-term seismic hazard assessment. We also found that the seismic
quiescence in a regional-scale area, detected by using
the RTM-algorithm and the Z-value method, was a common precursor to the 2016 Kumamoto earthquakes and2015 Off Satsuma Peninsula earthquake that occurred
5 months before the 2016 Kumamoto earthquakes. The
precursory durations of 1 year to the 2015 Off Satsuma
Peninsula earthquake were almost identical in both
approaches. The result allows us to interpret both eventsas the onset that occurred at a section along the tectonic
line from the Okinawa Trough through the Beppu–
Shimabara graben (Ishibashi 2016).

Additional file
Additional file 1. This file includes Appendixes 1–2, Figures A1–A12,
Table A1, and additional references.

Authors’ contributions
All the authors performed numerical simulations, analyzed data, and prepared
the figures, in particular b-value method (KZN, JI), RTM-algorithm (TN), Z-value
method (KZN), seismicity density analysis (YO, ST, HN, TO, RT, MK), and application of declustering programs (KZN, NF, TN). KZN and TN helped to draft themanuscript. KZN wrote the final manuscript. All authors read and approved
the final manuscript.
Author details
1
Global Center for Asian and Regional Research, University of Shizuoka, 3‑6‑1,
Takajo, Aoi, Shizuoka 420‑0839, Japan. 2 International Digital Earth Applied
Science Research Center, Chubu Institute for Advanced Studies, Chubu
University, 1200 Matsumoto‑cho, Kasugai, Aichi 487‑8501, Japan. 3 Department of Physics, Tokyo Gakugei University, 4‑1‑1, Nukuikita, Koganei, Tokyo184‑8501, Japan. 4 Institute of Oceanic Research and Development, Tokai
University, 3‑20‑1, Orido, Shimizu, Shizuoka 424‑0902, Japan.
Acknowledgements
The authors thank the editor (M. Hashimoto) and two anonymous referees
for their useful comments. This study was partly supported by the Ministry of
Education, Culture, Sports, Science and Technology (MEXT) of Japan, under
its Earthquake and Volcano Hazards Observation and Research Program (KZN,
JI, YO, MK, TN) and a Grant-in-Aid for Scientific Research (C), No. 26350483,
2014–2017 (MK), and the Collaboration Research Program of IDEAS, Chubu
University IDEAS201614 (JI).
Competing interests
The authors declare that they have no competing interests.
Received: 1 August 2016 Accepted: 29 October 2016

References
Enescu B, Ito K (2001) Some premonitory phenomena of the 1995 Hyogo-Ken
Nanbu (Kobe) earthquake: seismicity, b-value and fractal dimension.
Tectonophysics 338(3–4):297–314. doi:10.1016/S0040-1951(01)00085-3
Felzer KR, Brodsky EE (2006) Decay of aftershock density with distanceindicates triggering by dynamic stress. Nature 441:735–738. doi:10.1038/
nature04799

Page 9 of 10

Field EH (2007) Overview of the working group for the development ofregional earthquake likelihood models (RELM). Seismol Res Lett 78:7–16.
doi:10.1785/gssrl.78.1.7
Frohlich C, Davis SD (1990) Single-link cluster analysis as a method to evaluate
spatial and temporal properties of earthquake catalogues. Geophys J Int
100(1):19–32. doi:10.1111/j.1365-246X.1990.tb04564.x
Gutenberg B, Richter CF (1944) Frequency of earthquakes in California. Bull
Seismol Soc Am 34(4):185–188
Huang Q (2004) Seismicity pattern changes prior to large earthquakes: anapproach of the RTL algorithm. Terr Atmos Ocean Sci 15(3):469–491
Huang Q (2006) Search for reliable precursors: a case study of the seismic
quiescence of the 2000 western Tottori prefecture earthquake. J Geophys
Res 111:B04301. doi:10.1029/2005JB003982
Huang Q, Nagao T (2002) Seismic quiescence before the 2000 M = 7.3 Tottori earthquake. Geophys. Res. Lett. 29(12):19-1–19-4. doi:10.1029/200
1GL013835
International Commission on Earthquake Forecasting for Civil Protection(2011) Operational earthquake forecasting. State of knowledge andguidelines for utilization. Ann Geophys. doi:10.4401/ag-5350
Ishibashi K (2016) The 2016 Kumamoto earthquakes are not exceptional.
Kagau (Science) 86(6):532–540 (In Japanese)
Jordan T (2006) Earthquake predictability, brick by brick. Seismol Res Lett
77:3–6
Kanamori H (2003) Earthquake prediction: An overview. In: Lee WHK, Kanamori
H, Jennings PC, Kisslinger C (eds) International handbook of earthquake &engineering seismology. Academic Press, Amsterdam, pp 1205–1216
Katsumata K (2015) A long-term seismic quiescence before the 2004
Sumatra (Mw 9.1) earthquake. Bull Seismol Soc Am 105(1):167–176.
doi:10.1785/0120140116
Keilis-Borok V (2002) Earthquake predictions: state-of-the-art and emerging
possibilities. Ann Rev Earth Planet Sci 30:1–33. doi:10.1146/annurev.
earth.30.100301.083856
Lindh AG (2005) Success and failure at Parkfield. Bull Seismol Soc Am 76(1):3–6.
doi:10.1785/gssrl.76.1.3
Lippiello E, Marzocchi W, de Arcangelis L, Godano C (2012) Spatial organization of foreshocks as a tool to forecast large earthquakes. Sci Rep 2:846.
doi:10.1038/srep00846
Mogi K (1962) Magnitude-frequency relation for elastic shocks accompanying
fractures of various materials and some related problems in earthquakes.
Bull Earthq Res Inst 40:831–853
Nagao T, Takeuchi A, Nakamura K (2011) A new algorithm for the detection ofseismic quiescence: introduction of the RTM algorithm, a modified RTL
algorithm. Earth Planets Space 63:315–324. doi:10.5047/eps.2010.12.007
Nanjo KZ, Ishibe T, Tsuruoka H, Schorlemmer D, Ishigaki Y, Hirata N (2010) Analysis of the completeness magnitude and seismic network coverage ofJapan. Bull Seismol Soc Am 100(6):3261–3268. doi:10.1785/0120100077
Nanjo KZ, Hirata N, Obara K, Kasahara K (2012) Decade-scale decrease inb value prior to the M9-class 2011 Tohoku and 2004 Sumatra quakes.
Geophys Res Lett 39:L20304. doi:10.1029/2012GL052997
National Research Institute for Earth Science and Disaster Resilience (2016)
http://www.kyoshin.bosai.go.jp/kyoshin/topics/Kumamoto_20160416/
inversion/. Access 20 Sept 2016
Obara K, Kasahara K, Hori S, Okada Y (2005) A densely distributed high-sensitivity seismograph network in Japan: Hi-net by National Research Institutefor Earth Science and Disaster Prevention. Rev Sci Instrum 76:021301.
doi:10.1063/1.1854197
Richards-Dinger K, Stein RS, Toda S (2010) Decay of aftershock densitywith distance does not indicate triggering by dynamic stress. Nature
467:583–586. doi:10.1038/nature09402
Roberts NS, Bell AF, Main IG (2015) Are volcanic seismic b-values high, andif so when? J Volcanol Geotherm Res 308:127–141. doi:10.1016/j.
jvolgeores.2015.10.021
Scholz CH (1968) The frequency-magnitude relation of microfracturing in rockand its relation to earthquakes. Bull Seismol Soc Am 58(1):399–415
Scholz CH (2002) The mechanics of earthquakes & faulting, 2nd edn. Cambridge University Press, Cambridge
Scholz CH (2015) On the stress dependence of the earthquake b value. Geophys Res Lett 42:1399–1402. doi:10.1002/2014GL062863
Schorlemmer D, Wiemer S (2005) Microseismicity data forecast rupture area.
Nature 434:1086. doi:10.1038/4341086a

Nanjo et al. Earth, Planets and Space (2016) 68:187

Shi Y, Bolt BA (1982) The standard error of the magnitude-frequency b-value.
Bull Seismol Soc Am 72:1677–1687
Sobolev GA (2011) Seismicity dynamics and earthquake predictability. Nat
Hazards Earth Syst Sci 11:445–458. doi:10.5194/nhess-11-445-2011
Sobolev GA, Tyupkin YS (1997) Low-seismicity precursors of large earthquakes
in Kamchatka. Volcanol Seismol 18:433–446
Sobolev GA, Tyupkin YS (1999) Precursory phases, seismicity precursors, andearthquake prediction in Kamchatka. Volcanol Seismol 20:615–627
Sommerville P, Skarlatoudis A, Thio HK (2015) Source scaling relations of subduction earthquakes for strong ground motion and tsunami prediction.
http://earthquake.usgs.gov/research/external/reports/G13AP00028.pdf.
Accessed 20 Sept 2016
Tada T (1985) Spreading of the Okinawa Trough and its relation to the crustal
deformation in the Kyushu (2). Zisin 38(1):1–12 (In Japanese withEnglish abstract)

Page 10 of 10

Tormann T, Enescu B, Woessner J, Wiemer S (2015) Randomness of megathrustearthquakes implied by rapid stress recovery after the Japan earthquake.
Nat Geosci 8:152–158. doi:10.1038/ngeo2343
Wiemer S, Wyss M (1994) Seismic quiescence before the Landers (M = 7.5) and
Big Bear (M = 6.5), earthquakes. Bull Seismol Soc Am 84(3):900–916
Woessner J, Wiemer S (2005) Assessing the quality of earthquake catalogues:
Estimating the magnitude of completeness and its uncertainty. Bull
Seismol Soc Am 95:684–698. doi:10.1785/0120040007
Wyss M, Sobolev G, Clippard JD (2004) Seismic quiescence precursors to two
M7 earthquakes on Sakhalin Island, measured by two methods. Earth
Planets Space 56(8):725–740. doi:10.1186/BF03353081
Yagi Y, Okuwaki R, Enescu B, Kasahara A, Miyakawa A, Otsubo M (2016) Rupture
process of the 2016 Kumamoto earthquake in relation to the thermal
structure around Aso volcano. Earth Planets Space 68:118. doi:10.1186/
s40623-016-0492-3