Layout 6


1

ANNALS OF GEOPHYSICS, 64, 4, SE433, 2021; doi:10.4401/ag-8566 
O P E N  ACC E S S

Spatial and time variations of seismicity before 
strong earthquakes in the southern part of  
the Balkans 
Emil Oynakov, Emil Botev* 

Department of Seismology and Seismic Engineering National Institute of Geophysics, Geodesy and Geography 
Bulgarian Academy of Sciences (NIGGG-BAS) 

Article history: received September 1, 2020; accepted May 10, 2021 

Abstract  

A retrospective analysis of the spatial and time variations of three main statistical parameters of 
the seismicity before recent 4 stronger earthquakes (2015 – 2020) in the southern Balkans is 
presented. The modern extended software package ZMAP with various advanced seismological 
functions for earthquake catalog analysis is used for estimating the spatial-time variations in a-
value (seismic activity), b-value (slope of the recurrence graph) and z-value (parameter of the 
relative seismic quiescence). The catalog data from constantly updated catalog of the University 
of Athens for the period 1964-2020 and spatial window 32° - 44° N and 10° – 30° E are used for 
the various statistical interpretations. The main result of the whole analysis is that the abnormally 
low b-values and high z-values, defining the zones of relatively seismic quiescence, may be an 
indicator of imminent release of more significant stress in areas adjacent to the zones of relatively 
high a-values. Thus, the result of the proposed joint interpretation of the spatial-time variations 
of these three statistical parameters of seismicity could be considered as a kind of predictor of the 
stronger recent seismic events in the southern part of Balkans. 

Keywords: Earthquake catalog; b-value; a-value; z-value; Relatively seismic quiescence. 

1. Introduction

The territory of the Balkan Peninsula is characterized by active geodynamics. This is the most geodynamical
active region in Central and Eastern Europe. A number of dangerous geodynamic processes of endogenous 
(earthquakes, modern crustal movements, mud volcanoes) and exogenous origin (natural and man-made) are 
observed in the area. Today’s geodynamics of the Balkan region is controlled by the active tectonic processes 
in the Eastern Mediterranean: the collision of the Adriatic (Apulian) microplate with the Dinarides; the 
subduction of the oceanic Ionian and Levantine lithosphere under the Aegean arc system; the collision between 
the Eurasian and Arabian plates and the subsequent displacement of the Anatolian microplate to the west along 
the northern Anatolian fault. 

The main goal of this study is to map and analyze the spatial and time variations of three main statistical 
parameters of seismicity before 4 stronger earthquakes (Ml > 6.0) that occurred in the southern Balkans in the 

annagrazia
Cross-Out



last decade (2015-2020). To quantify the variations of the three parameters a-value (seismic activity), b-value 
(slope of the recurrence graph) and z-value (seismic quiescence), catalog seismic events from 1964 to 2020 are 
used. There are many publications considered geographical distribution of b-value (sometimes and a-value) in 
the seismic zones of southern Balkans [for instance: Hatzidimitriou et al., 1985; Papazachos, 1999; Chouliaras, 
2009; D’Alessandro et al., 2011; Popandopoulos and Baskoutas, 2011; Popandopoulos et al., 2016; Vamvakaris 
et al., 2016]. Some regularities in the time variations of b-values before one weak and one moderate earthquakes 
in two local zones in Greece are presented by Popandopoulos and Baskoutas, [2011]. The authors suggest the 
obtained features of the time series of b to be used for the expected time prediction of the local main 
earthquakes. 

In our study an alternative approach to calculate the regional spatial-time variations of the statistical 
parameters of seismicity (b-, a-, z-values) by means of the extended software package ZMAP [Westerhaus et al., 
2002] is used. The software package is open source software used to assess various seismic hazards [Westerhaus 
et al., 2002; Wiemer and Wyss, 1994, 2000]. ZMAP contains many advanced seismological functions for 
earthquake catalog analysis, including for estimating variations in b-, a-, z-values, time series analysis, 
temperature reduction, and more. As an example of the previous successful use of the capabilities of ZMAP for 
studying the spatial-time variations of the b-value before strong earthquakes, we can mention Damanik’s study 
[Damanik et al., 2010] for the region of the island of Java (106–115.5°, –0, 7 to –11.5°), which uses data from the 
US Geological Survey (USGS) earthquake catalog and the Engdahl catalog for very strong earthquakes (Мmin = 
4.7). We are not familiar with any examples of a combined study of the spatial-time variations of the above 3 
statistical parameters of seismicity (b-, a-, z-values) as predictors of nowadays strong seismic events in the 
southern Balkans. 

 
 

