Pure Appl. Geophys. Ó 2017 Springer International Publishing AG DOI 10.1007/s00024-017-1617-8

Pure and Applied Geophysics

Analysis of the 2005–2016 Earthquake Sequence in Northern Iran Using the Visibility Graph Method NAEEM KHOSHNEVIS,1

RICARDO TABORDA,1,2 SHIMA AZIZZADEH-ROODPISH,1,2 and LUCIANO TELESCA3

Abstract—We present an analysis of the seismicity of northern Iran in the period between 2005 and 2016 using a recently introduced method based on concepts of graph theory. The method relies on the inter-event visibility defined in terms of a connectivity degree parameter, k, which is correlated with the earthquake magnitude, M. Previous studies show that the slope m of the line fitting the k–M plot by the least squares method also observes a relationship with the b value from the Gutenberg–Richter law, thus rendering the graph analysis useful to examine the seismicity of a region. These correlations seem to hold for the analysis of relatively small sequences of earthquakes, offering the possibility of studying seismicity parameters in time. We apply this approach to the case of the seismicity of northern Iran, using an earthquake catalog for the tectonic seismic regions of Azerbaijan, Alborz, and Kopeh Dagh. We use results drawn for this region with the visibility graph approach in combination with results from other similar studies to further improve the universal relationship between m and b, and show that the visibility graph approach can be considered as a valid alternative for analyzing regional seismicity properties and earthquake sequences. Key words: Visibility graph, seismicity, northern Iran.

1. Introduction In a relatively recent study, Lacasa et al. (2008) introduced a simple and computationally inexpensive approach to convert time-series into a particular kind of mathematical graphs known as visibility graphs. In this approach, the connection between the nodes in a network—here understood as events in a time series—is based on the visibility of each node or event

1 Center for Earthquake Research and Information, The University of Memphis, Memphis, TN 38152, USA. E-mail: [email protected] 2 Department of Civil Engineering, The University of Memphis, Memphis, TN 38152, USA. 3 Institute of Methodologies for Environmental Analysis, National Research Council, Tito, Italy.

with respect to other previous and later events in a time sequence. Over the last decade, this method has been found to be applicable to multiple fields, including economics, medicine, and geophysics, using both directed and undirected, in one- and twodirectional approaches (e.g., Yang et al. 2009; Elsner et al. 2009; Telesca et al. 2012; Wang et al. 2012; Long 2013). For the particular case of seismological applications, previous studies have proposed to consider the nodes in the graph to be the earthquakes in a given time series of events, and analyzed the properties of the graph with respect to traditional seismicity parameters. The analysis of seismic sequences through the visibility graph approach for various tectonic seismic regions has proven to be a valid alternative to studying magnitude time series. Telesca and Lovallo (2012), for instance, studied the seismicity of Italy between the years 2005 and 2010 using the visibility graph method. Applying different threshold magnitudes to construct the graph and observing the collapsing effect of all distribution degrees, they observed that the properties of the visibility graph seemed to depend only on the magnitude values, and not on the threshold magnitude used in the analysis. Subsequently, Telesca et al. (2013) studied the seismicity of the Mexican subduction zone through the visibility graph approach and found that the properties of the graph were correlated to the seismic b value in the Gutenberg–Richter law (Gutenberg and Richter 1944). In particular, Telesca et al. (2013) extracted the characteristics of the visibility graph for five different tectonic seismic regions in the Mexican subduction zone and found that the slope m in the linear relationship between earthquake magnitude, M,

N. Khoshnevis et al.

