An application of earthquake prediction algorithm M8 in eastern Anatolia at the approach of the 2011

Van earthquake
Masoud Mojarab1, Vladimir Kossobokov2,3, Hossein Memarian1,∗ and Mehdi Zare4
1School of Mining Engineering, University of Tehran, Tehran, Iran.
2Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences,
Moscow, Russian Federation.
3Institut de Physique du Globe de Paris, Paris, France.
4Seismology Research Center, International Institute of Earthquake Engineering and Seismology (IIEES),
Tehran, Iran.
∗Corresponding author. e-mail: memarian@ut.ac.ir
On 23rd October 2011, an M7.3 earthquake near the Turkish city of Van, killed more than 600 people,
injured over 4000, and left about 60,000 homeless. It demolished hundreds of buildings and caused
great damages to thousand others in Van, Ercis, Muradiye, and C¸ aldıran. The earthquake’s epicenter is
located about 70 km from a preceding M7.3 earthquake that occurred in November 1976 and destroyed
several villages near the Turkey–Iran border and killed thousands of people. This study, by means of
retrospective application of the M8 algorithm, checks to see if the 2011 Van earthquake could have
been predicted. The algorithm is based on pattern recognition of Times of Increased Probability (TIP)
of a target earthquake from the transient seismic sequence at lower magnitude ranges in a Circle of
Investigation (CI). Specifically, we applied a modified M8 algorithm adjusted to a rather low level of
earthquake detection in the region following three different approaches to determine seismic transients.
In the first approach, CI centers are distributed on intersections of morphostructural lineaments recognized
as prone to magnitude 7+ earthquakes. In the second approach, centers of CIs are distributed on
local extremes of the seismic density distribution, and in the third approach, CI centers were distributed
uniformly on the nodes of a 1◦ × 1◦ grid. According to the results of the M8 algorithm application, the
2011 Van earthquake could have been predicted in any of the three approaches. We noted that it is possible
to consider the intersection of TIPs instead of their union to improve the certainty of the prediction
results. Our study confirms the applicability of a modified version of the M8 algorithm for predicting
earthquakes at the Iranian–Turkish plateau, as well as for mitigation of damages in seismic events in
which pattern recognition algorithms may play an important role.
1. Introduction
In the area under study (figure 1), the 23 October
2011, M7.3 Van earthquake had more than
600 victims and caused massive damages in several
cities, including Van, Ercis, Muradiye and C¸ aldıran
(Zare et al. 2011). The hypocenter was located near
the Tabanli village at a depth of 16 km (USGS
2011) (figure 2). The earthquake resulted from the
movement along the 28 km segment of the Van fault
Keywords. Earthquake prediction; pattern recognition; M8 algorithm; time of increased probability (TIP); circle of
investigation (CI); Van earthquake.
J. Earth Syst. Sci.
c
Indian Academy of Sciences
Masoud Mojarab et al.
Figure 1. The area under study. The study area is located on Google map with a black border.
Figure 2. Seismic settings in the study area. (1) Epicenter of the 23 November 2011 Van earthquake; (2) area under study;
(3) historical earthquakes which are distributed in the area under study; (4) epicenter of earthquakes reported by the
USGS/NEIC Global Hypocenters’ Data Base (GHDB) system in the period of 1900–2010; (5) Van fault with 28 km length.
with the north aspect slope and was felt not only
in eastern Turkey, but also as far as northwestern
Iran (Emre et al. 2011).
The seismic history of the eastern Anatolian
region before the 2011 Van earthquake is
summed up in table 1 (Ambraseys 1988, 2009).
The 24 November 1976, M7.3 C¸ aldıran-Muradiye
earthquake occurred 20 km northeast of Muradiye,
in the Van province of eastern Turkey, just about
70 km far from the epicenter of the 2011 Van earthquake.
The earthquake also had a maximum intensity
X on the Modified Mercalli Intensity (MMI)
scale. The area of severe damage, where over 80%
of the buildings were destroyed, covered 2000 km2.
Earthquake prediction algorithm M8 in eastern Anatolia
There were about 5000 victims (Gulkan et al. 1978).
The 6 May 1930, M7.2 Salmas earthquake occurred
further to the east from the 2011 epicenter. The
earthquake measured 7.2 on the Richter scale and
7.4 surface wave magnitude and resulted in 2500
direct fatalities. One foreshock occurred prior to
the rupture, and multiple aftershocks followed the
main event (Tchalenko and Berberian 1974).
Thus, it is evident that eastern Anatolian zone
shows potential of an earthquake with magnitude 7
and more. It raises a question whether the time of
such earthquakes in the region are predictable
by application of reliable algorithms? United
States National Research Council, Panel on Earthquake
Prediction of the Committee on Seismology
recommended the following definition of earthquake
prediction (Allen et al. 1976):
“An earthquake prediction must specify the
expected magnitude range, the geographical
area within which it will occur, and the time
interval within which it will happen with sufficient
precision so that the ultimate success
or failure of the prediction can readily be
judged. Only by careful recording and analysis
of failures as well as successes can the eventual
success of the total effort be evaluated
and future directions charted. Moreover, scientists
should also assign a confidence level to
each prediction.”
Table 1. Historical earthquakes in eastern Anatolian region (Ambraseys 1988, 2009).
No. Year Month Day Lat. Long. Ms Io Earthquake name
1 1983 Oct 30 40.35 42.18 6.7 VIII Narman
2 1976 Nov 24 39.1 44 7.1 IX Caldiran
3 1975 Sep 6 38.51 40.77 6.6 IX Lice
4 1972 Jun 16 38.25 43.35 4.9 VIII Van
5 1971 May 22 38.92 40.64 6.8 IX Bing˝ol
6 1966 Aug 19 39.17 41.56 6.8 IX Varto
7 1959 Oct 25 39.14 41.6 4.9 VIII Varto
8 1949 Aug 17 39.4 40.65 6.9 IX Elmalidere
9 1945 Nov 20 38.44 43.39 4.8 VII Van
10 1941 Sep 10 39.13 43.12 5.9 VIII Ercis
11 1941 Nov 12 39.85 39.35 5.9 VIII Erzincan
12 1939 Des 26 39.8 39.38 7.8 XI Erzincan
13 1930 May 6 38.25 44.6 7.2 X Salmas
14 1924 Sep 13 40.05 42.3 6.8 IX Horasan
15 1903 Apr 28 39.14 42.65 7 X Patnos
Table 2. Classification of earthquake prediction accuracy. Approximate temporal time
of earthquake occurrence being related to the rupture size L of the incipient earthquake.
Temporal, in years Spatial, in source zone size (L)
Long-term 10 Long range Up to 100
Intermediate-term 1 Middle range 5–10
Short-term 0.01–0.1 Narrow range 2–3
Immediate 0.001 Exact 1
Table 3. Effectiveness of the algorithm M8 (after Kossobokov 2013).
Target earthquakes Alarm volume p Confidence
Test period Predicted Total (%) level (%)
Prediction of earthquakes with magnitude M8.0+
1985–2012 16 21 32.84 99.99
1992–2012 14 19 29.80 99.99
Prediction of earthquakes with magnitude M7.5+
1985–2012 40 68 28.73 99.99
1992–2012 30 56 23.14 99.99
Masoud Mojarab et al.
The earthquake prediction algorithms, based on
pattern recognition methods, offer an answer to the
question stated above, by indicating the times of
increased probability, TIPs, in a specified range of
time, location, and magnitude. One of the important
time and location classifications is illustrated
in table 2. It should be noted that not only the
‘short-term exact’ class, but a wide variety of other
possible combinations might be useful for mitigation
of seismic risks. Moreover, having in mind the
complexities of the Earth lithosphere, in particular,
its blocks-and-faults structure and their evidently
non-linear dynamics, the ‘short-term exact’ prediction
accuracy might be unreachable, in principle
(Hough 2009).
According to table 2, the times of increased
probability (TIPs) diagnosed by the M8 algorithm
(briefly described in appendix A), whose original
version was published in 1987 (Keilis-Borok and
Kossobokov 1987) and is presently available in
the IASPEI software library (Kossobokov 1997),
belong to the class of intermediate-term middlerange
predictions. The algorithm has passed a
rigid testing in the ongoing global experiment
started in 1992 (Healy et al. 1992). The confidence
level achieved in more than 20 years of the
ongoing global testing (table 3), along with the
results of retrospective studies targeting magnitude
ranges other than M8.0+ and M7.5+ (Keilis-
Borok and Kossobokov 1990; Kossobokov 2013),
encourages the M8 algorithm real-time applications
on a regional scale. In this study, we attempt
to justify such an application in the region of the
Turkish–Iranian plateau.
2. Anatolian and Persia–Tibet plateau
Recognition of discrete plate tectonics between
Africa and Eurasia was reported by McKenzie
(1972). He combined seismicity, geology and topographic
data and proposed that Anatolian plate is
escaping towards west due to Arabia–Eurasia collision.
Recent GPS investigations proved McKenzie’s
model (Bird 2003). More recently, a global study
tried to identify tectonic boundaries for the entire
world; for which most of information was collected
from literatures (Bird 2003) (figure 3). A global
study shows that continental lithosphere in Middle
East is divided into Anatolian, Arabian, Eurasia,
and Persia–Tibet discrete plates, while many others
emphasized independency of Anatolian plate
(McKenzie 1972; Bird 2003; McClusky et al. 2000;
Reilinger et al. 2006). In the present study, the area
of investigation of seismic activity in advance the
2011 Van earthquake has been selected between
Anatolian, Arabia, and Tibet plates, which mostly
cover the Persia plate. Persia and Tibet are
the two main plateaus of the Alpine-Himalayan
collision system (Seng¨or and Kidd 1979; Dewey
Figure 3. Tectonic settings in the study area. (1) Epicenter of the 23 November 2011 Van earthquake; (2) area under
study; (3) I-rank lineaments (after Gelfand et al. 1974a); (4) active faults. NAFZ: North Anatolian fault zone. EAFZ: East
Anatolian fault zone.
Earthquake prediction algorithm M8 in eastern Anatolia
et al. 1986). The Persia plateau extends from
eastern Anatolia to eastern Iran, and typically has
elevations of 1.5–2 km, decreasing to about 500 m
in eastern Iran. The basement of Persia plateau
consists of micro-plates interspersed with zones of
ophiolites and melanges accreted to each other,
formed on this part of Eurasia by the late Cretaceous
or early Tertiary (Seng¨or 1990; Berberian
et al. 1982). Westward transport of the Anatolian
plate occurs by slip on the East and North Anatolian
faults (McKenzie 1972).
The left-lateral East Anatolian fault accommodates
motion between the Arabian and Anatolian
plates. Historical slip rates, based on seismic
moments of earthquakes along the East Anatolian
fault, are 6–10 mm yr−1 (Taymaz et al. 1991).
The GPS-derived slip rate is 9 ± 1 mm yr−1 and
slip rate for the North Anatolian fault is 24 ± 1
mm yr−1 (McClusky et al. 2000) (figure 3).
The focal mechanism of earthquakes in eastern
Anatolian plateau has a dominating northeast–
southwest trend, with a strong strike slip component.
The compressive component increases by approaching
to Zagros suture zone. The recent 23 October
2011 earthquake in Van was not an exception from
a series of compressive events assigned to the Persia
plateau. Based on field investigation in Van region,
this earthquake belongs to an exposed active reverse
fault with 28 km length (figure 2, zoomed panel),
known as the Van fault (Emre et al. 2011).
3. Applying M8 algorithm
The main objective of this paper is to demonstrate
predictability of major earthquakes on the
territory of the Iranian–Turkish plateau, exemplified
by a case-study of the 23 November 2011 Van
earthquake and the M8 algorithm. It should be
noted that it is of crucial importance to determine
the application areas before running any prediction
algorithm in such a way that the initial
conditions are satisfied. Specifically, the seismicity
rate in a Circle of Investigation (CI), and distribution
of the CI centers in the area under study
has a considerable impact on the results of the M8
algorithm application. If we determine CI centers
uniformly, without any guidance from tectonics,
seismic activity, pattern recognition, etc., our calculations
would be exhaustive, though not much
informative in many cases and may be even misleading
in some cases, e.g., where the actual size of
earthquake-prone locus inside a CI is much smaller
than its diameter. Hence, different approaches were
applied to design an appropriate distribution of CI
centers.
Figure 4. The D-intersections of morphostructural lineaments which have been selected as CIs centers, prone to earthquake
of magnitude 6.5+ and 7.0+. (1) Epicenter of the 23 November 2011 Van earthquake; (2) area under study; (3) I-rank
lineaments (after Gelfand et al. 1974a); (4) active faults; (5) D-intersections of the lineaments (CIs centers) recognized as
prone to earthquake of magnitude 6.5+ and 7+ (Gelfand et al. 1974b); (6–7) the 40 and 62.5 km outline of the D-intersections
prone to earthquake of magnitude 6.5+ and 7.0+ (Kossobokov and Rotwain 1977), respectively.
Masoud Mojarab et al.
The first approach makes use of pattern recognition
of intersections of morphostructural lineaments
with the potential of generating earthquakes of
magnitude 6.5+ and 7+ (figure 4). These intersections
were introduced by Gelfand et al. (1974a, b)
as dangerous, D-intersections, prone to magnitude
Figure 5. The empirical density of earthquake epicenter distribution. Earthquake density is the number of earthquakes
in total duration of a study area catalog. (1) Epicenter of the 23 November 2011 Van earthquake; (2) area under study;
(3) I-rank lineaments (after Gelfand et al. 1974a); (4) active faults; (5) centers of the circles of investigation associated with
the seismic density distribution patterns.
Figure 6. Regular grid of the centers of the CIs for the study area. (1) Epicenter of the 23 November 2011 Van earthquake;
(2) area under study; (3) I-rank lineaments (after Gelfand et al. 1974a); (4) active faults; (5) centers of the circles of
investigation.
Earthquake prediction algorithm M8 in eastern Anatolia
6.5+ earthquakes. Further investigations have shown
that some of these intersections have the potential
to produce larger earthquakes with magnitudes 7+
(Kossobokov and Rotwain 1977). In the mentioned
studies, the pattern recognition determination of
D(6.5+)- and D(7.0+)-intersections is based on
their geomorphologic description and geographic
parameters measured in circular areas of the radii
comparable to linear dimensions of M6.5+ and
M7.0+ earthquakes (i.e., 25–40 and 40–62.5 km,
respectively).
In the second approach, centers of CIs are distributed
on local extremes of the seismic density
distribution. The map of the associated frequency
of earthquakes (figure 5) has been prepared by
means of the Kernel method (Silverman 1986).
In our case, the Kernel method counts the number
of earthquakes in a circular area with a prefixed
search radius centered on a pixel and assigns
the obtained number to that pixel. In accordance
with pattern recognition results (Kossobokov and
Rotwain 1977), the search radius in the present
study is selected to be 62.5 km.
In the third approach, a regular grid of 1◦ × 1◦
was selected as centers of CIs, regardless of the seismicity
distribution and morphostructures (figure 6).
4. Data
Catalogs of earthquakes remain the most objective
record of seismic activity; that is why most
of the studies on seismic precursory phenomena
and earthquake prediction, are based on the analysis
of earthquake catalogs (Peresan et al. 2005).
The earthquake catalog of eastern Anatolian and
western Persia plateau regions at latitudes of 33◦–
44◦N and longitudes of 36◦–49.5◦E, was selected
from the USGS/NEIC Global Hypocenters’ Data
Base System, GHDB. Following numerical recipes
by Shebalin (1992), special searches for duplicates
and possible errors were performed. The analysis
of the frequency–magnitude graph for considering
the consistency of the catalog was accomplished
by inspection of the cumulative distribution of
earthquake frequency vs. years (figure 7) and the
Gutenberg–Richter graph (figure 8). The number
of earthquakes vs. years helps recognizing the stability
of earthquake occurrence in the area under
study. In the 1960’s, the quality and completeness
of the global catalogs were improved with an installation
of the Worldwide Seismic Network. Since
1963, all the earthquakes of magnitude 5 or greater
are located and in many areas catalogs are complete
to magnitude 4. As evident from figure 7, a
step in the number of earthquakes appear around
1965, confirming that the GHDB catalog in the
study territory is grossly incomplete at magnitude
4.5 before 1965 and reaches reasonable completeness
at this level in about one decade. Following
this observation, we decided to use the data
from 1975 for an application of the modified version
of M8 algorithm, although according to the
Gutenberg–Richter plots (figure 8), the level of
completeness is about magnitude 4.6 in this period
of time. (Here we do not investigate the nature of a
Figure 7. The annual number of earthquakes in the GHDB system in the area of investigation, 1900–2012. Each bound
corresponds to a half-magnitude range. These bounds are stacked from higher ranges. Vertical red lines indicate the interval
used for the Gutenberg–Richter plot.
Masoud Mojarab et al.
Figure 8. The Gutenberg–Richter plots for the area under study, 1975–June 2011. (1) Annual number of earthquake in
a certain magnitude bin; (2) cumulative number of earthquakes summed from higher magnitude ranges; (3) maximum
likelihood solution of Log(N) = a − b M (b = 1.11, a = 6.36). Possible choices of the magnitude of completeness are
indicated with a vertical line and an arrow.
remarkable bulge of the plots at above magnitude
6.0 which may signify the overcritical state of
earthquake activity on the territory considered in
the period associated with seismic extremes.)
5. Results
The modified version of M8 algorithm targeting
earthquakes from M7.0+ magnitude range was
applied to the entire study area following the three
approaches described above. In each case, the 23
October 2011, M7.3 earthquake was used to control,
retrospectively, the efficiency of the modified
version of M8 algorithm. To adjust an application
of M8 algorithm to the local level of the catalog
completeness, we have modified the two constants
that determine lower magnitude cut-offs (see
Appendix A) from ˜N = 20 and ˜N= 10 seismic
events per year to ˜N = 10 and ˜N = 5, respectively.
Note that the order of the M8 algorithm parameters
does not change the results of TIP diagnosis,
therefore, the modification differs by a single one
parameter of the original algorithm, i.e., ˜N = 20
changed for ˜N = 5.
5.1 First approach
Application of the modified M8 algorithm in CIs
centered at D-intersections with potential magnitude
7+ shows that a TIP has been diagnosed from
01/07/2008 to 02/07/2013 in a single one CI shown
in figure 9. The 2011 Van earthquake occurred in
this area of alarm and also in the 62.5 km outline
of the D-intersections prone to earthquake of magnitude
7.0+ (Kossobokov and Rotwain 1977). The
graphs of the functions of the M8 algorithm in the
period from 1985 to 2011 are shown in figure 10.
All the seven functions involved in the diagnosis
of TIP did rise to their anomalous values by mid-
2008. Accordingly, the 2011 Van earthquake has
been predictable using the first approach with an
accuracy of 3–4 years and about 200 km.
As illustrated in figure 11, the results of the
M8 algorithm application targeting M6.5+ earthquake
are formally confirmed by the 2011 Van
earthquake, although its epicenter is located away
from the triple intersection of CIs in TIPs. As evident
from figure 11, the CIs centered at D(6.5+)-
intersections nearest to the M7.3 epicenter did not
show a TIP. Both remarks may reflect the accuracy
of the M8 predictions and selectivity of this
approach in relation to magnitude range: the actual
magnitude of M7.3 fits better the range of M7.0+
than that of M6.5+.
5.2 Second approach
The same modified version of the M8 algorithm was
applied to seismic activity in CIs centered on top
of the peaks of earthquake occurrence density. In
application targeting M7.0+ earthquakes 8 out of
15 CIs were in state of TIP (figure 12). It is notable
that the 2011 Van earthquake epicenter falls into
the 8-fold intersection of all the alerted CIs. Thus,
following the second approach, the Van earthquake
Earthquake prediction algorithm M8 in eastern Anatolia
Figure 9. The results of the M8 algorithm application targeting M7.0+ earthquake with 281 km radius of investigation, as
of July 2011: morphostructural intersections. (1) Epicenter of the 23 November 2011 Van earthquake; (2) area under study;
(3) I-rank lineaments (after Gelfand et al. 1974a); (4) active faults; (5) D-intersections of the lineaments which are used
as CIs centers (after Gelfand et al. 1974b); (6–7) the 40 and 62.5 km outline of the D-intersections recognized as prone to
earthquake of magnitude 6.5+ and 7+, respectively (Kossobokov and Rotwain 1977). Area of TIP diagnosed for the second
half-year of 2011 is shaded gray. The center of CI with a TIP marked green.
Figure 10. Measures of the M8 algorithm targeting M7.0+ in a CI with a TIP. Functions are normalized to an arbitrary
unit ranging from minimal value of 0 to maximal 1. The anomalous values are marked with red dots.
Masoud Mojarab et al.
Figure 11. The results of the M8 algorithm application targeting M6.5+ earthquake with 192 km radius of investigation, as
of July 2011: morphostructural intersections. (1) Epicenter of the 23 November 2011 Van earthquake; (2) area under study;
(3) I-rank lineaments (after Gelfand et al. 1974a); (4) active faults; (5) D-intersections of the lineaments (after Gelfand et al.
1974b); (6–7) the 25 and 40 km outlines of the D-intersections recognized as prone to earthquake of magnitude 6+ and
6.5+, respectively. The centers of CIs with TIPs diagnosed for the second half-year of 2011 marked green; the area of TIPs
shaded gray.
Figure 12. The results of the M8 algorithm application targeting M7+ earthquake with 281 km radius of investigation, as
of July 2011: seismic density traces. (1) Epicenter of the 23 November 2011 Van earthquake; (2) area under study; (3) Irank
lineaments (after Gelfand et al. 1974a); (4) active faults; (5) centers of the circles of investigation associated with the
seismic density distribution patterns. The centers of CIs with TIPs diagnosed for the second half-year of 2011 are marked
green; the area of TIPs shaded gray. The darkest gray area is the intersection of CIs in TIPs.
Earthquake prediction algorithm M8 in eastern Anatolia
Figure 13. The results of the M8 algorithm application targeting M6.5+ earthquake with 192 km radius of investigation,
as of July 2011: seismic density traces. (1) Epicenter of the 23 November 2011 Van earthquake; (2) area under study;
(3) I-rank lineaments (after Gelfand et al. 1974a); (4) active faults; (5) centers of the circles of investigation associated
with the seismic density distribution patterns. The center of CIs with TIPs diagnosed for the second half-year of 2011 are
marked green; the area of TIPs shaded gray.
could have been predicted using the modified M8
algorithm with the intermediate-term accuracy in
time and about 300 km of this intersection.
A TIP targeting magnitude 6.5+ is diagnosed in
the CI centered at the peak of the seismic density
map nearest to and almost at the 2011 Van earthquake
epicenter. However, similar to the results of
the first approach, the epicenter falls off the 7-fold
intersection of the 7 CIs in state of magnitude 6.5+
alarm (figure 13).
5.3 Third approach
The modified M8 algorithm was also applied to a
set of CIs centered at the grid points of a regular
mesh that covers the territory under study. For 8
out of 42 CIs, the number of earthquakes reported
in GHDB is insufficient to run even the modified
version of the M8 algorithm (figure 14). The diagnosis
of TIPs targeting M7.0+ in the remaining 34
CIs is positive for the 11 CIs, which cover almost
the entire territory under study. By comparing
the results obtained in the first two approaches,
we have to conclude that running the M8 algorithm
on a regular grid that does not account for
tectonics either in terms of specific earthquakeprone
morphostructures or empirical earthquake
locations could lead to erroneous judgment
about the dynamic changes of seismic hazard and
associated risks. Same as in the first two approaches,
the results of the modified M8 algorithm
application targeting magnitude 6.5+ earthquakes
(not shown here) favour a trivial conclusion that
predicting a magnitude M7.3 event is better when
targeting the range of M7.0+ than that of M6.5+.
5.4 Stability of the results
The results of the modified M8 algorithm application
as of the 1st of July 2011, described above are
obtained with the prefixed starting point of analysis
To = 1975, the beginning Tb = 1987, and
the ending Te = 2011 of data used for diagnosis
of TIP. In this way the periods of 12 and 24 years
between To and Tb and between Tb and Te are
enough to obtain reliable estimates of the M8 functions
involved in the TIP diagnosis and their limits
of anomalously large values, as approved in course
of the M8 algorithm real-time applications worldwide
over 25 years (Kossobokov 2014). However,
this condition is satisfied just before the 2011 Van
earthquake. Following the practice of the first retrospective
applications of theM8 algorithm (Keilis-
Borok and Kossobokov 1990), we have checked
both a possibility of error in selecting To as
well as a stability of the modified M8 algorithm
applications by simulation of ‘seismic history’ with
Te changing from 1997 to 2011, To = 1975, and
Masoud Mojarab et al.
Figure 14. The results of the M8 algorithm application targeting M7.0+ earthquake with 281 km radius of investigation, as
of July 2011: uniform grid. Multiple intersection of CI’s with TIP’s. The centers of CIs with TIP’s are marked green. The
number of earthquakes in CI’s centered at white points is insufficient for running the modified version of the M8 algorithm.
Tb = 1987. The results of these experiments
confirm reliability of the TIP diagnostics and its
gradual improvement, which appears to stabilize after
2007. Obviously, such a short background of retrospective
monitoring in a relatively small region
does not allow assessing reliable estimates of statistical
significance of the M8 application. However,
the 2011 Van earthquake is predictable in all
the three approaches, and applications guided by
location of D-intersections, and density distribution
peaks appear more adequate than that on a
regular grid.
6. Conclusion
The algorithm used in this study is based on
the simultaneous quantitative analysis of a set of
precursory seismicity patterns (Keilis-Borok and
Kossobokov 1990). These integral seismic symptoms
of general activation of seismic static may be
interpreted as response of the lithosphere to the
tectonic stress. Direct analogies of such a response
exist in some other complex non-linear systems
including but not limited to solar flares, starquakes,
magnetic disturbances, outcomes of elections, starts
and ends of economic recessions, episodes of a
sharp increase in the unemployment rate, and surges
of homicides in a mega-city, etc. (Keilis-Borok 1996;
Keilis-Borok et al. 2001; Kossobokov and Soloviev
2008). The on-going global test of M8 has already
demonstrated predictability of the world’s strongest
earthquakes (Kossobokov et al. 1999; Kossobokov
and Shebalin 2003; Kossobokov 2011, 2013, 2014).
However, there are many challenging questions
related to extension of reliable earthquake predictions
to lower, yet, hazardous magnitude ranges,
as well as to reduction of their temporal and
spatial uncertainty.
Present study has demonstrated that the modified
version of M8 algorithm is able to respond efficiently
to the questions of earthquake forecasting in
eastern Anatolia. The alarm area could have been
diagnosed in advance of the 23 October 2011, M7.3
Van earthquake. Apparently, the results of retrospective
analysis suggest a possibility to increase
the accuracy of locating the incipient event by narrowing
down the spatial extent of prediction to
strong earthquake-prone sites or to the intersection
between the alerted areas, which is obviously
smaller than each of them. Our test of the three
different approaches to application of the M8 algorithm
shows that location of centers of CIs plays
a vital role emphasized from the very first studies
(Kossobokov 1986; Keilis-Borok and Kossobokov
1987). According to our test results, distribution
of CIs based on a natural approach accounting for
tectonics of the region (i.e., D-intersections and
extremes of earthquake spatial density) provides
Earthquake prediction algorithm M8 in eastern Anatolia
more significant and efficient predictions than the
one based on a regular grid.
Finally, the present study encourages using the
M8 algorithm, which can be applied with little
modifications, as a reliable tool to monitor the
times of increased probability of earthquake occurrence
all over the Iranian–Turkish plateau as well
as in other seismic regions worldwide (Keilis-Borok
and Kossobokov 1990). We consider the presented
case-study targeting predictability of the 2011 Van
earthquake as an essential part of learning before
its widespread application all over the Iranian
plateau for a retrospective check of the modified
version of M8 algorithm. In fact, such an application
expanded to the entire Iranian Plateau was
completed in advance to the two recent earthquakes
from M7.5+ range on both sides of the
Iran–Pakistan boundary. The results of application
in Iranian Plateau along with the evidence of
timely real-time prediction of the 16 April 2013,
M7.7 Saravan and the 24 September 2013, M7.7
Awaran earthquakes (Mojarab et al. 2014) confirm
reliability of the modification of M8 algorithm
presented here and its potential in a real-time
monitoring of TIPs in Iran and surroundings.
Acknowledgements
The authors appreciate constructive comments
and suggestions of an anonymous reviewer. VK
acknowledges the support from the Russian Foundation
for Basic Research (grants RFBR Nos.
13-05-91167 and 14-05-92691).
Appendix (M8 algorithm)
The M8 algorithm is one of the intermediateterm,
middle-range earthquake prediction methods.
It was designed by retrospective analysis of
the dynamics in seismicity preceding the great,
magnitude 8.0 or larger earthquakes worldwide
(hence its name). Its archetype (Keilis-Borok
and Kossobokov 1984) and the original version
(Kossobokov 1986; Keilis-Borok and Kossobokov
1987; Kossobokov et al. 1997) were tested retrospectively
in the vicinities of the 132 epicenters of
earthquakes of magnitude 8.0 or greater recorded
between 1857 and 1983. An experiment aimed
to test in a real-time prediction mode, the performance
of an intermediate-term, middle-range
earthquake prediction algorithm M8 at a global
scale started in 1992 (Healy et al. 1992). This
global test intended to predict the largest earthquakes
worldwide has been carried on routinely in
real-time for more than 20 years (Kossobokov et al.
1999, 2002; Kossobokov 2013, 2014). The results of
the global test prove predictability of earthquakes
in the M8.0+ and M7.5+ ranges.
A brief description of the M8 algorithm is as
follows:
M8 algorithm is aimed to predict the earthquakes
from magnitude range MM0+ = [M0,M0+
Δm] where M0 is a constant Δm < 1. If the data
are adequately complete, a number of target intervals
M0 ≤ M < (M0 +Δm), denoted hereinafter
as MM0+, could be distinguished for different values
of M0 with an increment of 0.5. In some cases,
the actual distribution of earthquake size suggests
a natural cut-off magnitude that determines characteristic
earthquakes at different levels of seismic
hierarchy.
Overlapping circles of investigation, CIs, of the
fixed size, proportional to a source of a target earthquake,
scan the territory of the seismic region under
study. The diameter D(M0) = exp(M0 − 56) + 1
in degrees of the Earth’s meridian is used for application,
i.e., about 5–10 larger than a source of a
target magnitude MM0+ event.
Within each CI, the sequence of main shocks is
considered with aftershocks removed (Gardner and
Knopoff 1974; Keilis-Borok et al. 1980). Each main
shock is described with a vector {ti, mi, hi, bi(e)},
i = 1, 2, . . . where i is the main shock number, ti
is its origin time, ti ≤ ti+1; mi is its magnitude,
hi is its focal depth, and bi(e) is the number of
aftershocks with magnitude equal or above prefixed
threshold Maft that happen during the first e days.
Sequences in different CIs are normalized to about
the same pre-fixed average annual number of earthquakes
˜N
by selecting the lower magnitude cut-off
M = Mmin( ˜N ).
Features of the seismic sequence in time are
quantified by computing several running averages
in the trailing time window (t−s, t) and magnitude
range M0 >Mi ≥ m. These counts include:
(i) the number of earthquakes N(t) of magnitude
m or greater in time window (t − s, t);
(ii) the deviation of N(t) from longer-term trend,
L(t);
(iii) linear concentration Z(t) estimated as the
ratio of the average source diameter to the
average distance between sources; and
(iv) the maximum number of aftershocks B(t).
Specifically, function N(t) = N(t|m, s) is the number
of earthquakes with M ≥ m in time interval
from (t − s) to t and represents the rate of
seismic activity. Function L(t) = L(tm, s, t0) =
N(t|m, s)−N(t−s|m, t−s−t0)×s/(t−s−t0) represents
deviation of seismic activity from a longerterm
trend over the period from t0 to t. Function
Z(t) = Z(t|m,M0 − g, s,α, β) represents the linear
concentration of earthquake sources, i.e., the
Masoud Mojarab et al.
average source size (n)−1 × Σ10β(Mi−α) of earthquakes
in time interval from (t − s) to t and
magnitude (m ≤ Mi < M0 − g), divided by the
average distance between them ∼ (n)−1/3. (Note
that in general Σ10β(Mi−α) may estimate different
properties of an earthquake depending on the
choice of β, e.g., β ≈ 0.5 for linear dimension, β ≈
1.0 for area, and β ≈ 1.5 for volume or energy
release.) B(t) characterizes clustering by accounting
the maximum of the number of aftershocks
bi(e) for earthquakes with M0 − p ≤ Mi <M0 − q
in time interval from (t − s ) to t.
Each of the functions N,L,Z is calculated twice,
with a different magnitude cut-off M = Mmin( ˜N ),
for ˜N = 20 and 10, respectively. As a result, the
earthquake sequence is represented by a robust
averaged description defined by seven functions:
N,L,Z (twice each), and B. ‘Very large’ values
are identified for each function by using the condition
that they are higher than Q percent of the
encountered values. A TIP is declared for 5 years,
when at least six out of seven functions, including
B, become ‘very large’ within a narrow time window
(t − u, t). To stabilize the prediction, this criterion
is required for two consecutive moments, t
and t + 0.5 years.
The following values of the parameters of N(t),
L(t) and Z(t) functions are fixed a priori in the
original version of algorithm M8 (Keilis-Borok and
Kossobokov 1987): D(M0) = {exp(M0 −5.6)+1}0
in degrees of the Earth meridian, which is 384, 560,
854, and 1333 km for M0 = 6.5, 7.0, 7.5, and 8,
respectively; s = 6 years, s = 1 year, g = 0.5,
p = 2, q = 0.2, u = 3 years, and Q = 75% for
B and 90% for the other six functions. Note that
such a dimension of the areas of investigation are
set proportional to the linear dimensions of the target
earthquake sources and are well within the limits
(about 3–4 times smaller) of their preparation
zone with an estimated radius of R = 100.43M km
(Dobrovolsky et al. 1979). This scaling is independently
confirmed by Bowman et al. (1998), who
claimed that the size of the critical region of accelerated
energy release scales with the magnitude of
the impending event according to log R ∼ 0.44 M.
Clearly, the magnitude scale, which is used in application,
should reflect spatial extent of an earthquake
source. For many catalogs, this magnitude
could be the maximum reported one, e.g., in this
paper the maximum magnitude is selected from the
four values provided by the National Earthquake
Information Center/US Geological Survey Global
Hypocenters’ Data Base System, i.e., average mb,
average MS, and the two authoritative magnitudes
MA1 and MA2. The performance of the M8 algorithm
in the on-going Global Test started in 1992
(Healy et al. 1992) is summed up in table 3.
Glossary
The circles of investigation, CIs: The M8
algorithm examines seismicity in circular regions.
In early testing of the algorithm on magnitude
8+ earthquakes, the best results were obtained
using a radius of six degrees of the Earth meridian,
i.e., 667 km. Instrumental observations had
established an empirical multiplicative scaling law
between magnitude M and source dimensions,
so that Dobrovolsky et al. (1979) estimate the
radius of the preparation zone as 100.43M km, i.e.,
about exp(M). After some testing and accounting
for an error in location of epicenter, the developers
of the M8 algorithm set the relationship
between the radius of the CI and the magnitude
of the predicted earthquake, M0, as: R(M0) = 55.5
(exp(M0−5.6)+1) km; which gives, 667 km for
magnitude 8.0, 427 km for 7.5, 281 km for 7.0 and
192 km for 6.5, since CIs are areas of investigation
in which M8 algorithm will run.
TIPs (time of increased probability): A time
of increased probability is declared for the CI when
five of the first six functions and the seventh function
have voted sometime in the preceding three
years, and when this condition is met in two successive
six-month evaluations. Once a TIP is declared,
it lasts for five years. As new data is added to the
catalog, the thresholds which cause the functions
to vote may change, so a TIP may be terminated
early or extended for more than five years. Then
TIP is a time of increased probability on each CI.
D-intersections: D-intersections are the morphostructural
intersections with the potential of generating
earthquakes of magnitude 6.5+ and 7+.
These intersections were introduced by Gelfand
et al. (1974a, b) as dangerous (D-intersections)
prone to strong earthquakes. According to the morphostructural
investigation (Gelfand et al. 1974a,
b), areas prone to earthquakes of 6.5+ and 7+
magnitude ranges are determined with circles 40
and 62.5 km radius, about one source size of
the target shocks, respectively (Kossobokov and
Rotwain 1977). In this study, we have used these
D-intersactions as the centers of CIs with R(M0) for
running M8 algorithm.
References
Allen C R, Edwards W, Hall W J, Knopoff L, Raleigh C B,
Savit C H, Toksoz M N and Turner R H 1976 Predicting
earthquakes: A scientific and technical evaluation with
implications for society; Panel on Earthquake Prediction
of the Committee on Seismology, Assembly of Mathematical
and Physical Sciences, National Research Council, US
National Academy of Sciences, Washington DC.
Earthquake prediction algorithm M8 in eastern Anatolia
Ambraseys N N 1988 Engineering seismology; Earthquake
Eng. Struct. Dynam. 17 1–105.
Ambraseys N N 2009 Earthquakes in the Mediterranean and
Middle East: A multidisciplinary study of seismicity up
to 1900; Cambridge University Press, UK, 968p.
Berberian F, Muir I D, Pankhurst R J and Berberian M 1982
Late Cretaceous and early Miocene Andean-type plutonic
activity in northern Makran and central Iran; J. Geol.
Soc. London 139 605–614.
Bird P 2003 An updated digital model of plate boundaries;
J. Earth Sci. 4(3) 1027.
Bowman D D, Ouillon G, Sammis C G, Sornette A and
Sornette D 1998 An observational test of the critical
earthquake concept; J. Geophys. Res. 103(24) 359–372.
Dewey J F, Hempton M R, Kidd F W S, Saroglu F and
Sengor A M C 1986 Shortening of continental lithosphere:
The neotectonics of eastern Anatolia, a young collision
zone; In: Collision Tectonics (eds) Coward M P and Ries
A C, Geol. Soc. Spec. Publ. 19 3–36.
Dobrovolsky I R, Zubkov S I and Myachkin V I 1979 Estimation
of the size of earthquake preparation zone; Pure
Appl. Geophys. 117 1025–1044.
Emre ¨ O, Duman T Y, ¨ Ozalp S and Elmacı H 2011 Site observations
and preliminary evaluation of source fault of 23
October 2011 Van earthquake; Active Tectonics Research
Group, MTA Publications, Ankara, Turkey (in Turkish).
Gardner J K and Knopoff L 1974 Is the sequence of earthquakes
in southern California, with aftershocks removed,
Poissonian? Bull. Seismol. Soc. Am. 64(1) 363–367.
Gelfand I M, Guberman Sh A, Zhidkov M P, Kaletzkaya
M S, Keilis-Borok V I, Ranzman E Ia and Rotwain I M
1974a Recognition of places where strong earthquakes
may occur. II: Four regions in the Asia Minor and Southeastern
Europe; In: Computer analysis of digital seismic
data (ed.) Keilis-Borok VI, Comput. Seismol. 7 3–40 (in
Russian).
Gelfand I M, Guberman Sh A, Zhidkov M P, Keilis-Borok
V I, Ranzman E Ia and Rotwain I M 1974b Recognition of
places where strong earthquakes may occur. III: The case
when the boundaries of disjunctive knots are unknown;
In: Computer analysis of digital seismic data (ed.) Keilis-
Borok VI, Comput. Seismol. 7 41–58 (in Russian).
Gulkan P, Gurpinar A and Celebi M 1978 Engineering
Report on the Muradiye-C¸ aldıran Turkey earthquake 24
November 1976; National Academy of Sciences, 59p.
Healy J H, Kossobokov V G and Dewey J W 1992 Test
to evaluate the earthquake prediction algorithm, M8; US
Geol. Surv. Open-File Report 92–401, 23p.
Hough S 2009 Predicting the unpredictable: The tumultuous
science of earthquake prediction; Princeton University
Press, 272p.
Keilis-Borok V I 1996 Intermediate term earthquake prediction;
Proc. Natl. Acad. Sci. USA 93 3748–3755.
Keilis-Borok V I and Kossobokov V G 1984 A complex of
long term precursors for the strongest earthquakes of the
world; Proc. 27th Geological Congress, Nauka, Moscow
61 56–66.
Keilis-Borok V I and Kossobokov V G 1987 Periods of
high probability of occurrence of the world’s strongest
earthquakes; Comput. Seismol. 19 45–53.
Keilis-Borok V I and Kossobokov V G 1990 Premonitory
activation of seismic flow: Algorithm M8; Phys. Earth
Planet. Inter. 61 73–83.
Keilis-Borok V I, Knopoff L, Rotvain I M and Sidorenko T M
1980 Bursts seismicity as long-term precursors of strong
earthquake; J. Geophys. Res. 85 803–811.
Keilis-Borok V I, Ismail-Zadeh A T, Kossobokov V G and
Shebalin P 2001 Non-linear dynamics of the lithosphere
and intermediate-term earthquake prediction; Tectonophys.
338 247–260.
Kossobokov V G 1986 The test of algorithm M8; In:
Algorithms of Long-Term Earthquake Prediction (ed.)
Sadovsky MA, CERESIS Lima, Peru, pp. 42–52.
Kossobokov V G 1997 User manual for M8; In: Algorithms
for Earthquake Statistics and Prediction (eds) Healy J
H, Keilis-Borok V I and Lee W H K, IASPEI Software
Library, vol. 6, Seismol. Soc. Am. El Cerrito, CA.
Kossobokov V G 2011 Are mega earthquakes predictable?
Izvestiya, Atmos. Oceanic Phys. 46(8) 951–961.
Kossobokov V G 2013 Earthquake prediction: 20 years
of global experiment; Nat. Hazards 69 1155–1177, doi:
101007/s11069-012-0198-1.
Kossobokov V 2014 Chapter 18. Times of increased probabilities
for occurrence of catastrophic earthquakes: 25 years
of hypothesis testing in real time; In: Earthquake hazard,
risk, and disasters (eds) Wyss M and Shroder J (London:
Elsevier), pp. 477–504.
Kossobokov V G and Rotwain I M 1977 Recognition of
places where strong earthquakes may occur. VI. Magnitude
M ≥ 7.0; In: Pattern recognition and spectral analysis
in seismology (ed.) Keilis-Borok VI, Comput. Seismol.
10 3–13 (in Russian).
Kossobokov V G and Shebalin P 2003 Earthquake prediction;
In: Nonlinear Dynamics of the Lithosphere
and Earthquake Prediction (eds) Keilis-Borok V I and
Soloviev A A, Springer, Heidelberg, pp. 141–207.
Kossobokov V G and Soloviev A A 2008 Prediction of
extreme events: Fundamentals and prerequisites of verification;
Russian J. Earth Sci. 10 ES2005, doi: 102205/
2007ES000251.
Kossobokov V G, Healy J H and Dewey J W 1997 Testing
an earthquake prediction algorithm; Pure Appl. Geophys.
149 219–232.
Kossobokov V G, Romashkova L L, Keilis-Borok V I and
Healy J H 1999 Testing earthquake prediction algorithms:
Statistically significant real-time prediction of
the largest earthquakes in the Circum-Pacific, 1992–1997;
Phys. Earth Planet. Inter. 111(3–4) 187–196.
Kossobokov V G, Romashkova L L, Panza G F and
Peresan A 2002 Stabilizing intermediate-term middlerange
earthquake predictions; J. Seismol. Earthq. Eng. 8
11–19.
McClusky S et al. 2000 Global Positioning System constraints
on plate kinematics and dynamics in the eastern
Mediterranean and Caucasus; J. Geophys. Res. 105
5695–5719.
McKenzie D P 1972 Active tectonics of the Mediterranean
region; Geophys. J. Astron. Soc. 30 109–185.
Mojarab M, Memarian H, Zare M and Kossobokov V 2014
Adjusting the earthquake prediction algorithm M8 for
application in Iranian Plateau with special reference to
the 16 April 2013, M7.7 Saravan and the 24 September
2013, M7.7 Awaran earthquakes (Manuscript).
Peresan A, Kossobokov V G, Romashkova L and Panza
G F 2005 Intermediate-term middle-range earthquake
predictions in Italy: A review; Earth-Sci. Rev. 69
97–132.
Reilinger R S, McClusky P, Vernant S, Lawrence S, Ergintav
R, Cakmak H, Ozener F, Kadirov I, Guliev R, Stepanyan
M, Nadariya G, Hahubia S, Mahmoud K, Sakr A,
ArRajehi D, Paradissis A, Al-Aydrus M, Prilepin T,
Guseva E, Evren A, Dmitrotsa S V, Filikov F, Gomez R
and Al-Ghazzi K G 2006 GPS constraints on continental
deformation in the Africa–Arabia–Eurasia continental
collision zone and implications for the dynamics of plate
interactions; J. Geophys. Res. 111 B05411, doi: 101029/
2005JB004051.
Seng¨or A M C 1990 A new model for the late Palaeozoic–
Mesozoic tectonic evolution of Iran and implications for
Masoud Mojarab et al.
Oman; In: Geology and tectonics of the Oman Region,
Geol. Soc. London, Spec. Publ. 49 797–831.
Seng¨or A M C and Kidd W S F 1979 Post-collisional tectonics
of the Turkish–Iranian plateau and a comparison
with Tibet; Tectonophys. 55 361–376.
Shebalin P N 1992 Automatic duplicate identification in
set of earthquake catalogues merged together; Open-File
Rep., US Geol. Surv., pp. 92–401 (Appendix II).
Silverman B W 1986 Density Estimation for Statistics and
Data Analysis; Chapman and Hall, New York.
Taymaz T, Jackson J A and McKenzie D 1991 Active tectonics
of the north and central Aegean Sea; Geophys. J. Int.
106 403–490.
Tchalenko J S and Berberian M 1974 The Salmas (Iran)
earthquake of May 6, 1930; Annali di Geofisica 27
1–2.
Zare M, Haghshenas E and Bastami M 2011 The Van,
Turkey Earthquake of 23 October 2011; A Preliminary
Report on the Reconnaissance Visit preformed by IIEES
Reconnaissance Team.
MS received 9 July 2014; revised 3 February 2015; accepted 12 February 2015