2. Preparation of catalog data  
 
The assessment of spatial-time variations in seismic activity in a given region is important for understanding 

the processes in the seismic-tectonic environment development, as well as for the seismic analysis of the danger 
of an imminent natural disaster. To quantify seismic variations in an area, we need information in the form of 
an earthquake catalog [Burton, 1990]. 

In this study a constantly updated catalog of the University of Athens for the period 1964-2020 is used 
(http://dggsl.geol.uoa.gr/en_index.html). The catalog of earthquakes covers data in a spatial window 32° - 44° 

N and 10° – 30° E. The total number of included earthquakes is 295029, as the depths vary from 0 ≤ h ≤ 220 km, 
and the magnitude estimates Ml are in the range of 2.0 < Мl < 7.0.  

Duplicate events were recognized and removed algorithmically, and later checked by visual assessment. The 
catalog of independent events is defined as final after the identification and removal of clustered events by the 
Gardner and Knopoff algorithm [Gardner and Knopoff, 1974]. Any earthquake that occurred within the spatial-
time window around and after each larger event is taken into account and is considered as a cluster event. The 
spatial-time window is larger in size for stronger previous events according to Gardnerand Knopoff dependency. 
This step is relevant merely for our catalog data comprising large events with long aftershock and 
foreshock,series. 

Figure 1 shows the distribution of the epicenters of the earthquakes after declustering. To guarantee the 
completeness of data, analysis comprise only events with magnitudes equal to or larger than the threshold 
magnitude Mc. So, the next step has been the determination of minimum value of magnitude Mc for catalog 
completeness. The value of Mc usually decreases with time, as the number of seismic stations and their sensitivity 
increase [Wiemer and Wyss, 2000]. The method of maximum curvature (Mmax curvature) is used to verify the 
estimates of Mc [Wiemer, et al., 1999; Wyss et al., 2004] (Figure 2 a, b). The very final catalog includes 24750 events 
with Mc ≥ 3.5 limited mainly to a depth of 40 km (shallow earthquakes) – maximum depth is 70 km (Figure 3). 

 

Emil Oynakov et al.

2



3

Seismic variations before some Balkan events

Figure 1. Map of epicenters of the independent seismic events obtained after declustering of the permanently updated 
earthquakes catalog data of the University of Athens. http://dggsl.geol.uoa.gr/en_index.html

Figure 2. a) Estimation of the minimum magnitude of catalog completeness by the method of Maximum curvature. The 
arrows indicate the breakpoints of the cumulative distribution, marking the minimum magnitude Mc for 
completeness of the catalog; b) spatial distribution of the magnitude threshold Mc. 

http://dggsl.geol.uoa.gr/en_index.html


Emil Oynakov et al.

4

3. Analysis of the b-value 
 
Analysis of the b-value in the studied region was performed before 4 stronger earthquakes with magnitude  

Ml > 6.0. These are the earthquakes from: 16.04.2015 with coordinates 26.82° E, 35.23° N; M1 = 6.1, h = 37 km and 
T0 = 18:07; 12.06.2017 with coordinates 26.36° E, 38.84° N; M1 = 6.1, h = 28 km and T0 = 12:28; 20.07.2017 with 
coordinates 26.56° E, 35.64° N; Ml = 6.2, h = 63 km and T0 = 22: 31 and 20.05.2020 with coordinates 25.69° E, 34.40° 
N; M1 = 6.0, h = 51 km and T0 = 12:51. A polygon locked in a spatial window with coordinates 24.20° - 29.98° E, 
34.40° - 41.00°  N presented in Figure 4, and containing 2637 events with Mc > 4 (Figure 4, b). 