and the inter-event visibility in the graph, defined as the connectivity degree parameter, k, was correlated with b. In a similar article, Telesca et al. (2014b) studied the seismicity of the 2002–2011 Pannonian region using the visibility graph method for two sub-catalogs of shallow and deep earthquakes. They extracted the visibility graph characteristics for each group of events and confirmed that there was a close correlation between the Gutenberg–Richter b value and the slope m of the k–M relationship obtained from the visibility graph. According to Telesca et al. (2014b), the high linear correlation coefficient value (close to 1.0) between the b and m indicates that this relationship exhibits a nearly universal character. This observation was reinforced experimentally by Telesca et al. (2014a), who investigated the behavior of a mechanical stick-slip system with different asperities (sandpaper of different grades). After collecting data for different experiments emulating young and mature faults, Telesca et al. (2014a) found again that b was linearly correlated with the value of m derived from the visibility graph of the experimental events. More recently, Telesca et al. (2016) used the visibility graph method to gain insight about the timespan between earthquakes. The authors defined the parameter \Tc[ as the window mean interval connectivity time, which provides information about the average of all the time intervals between visible events. Applying this concept to the 2003–2012 earthquake sequence in Kachchh, western India, Telesca et al. (2016) found that the variation of\Tc[in time exhibited a plausible relationship with the occurrence of earthquakes. In particular, they observed that the value of \Tc[ decreases significantly before the largest shock in the Kachchh catalog. In this article we use the visibility graph method to analyze the seismicity of northern Iran from 2005 to 2016. We focus our analysis on these years because in a previous study we observed that the latest decade of seismic records in the region offered the most complete catalog of events Khoshnevis et al. (2017). Upon a brief review of the basic concepts of the method and its applications, we describe the seismicity of northern Iran considering three dominant seismic regions: Azerbaijan, Alborz, and Kopeh Dagh. These are regional areas of considerable

Pure Appl. Geophys.

ground shaking activity, with average seismic rates varying between about 0.17 and 3.25 events per year (e.g., Nemati 2015). The considerable amount of earthquakes registered in a 11-year period facilitates the analysis done with the visibility graph approach within what would otherwise be considered a relatively short period of time, and serves as a good case study to test the capabilities of the method. We first summarize the dataset, including a description of the seismic catalog and its completeness based on a recent work on the seismicity of the region (Khoshnevis et al. 2017). We then present results obtained for the region regarding the relationship between the graph’s k–M slope and the seismicity b value for the three seismic areas, analyze the sensitivity of the results to the definition of the catalogs, the number of events and composition, and present results for the variation of \Tc[ over time.

2. Methodology Discrete mathematics and computer science often deal with different types of data abstractions and structures such as sets, trees, and graphs, where concepts like nodes and links are used to describe the topology of a collection of objects and manipulate their data (e.g., Skiena 2008). In graphs, the data— sometimes also referred to as the payload—resides at the vertices or nodes who are connected through links or edges, which may or may not have an associated direction. A visibility graph is a special type of graphs in which the links are straight lines connecting intervisible nodes; that is, straight lines that do not go across any obstacle while connecting nodes that can see each other in a physical space (Lozano-Pe´rez and Wesley 1979). Visibility graphs have been mostly used in robotics for navigation path planning (e.g., Huang and Chung 2004; Oommen et al. 1987) but have also seen applications in other fields including urban studies, interior architecture, medicine, and geosciences (e.g., Raman 2010; Ahmadlou et al. 2010; Varoudis and Psarra 2014; Phillips et al. 2015). In a relatively recent study, Lacasa et al. (2008) applied the concepts of visibility graphs to the representation and analysis of time series. Multiple applications have been found to this idea in fields

Analysis of the 2005–2016 earthquake sequence

yi  yp yi  yj [ ; tp  ti tj  ti

such as economics (Yang et al. 2009; Wang et al. 2012) and climatology (Elsner et al. 2009). In this study, we are particularly interested in the application of visibility graphs to the analysis of seismic sequences (Telesca and Lovallo 2012). In such an application, the nodes in the visibility graphs are considered to be seismic events distributed over time. For any seismic time series, two characteristics are attributed to each event: (a) its occurrence time, and (b) the value associated to the event, here considered as the magnitude. The obstacles in the time-magnitude space are the vertical lines (or sticks) between the time axis and the magnitude of the event. Two events are in connection, or visible to each other, whenever the magnitude stick of no other event interrupts their linear connection. In mathematical form, events i and j are visible to each other if they satisfy the inequality

ð1Þ

where y is the value associated with the event and t is the time of the event. The index p indicates any event occurring between events i and j. It follows from this that the visibility graph generated from a time series holds the following conditions. Connectivity: each event is visible to the two events to its immediate right and left sides, if there is any. Directivity: the graph is considered undirected by definition; that is, the algorithm explaining the connections between events is developed without defining a direction for the links between the events. And invariance: scaling or time-shifting the series does not change the resulting visibility graph, provided the transformation is done under affine conditions (Lacasa et al. 2008). Figure 1 shows an example of a visibility graph for a portion of one of the magnitude-time series to be Tce = 3.62 × 105 min TC = 4.40 × 105 min 3.99 × 105 2.46 × 105

4.5 8

4

4

3.5

3 3

Mw 3.48 06-09-2006 12:39:24

Mw 3.08

Mw 3.69

Mw 3.38

Mw 3.18

Mw 3.38 12-20-2005 11:12:00

2

3.0

2

Mw 2.98

3

Mw 3.99 09-05-2005 05:55:12

4 Mw 4.39 08-07-2005 23:44:54

Magnitude (Mw)

4.0

2.5 01

02

03

04

05

06

07

08

09

10

11

2005

12

01

02

03

04

05

06

2006

Time (year, month)

Figure 1 Illustration of the visibility graph method, as applied to a subset of events in a window of the magnitude-time series of the north Iranian seismic region of Kopeh Dagh. The events are identified by vertical thick lines (or sticks) on the horizontal time axis with size equal to the magnitude as indicated by the labels at the bottom of each stick. These sticks are considered to be the obstacles in the time-magnitude space where the graph resides. The visibility graph is composed by the nodes, represented with circles at the tip of each event stick, and the edges or links connecting them. The values inside each node correspond to the connectivity degree, k. This value indicates the number of straight lines that can be drawn between two (mutually visible) nodes without intersecting any obstacle (event stick). The values and time ranges shown at the top correspond to the time difference between events visible by the last event. These times are given in minutes and correspond to individual Tc values. The mean\Tce[value associated with this event is also shown at the top right. The mean value of all\Tce[for the events in the graph is the window mean interval connectivity time, \Tce[

N. Khoshnevis et al.

considered later in this study. In this application, the visibility graph analysis is one-dimensional (i.e., there is no actual physical space but simply events in a catalog), and the magnitude sticks play the role of virtual obstacles distributed along the time line as defined by the occurrence of the events, which represent the nodes of the graph. For each event i we compute the connectivity degree, k, which is the sum of the connections to all other events j visible by the ith event. Events are often categorized in small magnitude bins with DM ¼ 0:1, and we plot the connectivity degree as a function of the bin magnitude. Note that two earthquakes with the same magnitude can have different k values, which ultimately depend on the occurrence time of the event and whether the time-neighboring events are of a larger or smaller magnitude. It follows from this description that, for a given earthquake catalog, the visibility graph offers an alternative mathematical approach to study the correlations that exist between seismic events and their distribution (occurrence) in time. Such analysis can be explored in global terms (i.e., for all the events in a catalog), where long-range correlations might be more relevant in terms of the inter-visibility between large earthquakes, and the buildup of small earthquakes leading to greater magnitude events; or in local terms (i.e., for limited-time windows), where the occurrence of foreshocks and aftershocks may be more dominant. We explore such correlations in terms of how the properties of the visibility graph relate to the seismicity of the region. For instance, plotting k against magnitude results on a scattered set of points which have been shown to be acceptably represented by a linear regression. We refer to this correlation as the k–M relationship. This is of relevance because it happens that, for a significant number of earthquake catalogs from distinctive seismic zones, the slope m of the k–M relationships show a linear regression with the seismicity b value of the seismic zones under consideration. This observation seems to hold for a universal sample of m and b values from different regions, as suggested by Telesca et al. (2013, 2014b). Telesca et al. (2014b) also observed that it was reasonable to draw a relationship between the distribution of the events when considering whether these were connected (visible to

Pure Appl. Geophys.

each other) or not, and time. Let Tc be the interval connectivity time, which is nothing but the time difference between two inter-visible events (see Fig. 1); and \Tce[ the average of all Tc values computed for event e; then it is possible to compute the visibility graph mean interval connectivity time \Tc[, which is the mean value of all\Tce[in a given time sequence of events. Because \Tc[ can be obtained for any sub-graph within a larger graph, that is, for any time window in a larger magnitude-time series of events, then it is possible to investigate the evolution of \Tc[ (and that of m and b) for a moving window—with a constant number of events—sliding along the magnitudetime series. When doing so, Telesca et al. (2016) suggested to associate the values of each subsequence with the last event in the subsequence. Following this approach, Telesca et al. (2016) found that the variation of \Tc[ over time exhibited a possible relationship with the occurrence of earthquakes. This and the aforementioned relationship between m and b are aspects we explore for our region and yearsinterval of interest.

3. Seismicity of Northern Iran We are interested in applying the visibility graph method to analyze the seismicity of northern Iran. This region is part of the Iranian plateau, on the Himalayan-Alpine seismic belt. It is confined by the relative movements between the Arabian, Eurasian, eastern Asia-Minor, and Indian plates; and has a long history of large magnitude (M [ 7) earthquakes that are well documented dating back to the eighth century (e.g., Berberian 1981). According to the tectonic settings and geologic provinces of the plateau, the seismic activity in Iran has been categorized into different seismic zones. These vary between four to nine major seismic zones, in the more traditional models (e.g., Stocklin 1968; Takin 1972; Berberian and Arshadi 1976), and up to twenty to twenty-three seismotectonic provinces in the most elaborate ones (e.g., Nowroozi 1976; Tavakoli and Ghafory-Ashtiany 1999). We adopt the model proposed by Mirzaei et al. (1998) with modifications introduced by Karimiparidari et al. (2013). In it, Iran is divided

Analysis of the 2005–2016 earthquake sequence

and Yeats 1999), and a few M [ 7 earthquakes in 1042, 1721, and 1780 (Jones 1834). More recently, this region was struck by the Mw 6.1 1997 Ardabil earthquake near the city of Ardabil, and the Mw 6.4 2012 Tabriz earthquake northeast of the city of Tabriz. These earthquakes caused extensive damage and took the lives of more than 1500 people. The seismicity in the north-central region of Alborz is dominated by multiple fault systems, including the Talesh, Rubdar, North Alborz, and North Tehran faults, which also have a history of producing strong ground shaking. According to Ambraseys and Melville (1982), Tehran was devastated by severe M [ 7 earthquakes in 743, 958, 1177, 1665, and 1830. This region has also seen some significant recent seismic activity (Berberian and Yeats 1999), including the Mw 7.4 1990 Manjil– Rudbar earthquake, which caused numerous deaths and damage to the region in the south Caspian depression, south from the city of Rasht and northwest from Tehran. Last, there is the Kopeh Dagh seismic zone to the east and northeast. This region is dominated by the main Kopeh Dagh fault system, which exhibits active

into six seismic regions: Azerbaijan, Alborz, Kopeh Dagh, Zagros, Central-East Iran, and Makran. Each of these seismic regions, and in general, the whole Iranian plateau, have been the subject of various studies from different perspectives in the past (e.g., Zafarani and Farhadi 2017; Sedaghati and Pezeshk 2017; Golara 2014; Nemati 2015). Here, our focus is on the northern part of the country, for which we define a region of interest between longitudes 43.5 E and 61.5 E, and latitudes 34 N and 40 N, as shown in Fig. 2. Although this region encloses part of the Zagros and Central-East seismic zones, we limit our analysis to the seismicity of Azerbaijan, Alborz, and Kopeh Dagh. Northern Iran houses about 41% (32 million) of the total population of the country, and has suffered devastating earthquakes in the past (e.g., Mehrain 1990; Ghafory-Ashtiany 1999; Razzaghi and Ghafory-Ashtiany 2012). At the northwest, the seismic zone of Azerbaijan is strongly controlled by the North Tabriz Fault system in the vicinity of the city of Tabriz, shown in Fig. 2. Historical accounts document the occurrence of strong M [ 6 earthquakes in this region as far back as the ninth century (Berberian

32˚

40˚

48˚

56˚

64˚

72˚ 48˚

Kazakhstan Russia Black Sea

Georgia

Uzbekistan

Turkey

40˚

40˚

Turkmenistan

Syria Afghanistan Iraq

Azerbaijan Azerbaijan Ardabil Ardabil

Tabriz Tabriz

Caspian Sea

Pakistan Saudi Arabia

24˚

Ashgabat Ashgabat

38˚

Latitude

32˚

Iran

Rasht Rasht

Kopeh Kopeh Dagh Dagh Alborz Alborz

36˚

Mashhad Mashhad

Tehran Tehran

Zagros

Central-East 34˚ 44˚

46˚

48˚

50˚

52˚

54˚

56˚

58˚

60˚

Longitude Figure 2 Region of interest, highlighting the seismic zones that are the focus of this study: Azerbaijan, Alborz, and Kopeh Dagh. The top-right locationmap shows Iran and the selected area with respect to neighboring countries. The solid dark lines show the boundaries of the additional Zagros and Central-East seismic zones. This zonation follows the division proposed by Mirzaei et al. (1998) and later modified by Karimiparidari et al. (2013). The background shows fault lines and the shaded relief

N. Khoshnevis et al.

tectonic displacements along a distance of more than 500 km (Trifonov 1978). This fault is responsible for the Mw 7.3 1948 earthquake, which struck the capital city of Ashgabat in Turkmenistan and destroyed more than 30 villages in Iran. While no large magnitude events have been registered in recent years, historically, the Kopeh Dagh seismic zone is responsible for the Ms 7.1 earthquake in 10 A.D. (Berberian and Yeats 2001) near Ashgabat, and two significant earthquakes in 1209 and 1405 at the boundary between the Neyshabur and Binalud faults near the city of Mashhad (Berberian and Yeats 1999). It follows from this description that the region of interest is one of significant seismic activity, which goes back centuries. We are, however, interested in the most recent instrumental seismicity. In particular, we focus on events recorded in the last decade, between January 2005 and December 2016. We limit our analysis to this magnitude-time sequence period because in a previous study we observed that the latest decade of seismic records in northern Iran offered the most complete dataset of recorded events, especially for the case of small M\4 earthquakes (Khoshnevis et al. 2017). We compiled a catalog of all recorded earthquakes using data downloaded from the International Institute of Earthquake Engineering and Seismology, IIEES (http://www.iiees.ac.ir/en/eqcatalog/). In general, data collected from IIEES contains a mixture of earthquake magnitude scales, including: moment, Mw ; local, ML ; body wave, mb ; surface wave, Ms ; and duration MD magnitudes. However, for our time-period of interest (2005–2016), all the events had local magnitude (ML ) values, which we, in turn, and in order to maintain international consistency with other studies, converted to moment magnitude (Mw ) using the relationships defined in Zare et al. (2014). The whole dataset we downloaded from IIEES can be seen in Fig. 3a, which shows events scaled by their magnitude Mw at their epicenter location for the three seismic zones of interest. We refer to the catalog composed by the events shown in this figure as the whole catalog. We then determine the magnitude of completeness, Mc , for each region, and based on the value of Mc , we filter out the smaller events and obtain a reduced dataset of events to which we refer to as the complete catalog.

Pure Appl. Geophys.

There are different alternatives to obtain Mc . In this initial part of the analysis for the whole subcatalogs of the three regions, we determined Mc using the time-dependent MAXC method introduced and implemented in the ZMAP software by Wiemer and Wyss (2000). This method allows one to investigate the behavior of Mc as a function of time, and proposes to set Mc to a value near the observed maximum over time. Figure 4 shows the evolution of Mc for the three seismic zones, including its variability over time, and indicates the chosen value for each region. In turn, Fig. 3b shows the complete catalog that results after removing all the events with magnitudes smaller than Mc . We also considered a declustered catalog, where foreshocks and aftershocks were removed from the complete catalog following a Poissonian occurrence model using the method introduced by Gardner and Knopoff (1974). While this is consistent with previous catalog compilations available for the region (e.g., Zare et al. 2014), we recognize that this declustering approach is not calibrated for Iran. Nonetheless, we were interested in investigating the impact that declustering may have on the visibility graph analysis, an aspect we discuss later. Subsequently, we obtained the b value in the Gutenberg–Richter law for each region following the maximum likelihood estimation (Aki 1965), in which b¼

log ðeÞ  10 DM  ; M  Mmin  2bin

ð2Þ

where e is the mathematical constant or Euler’s number, M is the average magnitude, Mmin is the minimum magnitude in the sample, and DMbin is the bin size used to discretize the magnitude scale. Here, the value of Mmin is taken as the completeness magnitude Mc obtained or selected for each region’s sub-catalog. In this study we assumed DMbin ¼ 0:1. In addition to using Mc to determine b for each sub-catalog, we also use it as the threshold for limiting the minimum magnitude to be considered in the construction of the visibility graphs for each region. Table 1 shows the number of events of the whole catalog, along with the number of events for each region’s complete and declustered catalogs. It also includes the corresponding values of Mc , which happened to be the same for both catalogs, and the

Analysis of the 2005–2016 earthquake sequence

(a)

40˚

Azerbaijan Azerbaijan Kopeh Kopeh Dagh Dagh

Caspian Sea

Latitude

38˚

Alborz Alborz

36˚

34˚

(b)

40˚

Azerbaijan Azerbaijan Kopeh Kopeh Dagh Dagh

Caspian Sea

Latitude

38˚

Alborz Alborz Magnitude (Mw) 3.0 36˚

4.0 5.0 6.0

34˚ 44˚

46˚

48˚

50˚

52˚

54˚

56˚

58˚

60˚

Longitude Figure 3 Instrumental seismicity of northern Iran for the 2005–2016 period, considering only the events in the seismic zones of interest, namely Azerbaila, Alborz, and Kopeh Dagh, for the whole (a) and complete (b) catalogs. The size of symbols are proportional to the magnitude of the events, as indicated by the artificial scale shown on the right. The background shows fault lines and the shaded relief

4.0 Mc δMc

3.8 3.6 3.4

Mc

3.2 3.0 2.8 2.6 2.4 2.2 2.0

Azerbaijan ’05 ’06 ’07 ’08 ’09 ’10 ’11 ’12 ’13 ’14 ’15 ’16

Alborz

Kopeh Dagh

’05 ’06 ’07 ’08 ’09 ’10 ’11 ’12 ’13 ’14 ’15 ’16

’05 ’06 ’07 ’08 ’09 ’10 ’11 ’12 ’13 ’14 ’15 ’16

Time (Year)

Figure 4 Completeness magnitude, Mc , as a function of time (in years) for the three seismic zones in the region of interest, corresponding to the timemagnitude sequences of earthquakes during the 2005–2016 period

N. Khoshnevis et al.

Pure Appl. Geophys.

Table 1 Total number of earthquakes for the whole, complete, and declustered catalogs, and seismic parameters Mc and b (including standard deviation) for the complete catalogs of the three seismic zones of northern Iran Region

Azerbaijan Alborz Kopeh Dagh

Number of events in catalog

Mc

Whole

Complete

Declustered

1047 1898 640

202 575 117

102 440 96

separate values of b, including the standard deviation, rb , which is calculated as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 PN  i¼1 Mi  M 2 ð3Þ ; rb ¼ 2:3 b N ð N  1Þ where N is the number of events in each catalog (Shi and Bolt 1982). Figure 5 shows the final time-magnitude sequences for the complete catalogs of Azerbaijan, Alborz, and Kopeh Dagh. We omit here the sequences of the declustered catalogs, which are, in general, very similar to those shown in Fig. 5. Note from Table 1 that in the case of Azerbaijan and Kopeh Dagh, we consider events in the catalog with Mw  3:5, whereas in the case of Alborz all events are Mw  3:0.

4. Results We generate the visibility graphs for each subregion of interest using the time-magnitude sequences of both the complete and declustered catalogs. In the case of the complete catalog, that corresponds to the sequences shown in Fig. 5. In every case, the numbers of nodes in each graph are the same as those shown in Table 1; and the links between inter-visible events are established following the condition in Eq. (1), as in the example shown in Fig. 1. (We do not include a visualization of the graphs because the links are so many, that it is only practical to visualize the graphs of short sequences). We collect information about the number of inter-visibility links associated with each event (connectivity degree, k) and categorize the events in magnitude bins of size DMbin ¼ 0:1, as mentioned in the Methodology section. Figure 6 shows the scattered distribution of

3.5 3.0 3.5

b value Complete

Declustered

0.85 ± 0.05 0.82 ± 0.03 0.78 ± 0.05

0.88 ± 0.07 0.79 ± 0.03 0.73 ± 0.05

events in the magnitude-connectivity degree plane for each subregion, corresponding to the complete catalogs. The figure shows that, in general, the value of k increases with Mw ; and includes the results of obtaining linear k–M regressions for each dataset along with the values of the slope m for the three seismic zones, namely 10.05, 9.85, and 8.71 for the Azerbaijan, Alborz, and Kopeh Dagh regions. Similar results regressions were obtained for the declustered catalogs, which yielded m values of 8.92, 8.97, and 7.96, respectively. Next, we examine the relationship between the b values from Table 1 and the k–M slopes. Figure 7a, b show the scattered results of the m versus b for the complete and declustered catalogs, respectively. Along with the data-points obtained for our analysis of the three seismic regions of northern Iran, in each plot, we also include the data points corresponding to the magnitude-time sequences of the Mexican subduction zone provided by Telesca et al. (2013), the Pannoninan seismic zone from Telesca et al. (2014b), the experimental results obtained by Telesca et al. (2014a), and the result from a similar preliminary analysis done by Azizzadeh-Roodpish et al. (2017) for California. These figures also include the linear trends of the data points, as they are aggregated (i.e., each regression reflects the addition of new datapoints from different studies). In the figure, we indicate which data points were considered for each of the computed linear regressions and their corresponding correlation coefficient, R. Here, as well as in other figures, we present our results for both the complete and declustered catalogs to shed light upon the local differences for the case of northern Iran, and to contribute, globally, to the discussion of the appropriate treatment and selection of catalogs for the use of the visibility graph method.

Analysis of the 2005–2016 earthquake sequence 7.0

2012 Mw 6.4 East Azerbaijan

Azerbaijan

6.0 5.0 4.0

Magnitude (Mw)

3.0 7.0 2010 Mw 5.8 Damghan

Alborz

6.0 5.0 4.0 3.0 7.0

Kopeh Dagh

6.0

2012 Mw 5.4 Neyshabur

5.0 4.0 3.0

’05

’06

’07

’08

’09

’10

’11

’12

’13

’14

’16

’15

Year

Figure 5 Time-magnitude sequences for the northern Iran seismic zones of Azerbaijan, Alborz and Kopeh Dagh, during the period between 1 January 2005 and 31 December 2016. The event occurrences are represented by vertical bars or sticks of length equal to the Mw magnitude of each earthquake. These sequences correspond to the complete catalog. As shown in Fig. 1, the dots at the top of each stick represent the nodes of the visibility graphs. The edges or links are omitted for visual convenience. Highlighted in the figure are the largest events in each sub-catalog, namely the 2012 Mw 6.4 East Azerbaijan, 2010 Mw 5.8 Damghan, and 2012 Mw 5.4 Neyshabur earthquakes

Connectivity Degree (k)

50

60

(a) Azerbaijan

40 30

30

(b) Alborz

50

25

40

20

30

15

(c) Kopeh Dagh m = 8.71

m = 10.05

20

m = 9.85

20

10

10 5

10

0

0

3.0

3.5

4.0

4.5

5.0

5.5

Magnitude (Mw)

6.0

6.5

0

2.5

3.0

3.5

4.0

4.5

5.0

Magnitude (Mw)

5.5

6.0

3.2

3.6

4.0

4.4

4.8

5.2

5.6

Magnitude (Mw)

Figure 6 Scattered distribution of events in the magnitude-connectivity degree plane and linear regressions obtained for the k–M relationships for the three seismic regions in northern Iran. The value next to each of regression line corresponds to the slope of the line, which is referred here as the k–M slope

We note, for instance, that the data points in Fig. 7 from Telesca and Lovallo (2012) and Telesca et al. (2013) used whole catalogs, whereas those from Telesca et al. (2014b) and (2016) used declustered catalogs. Note that, overall, the linear trends shown in Fig. 7 hold independently of the differences in the studies in terms of number of events, seismic characteristics, and the treatment of the catalogs. Also

note that these trends tend to improve with increasing number of data points, as reflected by the increasing value of R. Regarding this, we should note that the data points in Fig. 7 cover what would otherwise be considered a large span of b values. While we recognize that most earthquake data lead to b values ranging between 0.8 and 1.2, it is also known that there are regions with low- and high-b values (e.g., Singh et al. 1983; Pacheco et al. 1992; Nakaya

N. Khoshnevis et al.

(a)

Pure Appl. Geophys.

(b) 2.5

2.0

R = 0.985 R = 0.952 R = 0.948

R = 0.986 R = 0.959 R = 0.948

R = 0.904

R = 0.768

b

1.5

1.0 Northern Iran (this study) So. California (Azizzadeh et al., 2017) Mexican Subduction (Telesca et al., 2013) Pannonian Region (Telesca et al., 2014b)

0.5

Synthetic Data - R1 (Telesca et al., 2014a) Synthetic Data - R5 (Telesca et al., 2014a)

0 0

5

10

15

20

25

0

5

10

15

20

25

m

m

Figure 7 Correlation between k–M slope (m) and the b value as drawn from the results of the present study for the region of northern Iranian using the complete (a) and declustered (b) catalogs, along with other previous studies including analysis of the Mexican subduction zone (Telesca et al. 2013), the Pannonian seismic zone (Telesca et al. 2014b), southern California (Azizzadeh-Roodpish et al. 2017), and results from two experiments (Telesca et al. 2014a). The lines represent linear regressions obtained to fit the different data points, considering different combinations. The values of the correlation coefficient, R, are indicated for each regression line

2006). In addition, Scholz (1968) reported that experimental microfractures in rocks, similar to those observed in earthquakes, exhibit b values ranging between 0.11 and 2.58. Considering that, by definition, b represents the proportion between the amount of large and small events in a region, the trends shown in Fig. 7 contribute to the suggestion made by Telesca et al. (2014b) about the universal character of the relationship between b and m. According to these results, a universal relationship between b and m can be expressed as b ¼ 0:073 þ 0:084m

ð4Þ

for the case in which we consider the complete catalogs for northern Iran, and b ¼ 0:077 þ 0:085m

ð5Þ

for the case in which we consider the declustered catalogs for northern Iran. This indicates that once combined with the data points of other regions, the option of considering either catalog leads to very similar results as, on average, the data points oscillate about the same universal relationship. Locally, however, for the case of northern Iran trends alone, it

seems that the declustered catalog leads to a closer fit with the universal relationship. Separately (not shown here for brevity), we also observed that, for the particular sequences under study, the results seem to be more sensitive to less conservative values of Mc , as the choice of a smaller Mc leads to larger datasets. Telesca and Lovallo (2012) and Telesca et al. (2013), however, pointed out that for sufficiently long sequences, the threshold magnitude have only small effects on the graph parameters. To further explore the sensitivity of the graph properties to the number of events in each catalog and the value of the minimum magnitude, we randomly picked a significant number of subsequences from within the initial catalog compiled for each region, and repeated the analysis for each subsequence. In total, for each region, we extracted 200 new subsequences from the whole catalog. The number of events in each subsequence, n, was varied randomly but chosen to be large enough to represent the seismic characteristics of each region. In particular, the minimum size of each subsequence was set to be n  200, and the maximum size in the sequence was set to be as large as the original whole catalog (see

Analysis of the 2005–2016 earthquake sequence

method to determine Mc systematically, instead of the MAXC method. The GFT method follows the recommendations of Wiemer and Wyss (2000) and Wiemer (2001), and leads to less conservative values of Mc , which is preferable to keep the subsequences large enough. Figure 8 shows the results for this analysis on random subsequences for both the complete and declustered catalogs. In particular, Fig. 8a, d show scattered points for all the individual m and b pairs for all the subsequences that were randomly picked from the complete and declustered catalogs of each region in northern Iran (small symbols), respectively. Figure 8b, e show the variability of m, and Fig. 8c, f

Table 1). We forced the random subsequences to progress positively in time without altering the natural occurrence of events. In other words, we randomly determined the initial event and the subsequence size (number of events to be considered), and then picked that number of events following the initial earthquake in the subsequence. Next, we determined the value of Mc and b to obtain the complete and declustered catalogs for each region, created the graph for all events with M  Mc , and extracted the connectivity degree of the events in each subsequence. Because this needed to be done for a large number of subsequences, in this part of the analysis we used the goodness-of-fit test (GFT) 1.2

(a) Complete

10

1.0

(b)

(c)

1.1 0.9

0.8 0.7 0.6 0.5 0.4 1.2

(d) Declustered

1.1

b

0.9

8.05 ± 0.76 (9.38%) 7.27 ± 0.45 (6.23%)

6

0.8 0.7 0.6 Azerbaijan (Az) Alborz (Al) Kopeh Dagh (Ko)

0.5

6

0.6

0.5

4

0.4

(e)

1.0 8.49 ± 0.61 (7.16%)

9

8 8.06 ± 0.36 (4.41%)

7

6

6.39 ± 0.42 (6.59%)

0.71 ± 0.04 (5.41%) 0.65 ± 0.04 (5.43%)

(f)

7

8

9

10

11

0.78 ± 0.03 (3.31%)

0.8

0.7

0.6

0.69 ± 0.04 (6.14%) 0.60 ± 0.03 (4.76%)

0.5

4 5

0.7

5

5

0.4 4

0.8

0.9

Mean ± STD (CV) of m

1.0

7

10

Mean of random catalogs Regional catalogs Universal data Random catalogs Regional catalogs

9.50 ± 0.70 (7.37%)

8

Mean ± STD (CV) of b

b

0.9

0.81 ± 0.04 (4.80%)

9

Mean ± STD (CV) of b

Mean ± STD (CV) of m

1.0

0.4 Az

Al

Ko

Az

Al

Ko

m

Figure 8 Correlation between the k–M slope (m) and b values for the random subsequences extracted from the catalogs of the three northern Iranian seismic regions considered in this study (scattered symbols, 200 each region, 600 total) for both the complete (a) and declustered (d) catalogs, including the mean values (empty symbols with thick border) and the data points obtained from the regional analyses (solid symbols), along with the linear regressions of each sample as indicated in the legend; and mean 1 standard deviation bar plots for m (b, e) and b (c, f) for the complete and declustered catalogs, respectively. In the latter, the values in the parentheses corresponds to the coefficients of variation

N. Khoshnevis et al.

show the variability of b. These variability plots show the mean values of each parameter for all the random subsequences and the amplitudes of 1 standard deviation, as well as the coefficients of variation (in percentage). Figure 8a, d also include the data points corresponding to the regional analysis of the catalogs, the points corresponding to the mean values of m and b, and the universal linear regression from Eqs. (4) and (5), as well as the linear regressions for northern Iran when using the full-sequence catalogs (dashed thin line) and the random subsequence catalogs (continuous thin line). In essence, the results in Fig. 8 show that, despite the variability in the random subsequences and the differences in the computation of Mc , both of which are sources of additional uncertainty, the analysis produces an outcome in fair agreement with the universal relationships. Nonetheless, the analysis of the random subsequences in the complete catalogs produces some outliers, whereas the processing of the declustered subsequences tends to concentrate more evenly around the universal trend. Note that, likely due to the removal of foreshocks and aftershocks, the declustered random subsequences exhibit a clearer distinction between the three seismic zones. While these sequences do not represent the dominant seismicity pattern of the region, the comparison between the mean values of the random sequences analysis versus the complete catalogs (large, empty and solid symbols in Fig. 8) indicates that there exists only a small bias, which is well within the standard deviation of the different value samples. Note also that the values of m in Fig. 8 are slightly smaller when obtained with the random sequence analysis than when done for the regional sequences, whereas the b values seem more stable. This explains why the regressions of the random subsequences are less steep than the universal regressions, in both cases. In a separate initial analysis, not included here for brevity, we observed that lower, less conservative selections of Mc can lead to random subsequences that yield regressions in better agreement with the universal trend. Using lower values of Mc , however, may be misleading. In this sense, larger catalogs in better instrumented regions may shed light on the stability of the m–b relationship, which is the reason we included the data point

Pure Appl. Geophys. Figure 9 c Variation of \Tc[, k–M slope (m) and the b value with respect to time for the seismic regions of Azerbaijan (a) and Alborz (b) considering a moving window analysis of the visibility graphs of subsequences of n ¼ 50 consecutive events, along with the magnitude and the connectivity degree (k) of events in the 2005–2016 period. Black triangles indicate the occurrence of significant earthquakes in each of the regions, with the corresponding magnitude at the top of each symbol

from Azizzadeh-Roodpish et al. (2017) for southern California in Fig. 7, which falls well in line with the universal regression. We now investigate the relationship of the parameters obtained through the visibility graph analysis over time. Here, the visibility graph analysis is done by windows of equal number of events moving in time along the catalog sequence. The number of events in the window is kept fixed independently of the time between them, and the results are associated with the last event in the window. In this case, we are interested in using a small number of events to capture the relevance of each new event as the window moves in time. We tried different window sizes between 20 and 100 events per window for each seismic zone. As it is natural to expect, there are trade-offs between different window-sizes. Smaller window sizes emphasize the local changes in the time series, but are less reliable when it comes to the seismicity parameters (e.g., b value). Larger window sizes, on the other hand, offer less insight into the time dependence of the window sequence characteristics. We chose a window size n ¼ 50, which is consistent with the suggested minimum number of events acceptable to estimate b according to Woessner and Wiemer (2005). In this windowing analysis, in addition to b and m, we are interested in computing the value of \Tc[ explained in the Methodology section. Figure 9 shows the earthquake sequence for the regions of Azerbaijan and Alborz, and the variation of k, m, b, and \Tc[ with time. Although we performed the analysis for all three seismic zones, we concentrate in this two regions because the windowed sequence of Kopeh Dagh was too short and did not add value to the discussion. In the two cases in Fig. 9 we can see here how larger values of k are associated with the

Analysis of the 2005–2016 earthquake sequence

(a) Azerbaijan (min. × 104)

2012 Mw 6.4 Tabriz

10 8 6 4 2 0 10

1.2 b m

1.0

b

m

8 6

0.8

4

0.6

60

k

40 20 0 7.0

Mw

6.0 5.0 4.0 3.0

(b) Alborz (min. × 104)

5.0

5.4

2010 Mw 5.8 Damghan

5.2

3 2 1 0 10 8

1.3

b m

1.0

6

b

m

5.1

4

4 0.8

2 60

k

40 20 0 7.0

Mw

6.0 5.0 4.0 3.0

’05

’06

’07

’08

’09

’10

’11

’12

’13

’14

’15

’16

Year

larger magnitude earthquakes, and how the resulting values of m for each 50-event sequence compares with b. For the case of Alborz, we observe a good

correlation between m and b for the period prior to the 2010 Mw 5.8 Damaghan earthquake, but after this event the correlation is not as clear. For Azerbaijan,

N. Khoshnevis et al. 1.2

Pure Appl. Geophys.

(a) Azerbaijan

1.1

(b) Alborz

b = 0.60 + 0.037 m R = 0.70

1.0 b = 0.40 + 0.061 m R = 0.74

b

0.9 0.8 0.7 0.6 0.5 0.4 3

4

5

6

7

8

9

10

m

11

3

4

5

6

7

8

9

10

11

m

Figure 10 Correlation between m and b based on the fixed-size windowing analysis done for the regions of Azerbaijan (a) and Alborz (b), along with the linear regressions of each case for the complete catalog. While a direct comparison is not appropriate, the universal regression is shown here for reference

the relationship between m and b is less evident, except for the timespan around the occurrence of the 2012 Mw 6.4 Tabriz earthquake, in which case both all three parameters (b, m and \Tc[) show a clear drop in their values. By contrast, in the case of Alborz, we do not observe any clear correlation between the behavior of \Tc[ and the occurrence of the largest (Mw  5) earthquakes in this region. We notice that while the occurrence of the larger events does seem to coincide with drops in m and b, there are other similar changes in the values of m and b that are not associated with any significant event, making it difficult to draw strong conclusions solely based on the results for the region of Azerbaijan and the 2012 Tabriz earthquake. It is possible, however, that regions with larger events or catalogs with lower Mc lead to better insight on this matter. The better correlation observed between m and b for Alborz than for Azerbaijan is further explored in Fig. 10, which presents the condensed (scattered data-point) results for all the values k and b from Fig. 9. The plots in this figure confirm that the values of m and b correlate better for Alborz. Also included in Fig. 10, for reference, is the universal regression. While a direct comparison with the local regressions is not appropriate, the closer proximity of the Alborz

trend line to the universal regression line suggests that the larger number of data points improves the quality of the regression, bringing it closer to the universal correlation line. The results for Kopeh Daght, not included here for brevity, show the poorest correlation, likely due to the least amount of data points, which reinforces the point made about the quality of the regression being dependent on a sufficiently large number of data points.

5. Conclusions We studied the seismicity of the three main seismic regions in northern Iran in the time period between January 2005 and December 2016 using an approach based on the visibility graph method. We tested the applicability of this method for the specific region of interest and in reference to previous results from similar studies. The results confirm previous observations about the correlation that exists between the connectivity parameter k, the magnitude of the events in the sequence, the slope m of the linear regression adjusted to the relationship between k and M, and the correlation of the b value from the Gutenberg–Richter law with m. We used the data

Analysis of the 2005–2016 earthquake sequence

points computed from the visibility graph analysis of the seismicity of northern Iran to obtain updated mathematical expressions for a universal relationship between b and m, including data collected from previous studies done for other regions as well as data from experiments. We found the relationships to have a good level of similarity, independently of the alternative manipulations of the seismic catalog, but noted that the results were somewhat sensitive to the completeness magnitude. Overall, our results are indicative of the general nature of the relationship between m and b. We also explored the relationships between the visibility graph properties and the seismicity parameters for a random selection of earthquake subsequences and for a progression of subsequences in time. In these cases, we found the relationship between b and m to hold for subsequences with a sufficiently large number of events (n  200), but noticed they are in better agreement with the universal regressions when larger, entire regional sequences are considered. For the particular regions considered in this study, we were not able to detect conclusive evidence of other relationships between the visibility graph properties and the temporal occurrence of large earthquakes, as it has been suggested in previous similar studies. We believe this is due to the lack of larger magnitude events (Mw [ 6) or relatively high values of Mc in the selected catalogs. While the connection between the topological properties of the visibility graph of an earthquake sequence and the seismicity parameters of a region is an attractive finding, it remains to be seen how this can be put to use in a more practical scientific sense. We believe future efforts in understanding the mathematical nature of the observed relationships is a natural first step to follow. It will also be necessary to test the relationships in additional seismic regions with uncharacteristic b values, especially in the ranges outside the typical b values (0.8–1.2). It would also be of interest to build visibility graphs based on other earthquake parameters (e.g., moment) or considering directionality, and investigate other relationships. At this point, nonetheless, the method as used here seems to provide an alternative and interesting approach to the analysis of the seismicity of a region.

Acknowledgements We are thankful for the comments and suggestions of two anonymous reviewers, which helped us improve the original manuscript. The data used to compile the seismic catalog employed in this analysis was obtained from the International Institute of Earthquake Engineering and Seismology, IIEES (http:// www.iiees.ac.ir/en/eqcatalog/). This research was possible thanks to support from the Center for Earthquake Research and Information (CERI) at the University of Memphis. CERI is designated as a Center of Excellence by the Tennessee Board of Regents and is funded in part by the State of Tennessee under State Sunset Laws (SB 1510 and HB 1608, 2015–2016). REFERENCES Ahmadlou, M., Adeli, H., & Adeli, A. (2010). New diagnostic EEG markers of the Alzheimers disease using visibility graph. Journal of Neural Transmission, 117(9), 1099–1109. doi:10.1007/ s00702-010-0450-3. Aki, K. (1965). Maximum likelihood estimate of b in the formula log N ¼ a  bM and its confidence limits. Bulletin of the Earthquake Research Institute Tokyo University, 43, 237–239. Ambraseys, N. N., & Melville, C. P. (1982). A history of Persian earthquakes. Cambridge: Cambridge University Press. Azizzadeh-Roodpish, S., Khoshnevis, N., & Cramer, C. H. (2017). Visibility graph analysis of southern California. In Proceedings annual meeting of the seismological society of America, Denver, Colorado. doi:10.13140/RG.2.2.24083.48160. Berberian, M. (1981). Active faulting and tectonics of Iran. In H. K. Gupta & F. M. Delany (Eds.), Zagros hindu kush himalaya geodynamic evolution. Washington, DC: American Geophysical Union. doi:10.1029/GD003p0033. Berberian, M., & Arshadi, S. (1976). On the evidence of the youngest activity of the North Tabriz Fault and the seismicity of Tabriz city. Geological Survay Iran Report, 39, 397–418. Berberian, M., & Yeats, R. S. (1999). Patterns of historical earthquake rupture in the Iranian Plateau. Bulletin of the Seismological Society of America, 89(1), 120–139. Berberian, M., & Yeats, R. S. (2001). Contribution of archaeological data to studies of earthquake history in the Iranian Plateau. Journal of Structural Geology, 23(2), 563–584. doi:10. 1016/S0191-8141(00)00115-2. Elsner, J. B., Jagger, T. H., & Fogarty, E. A. (2009). Visibility network of United States hurricanes. Geophysical Research Letters, 36(16), L16,702. doi:10.1029/2009GL039129. Gardner, J., & Knopoff, L. (1974). Is the sequence of earthquakes in southern California, with aftershocks removed, Poissonian? Bulletin of the Seismological Society of America, 64(5), 1363–1367. Ghafory-Ashtiany, M. (1999). Rescue operation and reconstruction of recent earthquakes in Iran. Disaster Prevention and

N. Khoshnevis et al. Management: An International Journal, 8(1), 5–20. doi:10.1108/ 09653569910258174. Golara, A. (2014). Probabilistic seismic hazard analysis of interconnected infrastructure: A case of Iranian high-pressure gas supply system. Natural Hazards, 73(2), 567–577. Gutenberg, B., & Richter, C. F. (1944). Frequency of earthquakes in California. Bulletin of the Seismological Society of America, 34(4), 185–188. Huang, H. P., & Chung, S. Y. (2004). Dynamic visibility graph for path planning. Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 3, 2813–2818. doi:10.1109/IROS.2004.1389835. Jones, H. (1834). An account of the transactions of his majesty’s mission to the court of Persia: in the years 1807-11, vol 1. James Bohn. Karimiparidari, S., Zare´, M., Memarian, H., & Kijko, A. (2013). Iranian earthquakes, a uniform catalog with moment magnitudes. Journal of Seismology, 17(3), 897–911. doi:10.1007/s10950-0139360-9. Khoshnevis, N., Taborda, R., Azizzadeh-Roodpish, S., & Cramer, C. H. (2017). Seismic hazard estimation of northern iran using smoothed seismicity. Journal of Seismology, 21(4):941–964. doi:10.1007/s10950-017-9645-5. Lacasa, L., Luque, B., Ballesteros, F., Luque, J., & Nuno, J. C. (2008). From time series to complex networks: The visibility graph. Proceedings of the National Academy of Sciences, 105(13), 4972–4975. doi:10.1073/pnas.0709247105. Long, Y. (2013). Visibility graph network analysis of gold price time series. Physica A: Statistical Mechanics and its Applications, 392(16), 3374–3384. doi:10.1016/j.physa.2013.03. 063. Lozano-Pe´rez, T., & Wesley, M. A. (1979). An algorithm for planning collision-free paths among polyhedral obstacles. Communications of the ACM, 22(10), 560–570. doi:10.1145/359156. 359164. Mehrain, M. (1990). Reconnaissance report on the northern Iran earthquake of June 21, 1990. Technical Report NCEER-90-0017, National Center for Earthquake Engineering Research (NCEER). http://hdl.handle.net/10477/651. Mirzaei, N., Mengtan, G., & Yuntai, C. (1998). Seismic source regionalization for seismic zoning of Iran: Major seismotectonic provinces. Journal of Earthquake Prediction Research, 7(4), 465–495. Nakaya, S. (2006). Spatiotemporal variation in b value within the subducting slab prior to the 2003 Tokachi-oki earthquake (M8:0), Japan. Journal of Geophysical Research, Solid Earth, 111(B03311), 13. doi:10.1029/2005JB003658. Nemati, M. (2015). Intermediate-term variations in 200 years seismicity of south of Iran. Geomatics, Natural Hazards and Risk, 7(3), 1065–1080. doi:10.1080/19475705.2015.1030785. Nowroozi, A. A. (1976). Seismotectonic provinces of Iran. Bulletin of the Seismological Society of America, 66(4), 1249–1276. Oommen, B., Iyengar, S., Rao, N., & Kashyap, R. (1987). Robot navigation in unknown terrains using learned visibility graphs. Part I: The disjoint convex obstacle case. IEEE Journal on Robotics and Automation, 3(6), 672–681. doi:10.1109/JRA.1987. 1087133. Pacheco, J. F., Scholz, C. H., & Sykes, L. R. (1992). Changes in frequency? Size relationship from small to large earthquakes. Nature, 355(6355), 71–73. doi:10.1038/355071a0.

Pure Appl. Geophys. Phillips, J. D., Schwanghart, W., & Heckmannc, T. (2015). Graph theory in the geosciences. Earth-Science Reviews, 143, 147–160. doi:10.1016/j.earscirev.2015.02.002. Raman, S. (2010). Designing a liveable compact city: Physical forms of city and social life in urban neighbourhoods. Built Environment, 36(1), 63–80. doi:10.2148/benv.36.1.63. Razzaghi, M. S., & Ghafory-Ashtiany, M. (2012). A preliminary reconnaissance report on August 11th 2012, Varzaghan-Ahar twin earthquakes in NW of Iran. Technical report, International Association of Seismology and Physics of the Earth’s Interior, and Iranian Earthquake Engineering Association. Scholz, C. H. (1968). The frequency–magnitude relation of microfracturing in rock and its relation to earthquakes. Bulletin of the Seismological Society of America, 58(1), 399–415. Sedaghati, F., & Pezeshk, S. (2017). Partially nonergodic empirical groundmotion models for predicting horizontal and vertical pgv, pga, and 5 acceleration response spectra using data from the iranian plateau. Bulletin of the Seismological Society of America, 107(2), 934–948. doi:10.1785/0120160205. Shi, Y., & Bolt, B. A. (1982). The standard error of the magnitudefrequency b value. Bulletin of the Seismological Society of America, 72(5), 1677–1687. Singh, S. K., Rodriguez, M., & Esteva, L. (1983). Statistics of small earthquakes and frequency of occurrence of large earthquakes along the Mexican subduction zone. Bulletin of the Seismological Society of America, 73(6A), 1779–1796. Skiena, S. S. (2008). The algorithm design manual (2nd ed.). London: Springer. doi:10.1007/978-1-84800-070-4. Stocklin, J. (1968). Structural history and tectonics of Iran: A review. AAPG Bulletin, 52(7), 1229–1258. Takin, M. (1972). Iranian geology and continental drift in the Middle East. Nature, 235(5334), 147–150. doi:10.1038/ 235147a0. Tavakoli, B., & Ghafory-Ashtiany, M. (1999). Seismic hazard assessment of Iran. Annals of Geophysics, 42(6), 1013–1021. doi:10.4401/ag-3781. Telesca, L., & Lovallo, M. (2012). Analysis of seismic sequences by using the method of visibility graph. Europhysics Letters, EPL, 97(5), 50002. doi:10.1209/0295-5075/97/50002. Telesca, L., Lovallo, M., & Pierini, J. O. (2012). Visibility graph approach to the analysis of ocean tidal records. Chaos, Solitons and Fractals, 45(9), 1086–1091. doi:10.1016/j.chaos.2012.06.003. Telesca, L., Lovallo, M., Ramirez-Rojas, A., & Flores-Marquez, L. (2013). Investigating the time dynamics of seismicity by using the visibility graph approach: Application to seismicity of Mexican subduction zone. Physica A: Statistical Mechanics and its Applications, 392(24), 6571–6577. doi:10.1016/j.physa.2013. 08.078. Telesca, L., Lovallo, M., Ramirez-Rojas, A., & Flores-Marquez, L. (2014a). Relationship between the frequency magnitude distribution and the visibility graph in the synthetic seismicity generated by a simple stick-slip system with asperities. PLoS One, 9(8), e106233. doi:10.1371/journal.pone.0106233. Telesca, L., Lovallo, M., & Toth, L. (2014). Visibility graph analysis of 2002–2011 Pannonian seismicity. Physica A: Statistical Mechanics and its Applications, 416, 219–224. doi:10.1016/ j.physa.2014.08.048. Telesca, L., Lovallo, M., Aggarwal, S., Khan, P., & Rastogi, B. (2016). Visibility graph analysis of the 2003–2012 earthquake sequence in the Kachchh region of western India. Pure and

Analysis of the 2005–2016 earthquake sequence Applied Geophysics, 173(1), 125–132. doi:10.1007/s00024-0151034-9. Trifonov, V. G. (1978). Late quaternary tectonic movements of western and central Asia. Geological Society of America Bulletin, 89(7), 1059–1072. doi:10.1130/0016-7606(1978) 89\1059: LQTMOW[2.0.CO;2. Varoudis, T., & Psarra, S. (2014). Beyond two dimensions: Architecture through three dimensional visibility graph analysis. The Journal of Space Syntax, 5(1), 91–108. Wang, N., Li, D., & Wang, Q. (2012). Visibility graph analysis on quarterly macroeconomic series of China based on complex network theory. Physica A: Statistical Mechanics and its Applications, 391(24), 6543–6555. doi:10.1016/j.physa.2012.07.054. Wiemer, S. (2001). A software package to analyze seismicity: ZMAP. Seismological Research Letters, 72(3), 373–382. doi:10. 1785/gssrl.72.3.373. Wiemer, S., & Wyss, M. (2000). Minimum magnitude of completeness in earthquake catalogs: Examples from Alaska, the

western United States, and Japan. Bulletin of the Seismological Society of America, 90(4), 859–869. doi:10.1785/0119990114. Woessner, J., & Wiemer, S. (2005). Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty. Bulletin of the Seismological Society of America, 95(2), 684–698. doi:10.1785/0120040007. Yang, Y., Wang, J., Yang, H., & Mang, J. (2009). Visibility graph approach to exchange rate series. Physica A: Statistical Mechanics and its Applications, 388(20), 4431–4437. doi:10. 1016/j.physa.2009.07.016. Zafarani, H., & Farhadi, A. (2017). Testing ground-motion prediction equations against small-to-moderate magnitude data in Iran. Bulletin of the Seismological Society of America, 107(2), 912–933. doi:10.1785/0120160046. Zare, M., Amini, H., Yazdi, P., Sesetyan, K., Demircioglu, M. B., Kalafat, D., et al. (2014). Recent developments of the Middle East catalog. Journal of Seismology, 18(4), 749–772. doi:10. 1007/s10950-014-9444-1.

(Received January 1, 2017, revised June 26, 2017, accepted July 6, 2017)