The value of b for the tested region is calculated using the maximum likelihood method [Aki, 1965] (Figure 5), 
then mapped before each of the main studied events using the usual technique with the grid with a cell length of 5 
km and a fixed number of earthquakes (Ni = 50) [Wiemer and Wyss, 2000], as shown in figs. 6a, b, c, d. In the graph 
of the value of b (Figure 5) we can distinguish four periods: 1964-1987 average value of b ≈ 1.42; 1987-2003 average 
value of b ≈ 1.9; 2003-2009 average value of b≈1.4 and 2009-2020 average value of b ≈ 1.0. The last period, which 
includes the studied events, has the highest average value of the b-parameter, which implies increased heterogeneity 
and reduction of stresses [Görgün et al., 2009]. The significant decrease in the value of b can be associated with the 
increasing effective level of stress before major earthquakes [Schorlemmer et al., 2004; Wu, 2008]. 

The spatial variations of b-value for the studied region are estimated for the period as follows: from 1964 to 
01.03.2015 one month before the earthquake on 16.04.2015 (Ml 6.1) (Figure 6, a); from 1964 to 01.05.2017 one 
month before the earthquake on 12.06.2017 (Ml 6.1); from 1964 to 01.07.2017 one month before the earthquake on 
20.07.2017 (Ml 6.2) and from 1964 to 01.04.2020 one month before the earthquake on 02.05.2020 (Ml 6.0). In the 
figures (Figure 6, a, b, c, d) in ascending order are marked from low to high values of the b-parameter. The spatial 
differences in the value of b illustrate variability in plan, as relatively low values can determine the places where 
an earthquake would most likely occur [Schorlemmer et al., 2003; Westerhaus et al., 2002]. The zones with relatively 
low values of b (0.9 - 1.2) are clearly delineated and the epicenter of the earthquake falls into them. According to 
Wyss et al. [2004] and Motaghi et al. [2010] low values of b indicate that stresses accumulate in these zones until 
the activation of the main event. We note that the zones with a low value of b in which the epicenters fall are very 
close to zones with relatively high values. 

 
 

Figure 3. Histogram of the depth distribution of the earthquake hypocenters. 



5

Seismic variations before some Balkan events

Figure 4. a) Epicentral map for the tested polygon; b) estimation of the minimum magnitude Mc for completeness of the 
catalog for the selected polygon.

Figure 5. Time variations of the b-value of earthquakes with Ml > 4 within the studied polygon. The moments of the 
studied earthquakes from 2015 and 2017 are marked with red squares.



4. Analysis of the a-value 
 
The spatial distribution of the a-value characterizing the seismic activity before the studied strong earthquakes 

shows that the epicenters of strong earthquakes fall in areas with relatively low values of seismic activity (a = 6-7) 
(Figure 7 a, b, c, d). The estimates were performed for the same periods before the earthquakes as in the b-value studies. 

Without exception, however, the epicenters of strong earthquakes are located in close proximity to areas with 
relatively high activity (a = 10-12). This result is consistent with the general idea of creating non-equilibrium 
stresses in the vicinity of the zones of high seismic activity and their subsequent redistribution with migration of 
stronger events in the direction of the low active seismic zones. 

Emil Oynakov et al.

6

Figure 6. Spatial distribution of the b- value before the earthquake: a) on 16.04.2015 (Ml 6.1); b) on 12.06.2017 (Ml 6.1); 
c) on 20.07.2017 (Ml 6.2); d) on 02.05.2020 (Ml 6.0);  - epicenters of earthquakes.



 
 

5. Analysis of the z-value 
 
The method of the z-test of seismic catalogs for determining the spatial-time areas of seismic quiescence before 

strong earthquakes is based on the research of Wyss and Habermann and [Wiemer and Wyss, 2002; Habermann, 
1987]. This method is focused on the determination of spatial-time blocks in the seismically active zone with a 
significant change in the intensity (rate) of the seismic activity in the selected energy range. The analysis is based 
on the statistical function of the z-test. To determine the seismic quiescence, the study area is covered with a fixed 
node grid. For each network node at a given moment, the function z(t) is calculated 

      z(t) =  ,  [1] 
 
 

where R1 and R2 are respectively the average values of the velocity of earthquakes flow (number of earthquakes in 

R₁ – R₂
σ₁2  +σ₂2  n₁    n₂�
———

7

Seismic variations before some Balkan events

Figure 7. Spatial distribution of the a-value before the earthquake: a) on 16.04.2015 (Ml 6.1); b) on 12.06.2017 (Ml 6.1); c) 
on 20.07.2017 (Ml 6.2); d) on 02.05.2020 (Ml 6.0);  - epicenters of earthquakes.



a time window 30 days) for two-time intervals (one with a duration of not less than a year and tied to the test date, 
the other - includes all other time intervals); σ1 and σ2 - standard deviations of R1 and R2 respectively for the first 
and second time interval, and n1 and n2 - number of earthquakes for the respective period. Earthquakes with a 
hypocenter depth of up to 70 km are considered for each node of the network. High (positive) values of z indicate 
a decrease in the rate of the seismic activity, and respectively low z-values indicate increase [Maeda K., Wiemer S., 
1999]. All nodes of the network with z > 3 correspond to 99% reliability of determination of seismic quiescence 
[Wiemer, 2001; Wiemer and Wyss, 1994]. The calculated z-values in the nodes of the network are combined on the 
principle of spatial-time adjacency and determine the areas with seismic quiescence. The spatial distribution of 
the z-value before the earthquakes (Figure 8 a, b, c) are calculated for the same polygon, comparing two time periods 
as follows: t1 = 30.03.2012-30.09.2013 / t2 = 30.09. 2013-30.03.2015 for the earthquake on 16.04.2015 (Ml 6.1) 
(Figure 6, a); t1 = 30.05.2014-30.11.2015 / t2 = 30.11.2015-30.05.2017 for the earthquakes on 12.06.2017 (Ml 6.1), 
and on 20.07.2017(Ml 6.2) (Figure 6, b) and t1 = 28.02.2017-30.08.2018 / t2 = 30.08.2018-30.04.2020 for the 
earthquake on 02.05.2020 (Ml 6.0), (Figure 6, c). 

The high (positive) z-values of the maps can be interpreted as a decrease in the rate of seismic activity (seismic 
quiescence) compared to the first period, and the low (negative) z-values represent an increase in the seismic rate. 
Earthquake density and distribution is a critical factor in interpreting z-value variations. Large areas of constant 
color (value) could show the same density of earthquakes for different periods of time, i.e., may show a homogeneous 
degree of seismicity in this area. 

From the figures below (Figure 8 a, b, c) it is observed that the epicenters of the stronger events fall in areas 
with relatively high values of the z-value (4-5), which means that the selected periods before the earthquakes are 
periods of relative seismic quiescence. The relatively high values of z = 4 – 5 > 3 show 99% reliability of the results. 

 
 

6. Conclusions 
 
The time variations of the b-value for the period 1964 - 2020 shows a minimum of the average value of b for the 

period 2013-2020 (b ≈1), which includes the four studied stronger earthquakes. 
Large decreases in the value of b may be associated with increasing effective levels of stress before major 

earthquakes. The increase in the b-value after these earthquakes may indicate an increase in the heterogeneity of 
the earth’s crust and a decrease in shear stress, but after the study period the value of b does not increase 
significantly which may indicate imminent strong seismic events in the study area. 

The change in the spatial distribution of the b-value before the stronger earthquakes shows that the area with 
an abnormally low value of b covers the epicenters of the studied stronger earthquakes. These low values of b can 
be interpreted as potentially locked areas or areas with accumulated high stress before major earthquakes. 

The spatial distribution of the α-value characterizing the seismic activity before the studied stronger earthquakes 
shows that the epicenters of these earthquakes fall in areas that have relatively low values of seismic activity but 
are close to areas with relatively high activity.  

The results of the conducted research show that the epicenters of the selected stronger earthquakes are located 
in zones of relatively high z-value (z ≈ 4.2 > 3). This fact definitely indicates a statistically reliable marking of an 
area with a relatively seismic quiescence before the corresponding stronger earthquake. 

As a result of the whole analysis we can summarize that the abnormally low b-values and high z-values, defining 
the zones of relatively seismic quiescence, may be an indicator of imminent release of more significant stress in areas 
adjacent to the zones of relatively high a-values. Thus, the joint interpretation of the spatial-time variations of 
these three statistical parameters of seismicity can be interpreted as predictor of strong seismic events and also 
could be used for the assessment of seismic hazard in the southern part of the Balkan. 
 
 
Acknowledgments. This study is supported by the Bulgarian National Science Fund, Contract DN 14-1/11.12.2017 
“Exploration of changes in some geophysical fields preceding the occurrence of earthquakes in the Balkans” and 
particularly supported by the Ministry of Education and Science (MES) of Bulgaria (Agreement No Д01-322/18.12.2019). 
 
 

Emil Oynakov et al.

8



9

Seismic variations before some Balkan events

References 
 
Aki K. (1965). Maximum likelihood estimate of b in the formula log N = a-bM and its confidence limits, Bull. Earthq. 

Res. Inst., Tokyo Univ., 43, 237-239. 
Chouliaras G. (2009). Investigating the earthquake catalog of the National Observatory of Athens, Nat. Hazards 

Earth Syst. Sci., 9, 905–912, https://doi.org/10.5194/nhess-9-905-2009. 
Burton P.W. (1990). Pathways to seismic hazard evaluation: Extreme and characteristic earthquakes in areas of low 

and high seismicity, Nat. Hazards, 3, 275–291. 
D’Alessandro A., D. Papanastassiou, I. Baskoutas (2011). Hellenic Unified Seismological Network: an evaluation of 

its performance through SNES method, Geophys J. Int., 185, 1417–1430. 
Damanik R., H.E. Andriansyah Putra, M.T. Zen (2010). Variations of b-values in the Indian Ocean – Australian plate 

subduction in south Java Sea; In: Proceedings of the Bali 2010 International Geosciences Conference and 
Exposition, Bali, Indonesia, 19–22.  

Figure 8. Spatial distribution of the z- value (parameter of the relative seismic quiescence) before the earthquake: a) on 
16.04.2015 (Ml 6.1); b) on 12.06.2017 (Ml 6.1) and on 20.07.2017 (Ml 6.2); c) on 02.05.2020 (Ml 6.0);  - epicenters 
of earthquakes.



Gardner J.K. and L. Knopoff (1974). Is the sequence of earthquakes in southern California, with aftershocks removed, 
Poissonian? Bull. Seism. Soc. Amer. 64, 5, 1363-1367. 

Görgün E., A. Zang, M.Bohnhoff, C. Milkereit, G. Dresen (2009). Analysis of Izmit aftershocks 25 days before the 
November 12th 1999 Düzce earthquake, Turkey, Tectonophysics, 474, 3-4, 507-515. 

Habermann, R.E. (1987). Man-made changes of seismicity rates, Bull. Seism. Soc. Am., 77, 141–159.  
Hatzidimitriou P.M., E.E. Papadimitriou, D.M. Mountrakis, B.C. Papazachos (1985). The seismic parameter b of the 

frequency–magnitude relation and its association with the geological zones in the area of Greece, 
Tectonophysics, 120, 141–151. 

Maeda K., S. Wiemer (1999). Significance test for seismicity rate changes before the 1987 Chiba-toho-oki earthquake 
(M6.7), Japan, Ann. Geofis, 42, 5, 833–850.  

Motaghi K., K. Hessami, M. Tatar (2010). Pattern recognition of major asperities using local recurrence time in 
Alborz Mountains, Northern Iran, J. Seismol., 14, 787–802. doi:10.1007/s10950-0109201-z. 

Papazachos C. (1999). An alternative method for a reliable estimation of seismicity with an application in Greece 
and the surrounding area, Bull. Seismol. Soc. Am., 89, 1, 111–119. 

Popandopoulos G.A., I. Baskoutas (2011). Regularities in the time variations of seismic parameters and their 
implications for prediction of strong earthquakes in Greece, Izv., Phys. Solid Earth, 47, 11, 974–994. 

Popandopoulos, G.A., I. Baskoutas, E. Chatziioannou (2016). The spatiotemporal analysis of the minimum magnitude 
of completeness Mc and the Gutenberg–Richter law b-value parameter using the earthquake catalog of Greece. 
Izv., Phys. Solid Earth, 52, 195–209. 

Schorlemmer D., G. Neri, S. Wiemer, A. Mostaccio (2003). Stability and significance tests for b-value anomalies: 
example from the Tyrrhenian Sea, Geophys. Res. Lett., 30, 16, 1835, doi:10.1029/ 2003GL017335. 

Schorlemmer D., S. Wiemer, M. Wyss (2004). Earthquake statistics at Parkfield, Stationarity of b values, J. Geophys. 
Res., 109, B12307, doi:10.1029/2004JB003234. 

Stiphout T., D. Schorlemmer, S. Wiemer (2011). The effect of uncertainties on estimates of background seismicity 
rate, Bull. Seism. Soc. Am., 101, 2, 482–494, doi:10.1785/0120090143. 

University of Athens, Catalog of earthquakes data, http://dggsl.geol.uoa.gr/en_index.html.  
Vamvakaris D.A., C.B. Papazachos, Ch.A. Papaioannou, E.M. Scordilis, G.F. Karakaisis (2016). A detailed seismic 

zonation model for shallow earthquakes in the broader Aegean area, Nat. Hazards Earth Syst. Sci., 16, 55–84, 
doi:10.5194/nhess-16-55-2016. 

Westerhaus M., M. Wyss, R. Yilmaz, J. Zschau (2002). Correlating variations of b values and crustal deformation 
during the 1990’s may have pinpointed the rupture initiation of the Mw = 7.4 Izmit earthquake of 1999 August 
17, Geophys. J. Int.,148, 139–152.  

Wiemer S. (2001). A program to analyse seismicity: ZMAP, Geophys. Res. Lett., 72, 373–382.  
Wiemer S., K. Katsumata, (1999). Spatial variability of seismicity parameters in aftershock zones, J. Geophys. Res., 

104, 13, 135–151. 
Wiemer S., M. Wyss (1994) Seismic quiescence before the landers (M = 7.5) and big bear (M = 6.5), 1992 earthquakes, 

Bull. Seism. Soc. Am., 84, 3, 900–916. 
Wiemer S., M. Wyss (1997). Mapping the frequency-magnitude distribution in asperities: an improved technique to 

calculate recurrence times? J. Geophys. Res., 102, 15, 115–128. 
Wiemer S., M. Wyss (2000). Minimum magnitude of completeness in earthquake catalogs: examples from Alaska, 

the western United States, and Japan, Bull. Seism. Soc. Am., 90, 4, 859–869. 
Wiemer S., M. Wyss (2002). Mapping spatial variability of the frequency–magnitude distribution of earthquakes, 

Adv. Geophys., 45, 259–302. 
Wyss M., R.E. Habermann (1988). Precursory seismic quiescence, Pure Appl. Geophys., 126, 319–332.  
Wyss M., D. Schorlemmer, S. Wiemer (2000). Mapping asperities by minima of local recurrence time: the San Jacinto-

Elsinore fault zones, J. Geophys. Res., 105, B4, 7829–7844. 
Wyss M, R. Stefansson (2006). Nucleation points of recent main shocks in southern Iceland mapped by b-values, Bull. 

Seism. Soc. Am., 96, 599–608, doi:10.1785/0120040056.  
Wyss M., G. Sobolev, J.D. Clippard (2004). Seismic quiescence precursors to two M7 earthquakes on Sakhalin Island, 

measured by two methods, Earth Planets Space, 56, 725–740. 
Wyss M., F. Pacchiani, A. Deschamps, G. Patau (2008). Mean magnitude variations of earthquakes as a function of 

depth: different crustal stress distribution depending on tectonic setting, Geophys. Res. Lett., 35, L01307, 

Emil Oynakov et al.

10



11

Seismic variations before some Balkan events

doi:10.1029/2007GL031057. 
Wu Y.M., L.Y. Chiao (2006). Seismic quiescence before the 1999 Chi-Chi, Taiwan Mw 7.6 earthquake, Bull. Seism. 

Soc. Am., 96, 321 – 327.  
Wu Y.M., C.H. Chang, L. Zhao, T.L. Teng, M. Nakamura (2008). A comprehensive relocation of earthquakes in Taiwan 

from 1991 to 2005, Bull. Seism. Soc. Am., 98, 3, 1471-1481. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 
*CO R R E S PO N D I N G A U T H O R: Emil B OT E V,  

Department of Seismology and Seismic Engineering 

National Institute of Geophysics, Geodesy and Geography,  

Bulgarian Academy of Sciences, (NIGGG-BAS),  

e-mail: ebotev@geophys.bas.bg 

© 2021 the Author(s). All rights reserved. 

Open Access. This article is licensed under a Creative Commons Attribution 3.0 International