Insight into earthquake sequencing : analysis and interpretation of time-series constructed from the directed graph of the Markov chain model

Introduction Conclusions References


Introduction
Earthquake sequencing has been the subject of detailed research (Nava et al., 2005;Ünal andÇelebioğlu, 2011, 2014;Telesca et al., 2001Telesca et al., , 2008Telesca et al., , 2009Telesca et al., , 2011;;Cavers andVasudevan, 2013, 2015;Vasudevan andCavers, 2012, 2013) both in the regional and global sense in recent years.Nava et al. (2005)  the characteristics of the state-to-state transitions where each state is described by the earthquake occupancy of the zones.Thus, states can fall into zones of no occupancy to full occupancy at the extreme and into zones where some are occupied and some are not, and can thus be presented in binary format with "0" representing non-occupancy and "1" occupancy.The approach of Nava et al. (2005) was immediately extended to other regions (Herrera et al., 2006;Ünal andÇelebioğlu, 2011, 2014).Cavers and Vasudevan (2013) adapted the method of Nava et al. (2005) to a global catalogue which was partitioned into zones on the basis of the tectonic boundaries (Kagan et al., 2010).
The existing Markov chain model was refined by incorporating the record-breaking recurring events for each event in the catalogue under certain constraints.A directed graph representation of the modified Markov chain model was then subjected to detailed analysis for forecasting purposes (Cavers and Vasudevan, 2015).
One consequence of the approach taken by Cavers and Vasudevan (2015) and Vasudevan and Cavers (2013) is that it results in a time-series of state-to-state transition frequencies of the modified Markov chain model, x sstf (t).This time-series is for an optimized time-interval, ∆t.The fluctuations in state-to-state transitions are ∆t sampled.The time-series is a comprehensive representation of earthquake sequencing in which interaction of seismic events within and among zones are considered.Therefore, it can be subjected to a detailed analysis.
Earthquake sequencing may be considered a non-linear and non-stationary process (Kanamori, 2003;Telesca et al., 2001Telesca et al., , 2008Telesca et al., , 2009Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2012).In earthquake sequencing, earthquakes are viewed as part of a point process, with earthquake events occurring at some random locations in time.This means that the earthquake sequencing is dictated by the set of event times, and can also be expressed by the set of time-intervals between events.The time-series of earthquakes for any time-interval can be analyzed in many ways (Telesca et al., 2001(Telesca et al., , 2008(Telesca et al., , 2009(Telesca et al., , 2011)).
We postulate here that the non-linear and non-stationary behavior in the time-series should also be present in the time-series of the state-to-state transition frequencies Introduction

Conclusions References
Tables Figures

Back Close
Full derived from earthquake sequencing.Hence, we consider the approaches of Telesca et al. (2001Telesca et al. ( , 2008Telesca et al. ( , 2009Telesca et al. ( , 2011) ) to be appropriate for a study here.
Non-linear and non-stationary time-series have been examined in recent years with a method known as empirical mode decomposition (EMD) and the intrinsic mode functions derived from this are useful in this regard (Huang et al., 1998).The present timeseries of state-to-state transition frequencies is suited for such a study.
In general, the time-series has non-zero amplitudes for the state-to-state transition frequencies (Cavers and Vasudevan, 2015).In this particular case, there are instances where there are no earthquakes exceeding the magnitude of 5.6 in all zones for one ∆t or for successive ∆t s.This introduces "intermittency" in the time series.
However, because of the presence of intermittency in it, an ensemble approach to empirical mode decomposition, EEMD (Wu and Huang, 2004Huang, , 2009;;Flandrin et al., 2004Flandrin et al., , 2005) ) is applied here.The intermittency problem is handled with the addition of random noise to the time-series before carrying out the EEMD.We examine the criteria used for the selection of the added noise and the ensemble number for the EEMD.
Another aspect of the study here is to ask a question if the time-series resulting from a directed graph representation of the Markov chain model of earthquake sequences exhibits power-law statistics similar to a description of fractal stochastic point processes (Telesca et al., 2001(Telesca et al., , 2009(Telesca et al., , 2011) ) to model the time-occurrence-sequence of seismic events.Quantifying the earthquake sequencing in terms of its fractal properties was done by means of the Fano factor and the Allan factor (Allan, 1966;Barnes and Allan, 1966;Lowen andTeich, 1993a, 1995;Thurner et al., 1997;Telesca et al., 2001Telesca et al., , 2009Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2012;Serinaldi and Kilsby, 2013).Since the fractal properties of the time-series studied here has never been investigated, we calculate the Fano factor and the Allan factor for the purpose of quantitative analysis.
The remainder of the paper is divided into three sections.In the next section, we show how the time-series of the state-to-state transition frequencies for a modified Markov chain model as described in Cavers and Vasudevan (2015) is generated.In the following section, we describe the EEMD procedure used and the analysis of the results Introduction

Conclusions References
Tables Figures

Back Close
Full that accrue from this procedure.We extend the approaches of Telesca et al. (2001Telesca et al. ( , 2008Telesca et al. ( , 2009Telesca et al. ( , 2011) ) to calculate the Fano factor and the Allan factor with a view to study the fractal properties of the time-series.In the last section, we discuss the results of the analysis methods and draw certain inferences about the state-to-state transition frequencies.

Directed graph representation of earthquake sequencing
A Markov chain is a discrete-time stochastic process consisting of a collection of random variables {X 1 , X 2 , X 3 , . ..} indexed by time, where for each n, the state of X n+1 is independent of the past states X 1 , X 2 , . .., X n (Çınlar, 1975).Each X i takes values in a finite set, S, called the states of the system.To build a Markov chain model we first partition the region, either local or global, into zones.Typically these zones are made up of rectangles that divide the region (Nava et al., 2005;Ünal and Çelebioğlu, 2011).
Recently, other partitions have been used.In particular, Cavers and Vasudevan (2015) used a simplified 5-zone plate boundary template as given by Kagan et al. (2010) to study global seismicity, while Ünal et al. (2014) used a seismic zones map that uses geographic information system analysis to divide Turkey into regions.For this particular study, we used the five-zone model described in Cavers and Vasudevan (2015) and give an overview of its construction here.Kagan et al. (2010) partitioned the shallow (≤ 70 km depth) events with moment magnitude, M w > 5.6 from the Global CMT catalogue (1 January 1982-31 March 2008) into 5 zone sub-catalogues using their grid-assignment schemes (Table 1).The selected catalogue consists of 6752 earthquakes with 4407 from Zone 4 (Trenches), 723 from Zone 3 (Fast-spreading ridges), 487 from Zone 2 (Slow-spreading ridges), 898 from Zone 1 (Active continent), and 237 from Zone 0 (Plate interior) respectively.For these five zones, we express a state, corresponding to a time interval ∆t, as a concatenation occurrence in zone L during the specified time interval ∆t.We use Θ = [θ i j ] to denote the transition frequency matrix, where θ i j is the number of occurrences from state i to state j .Letting s(n) represent the state for interval number n, the probability transition matrix, P = [p i j ], consists of transition probabilities, p i j , given as For a Markov chain structure given earlier for the five zones, the computation of transition frequencies and hence, transition probabilities, depend on the chosen timeinterval ∆t.Following the selection rules given elsewhere (Nava et al., 2005;Ünal and Çelebioğlu, 2011;Cavers and Vasudevan, 2015), we used a ∆t value of 9 days for the construction of the Markov chain of transition probabilities.
A finite-state Markov chain can be depicted using a digraph representation, G, where the set of possible states (binary strings of length 5) are the nodes, and an arc (i , j ) connects two states i and j if and only if p i j > 0 (Jarvis and Shier, 1996).Figure 1 shows an example of a digraph representing a Markov chain with a three zone partition, hence, there there are 2 3 = 8 states {000, 001, 010, 011, 100, 101, 110, 111} that we write in decimal format {0, 1, 2, 3, 4, 5, 6, 7}, respectively.In this figure, we do not show all of the possible transitions between states and typically an arc (i , j ) is omitted when p i j = 0. We follow the same decimal state labelling format as in Fig. 1 for our 2 5 = 32 states, that is, state "0" (representing 00000 in binary) corresponds to no earthquake occurrence in all five zones in the chosen time interval, ∆t, and state "31" (representing 11111 in binary) points to earthquake occurrences in all five zones.Table 2 shows details for defining all other states, "1" to "29".
The combinatorial structure of a digraph representation of the Markov chain model contains important information for earthquake sequencing (Cavers and Vasudevan, 2015).It is often useful to use a weight w i j for each arc (i , j ) of the digraph to get a weighted digraph.The weights have the form w i j = θ i j , w i j = p i j , or can be empirically derived from the Markov chain.To introduce spatial-temporal complexity into the 404 Introduction

Conclusions References
Tables Figures

Back Close
Full model so that transitions with earthquake occurrences at large distances have less of an impact on our model than transitions with earthquake occurrences at short distances, we follow the approach by Cavers and Vasudevan (2015) to modify the weights w i j in the weighted digraph by considering recurrences.Each earthquake (event) in a zone may have several recurring events in the record-breaking sense.The recurring events for one event in a given zone may fall into other zones or may be in the same zone.This flexibility adds to the possibility of interactions among zones.We first form the network of recurrences as described by Davidsen et al. (2008).The weight applied to each arc in the network of recurrences is derived empirically by using a total count of record breaking events between the corresponding earthquake zones and the distance involved (Cavers and Vasudevan, 2015;Vasudevan and Cavers, 2013).Each recurrence from an earthquake a to an earthquake b in the sequence is given a weight between 0 and 1, with a weight equal to 1 if the distance between a and b is less than 50 km.If the distance is r with r > 50 km and earthquakes a and b occur in Zones j and k respectively, a weight of is given, where L j k (r) defined by Cavers and Vasudevan (2015) is the number of record-breaking events from zone j to zone k at distance at most r in the network of recurrences.The function in Eq. ( 3) is a decreasing function in r giving a weight close to 0 when the distance r is large.Note that for r = 50 km, an output of 1 is given while for r = 20 000 km, an output of 0 is given.As described by Cavers and Vasudevan (2015), a Markov chain with the inclusion of spatio-temporal complexity of recurring events is derived by summing the weights of the recurrence arcs corresponding to occurrences from state i to state j in consecutive time-intervals.Here, we calculated the time-series of the resulting state-to-state sequence (Fig. 2a) and the corresponding transition probability matrix (Fig. 2b).Since it is obtained from the non-linear, non-stationary global earthquake sequence, we consider it non-linear and non-stationary as well, and hence, Introduction

Conclusions References
Tables Figures

Back Close
Full can be subjected to analysis methods.Although it is not shown here, the approach equally applies to earthquake catalogues from localized seismogenic zones.

Analysis methods
Each sample in the time-series shown in Fig. 2a represents a "zone-configuration" state (Table 2).By definition, a zone-configuration has no zone or some zones or all zones highlighted by an earthquake or more in the optimally chosen time-interval.Going from one sample to the next does not only represent going from one state to the next but also shows the amplitude fluctuation between them.The adjacent states could represent the same zone-configuration or different zone-configurations.The time-series deduced from using the present approach with the five-zones marks the state-to-state fluctuations arising out of the fluctuations of oscillations or earthquake occurrences in the five-zones.We present in the following two analysis methods to glean an insight into the characteristics of the time-series.

Ensemble empirical mode decomposition as applied to state-to-state transition frequency sequence
For non-linear and non-stationary time-series, the method of empirical mode decomposition (EMD) has been recently proposed as an adaptive time-frequency analysis method (Huang et al., 1998(Huang et al., , 1999)).Since the process that leads to the state-tostate transition frequency sequence or time-series is inherently non-linear and nonstationary, it is appropriate to apply the EMD to this data to understand the behavior of the intrinsic mode functions.The time-series (Fig. 2a) reveals the fluctuations in the state-to-state transition frequencies arising out of varying occupancy of the zones from one time interval to the next.A situation would easily arise when two or three successive state-to-state transitions do not have earthquake occurrences in any of the zones studied.This would translate into intermittency in the time-series.Recent studies Introduction

Conclusions References
Tables Figures

Back Close
Full  (Flandrin et al., 2004(Flandrin et al., , 2005;;Gledhill, 2003;Wu andHuang, 2004, 2009) support the idea of carrying out noise-added analyses with the EMD.The noise added analyses involves multiple realization of added noises to the time-series in question, leading to the ensemble EMD (EEMD), as proposed by Wu and Huang (2004Huang ( , 2009)).
In the EEMD, the signal or the time-series in question with the added white noise, denoted as one trial, would populate the whole time-frequency space uniformly with the constituting component of different scales.Since the noise added in each trial is different, the ensemble mean of the noise cancels out and, hence, the signal resides in the intrinsic mode functions generated from the EEMD (Wu and Huang, 2009).
The time-series of state-to-state transition frequencies of the modified Markov chain model, x sstf (t), is taken as the signal.In each realization of the experiment, white noise, w(t), is added to the signal.One might interpret the added white noise as the possible random noise that would be encountered in the measurement process or in certain restrictions applied to the calculation of edge weights in the modified Markov chain.So, for the i th realization, For each realization, we decompose the data with the added white noise into intrinsic mode functions (IMFs).We consider the ensemble means of the IMFs of the decompositions as the final result.Wu and Huang (2009) recommended that the ensemble size should be kept large and the amplitude of the added noise should not be small.We set the ensemble number for the number of realizations in EEMD large such that the noise series cancel each other in the final mean of the corresponding IMFs.For the two parameters, we used an ensemble size of 1000 and added noise with an amplitude of 0.2 times the SD of the original data.We assume that the IMFs resulting from the EEMD truly represent the true IMFs.EEMD results are summarized in Fig. 3a-t with each intrinsic mode function followed by its state-to-state transition frequency matrix.The Hilbert-Huang amplitude spectrum of the time series, shown in Fig. 4, reveals at least two important 407 Introduction

Conclusions References
Tables Figures

Back Close
Full features: (1) the temporal fluctuations in amplitudes occur in packets, each packet containing a set of zone to zone interactions.The oscillatory behaviour of packets contains certain periodicity within the earthquake sequence.
(2) The frequency-dependence of amplitude packets encapsulates the relative importance of the interaction among multiple zones over different time intervals.

Evaluation of fractality in a state-to-state transition frequency sequence
Earthquake occurrences have been modelled to be stochastic point processes (Thurner et al., 1997;Telesca et al., 2001Telesca et al., , 2005Telesca et al., , 2009Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2012).One representation of the point process is to examine the inter-event time-intervals.The resulting inter-event interval probability density function says something about the behavior of the times between events.We do not know anything about the information contained in the relationships among these items.Since successive events do not occur in constant time-intervals, another representation of a point process is given by dividing the time-axis into equally spaced contiguous counting windows of duration τ, and producing a sequence of counts that fall within each time-window.For example, for the kth time-window, the expression for the number of counts, N k (τ), is given by where N k (τ) is the number of earthquakes in the kth window (Fig. 5a-d).The correlation in the process {N k (τ)} is the correlation in the underlying point process (Lowen and Teich, 1993a, b;Thurner et al., 1997;Telesca et al., 2001Telesca et al., , 2005Telesca et al., , 2009Telesca et al., , 2011) ) have accessed such a representation of the point-processes to underscore the existence or non-existence of fractality in them.They have two calculable measures, Fano factor (FF) and Allan factor (AF), to quantify the fractality of the process (Lowen and Teich, Introduction

Conclusions References
Tables Figures

Back Close
Full  1993a, b;Thurner et al., 1997;Telesca et al., 2001Telesca et al., , 2005Telesca et al., , 2009Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2012).The Fano factor is a measure of correlation over different timescales (Thurner et al., 1997).It is defined as the ratio of the variance of the number of events in a specified counting time τ to the mean number of events in the counting time, as is given by where denotes the expectation value.Lowen and Teich (1995) point out that the FF of a fractal point process follows a power law with the power-law exponent, α, obeying 0 < α < 1.In other words, the FF is always greater than 1.For Poisson processes, the FF is always near unity for all counting times, and the fractal exponent is approximately equal to zero.
The Allan factor is a relation with the variability of successive counts (Allan, 1996;Barnes and Allan, 1966).It is the ratio of the variance of successive counts for a specified counting time τ divided by twice the mean number of events in the counting time.
The expression of AF is given as Similar to the FF, the AF assumes values near unity for Poisson processes.Telesca et al. (2009Telesca et al. ( , 2011) ) and Flores-Marquez and Valverde-Esparza (2012) have shown the power-law exponent for the AF to be 0 < α < 1.
In this paper, we examine a new representation of the point process with a Markov chain model.For the working model, we compute the state-to-state transition frequencies as described by Nava et al. (2005) and as applied to global seismicity (Vasudevan and Cavers, 2012;Cavers and Vasudevan, 2013).An expression similar to Eq. ( 7) can 409 Introduction

Conclusions References
Tables Figures

Back Close
Full be derived if we know the optimal time-interval for the Markov chain model.Since we know the optimal time-interval, we introduce a sequence of state-to-state transition frequencies, {N sstf,k (τ)}, with N sstf,k (τ) referring to the number of state-to-state transitions over the kth window for the optimal time-interval, as is shown in Fig. 5e.
There are a few observations to be made.First, N sstf,k (τ) is not necessarily an integer number for any kth window.Following the definition of a state, in the context of a directed graph of a Markov chain model, a state-to-state transition refers to an edge of a graph.It is the weight associated with the edge of the directed graph that plays an important role.Since we have used a modified Markov chain model which includes the influence of the event recurrences in the record-breaking sense, the above expression includes their weights as well in the computation of N sstf (k).The sequence of state-tostate transition frequencies, {N sstf (k)}, yields a time-series.This time-series is the new expression of the point-process where the weighted edges of directed graph of the modified Markov chain represent the significance of the earthquakes between states.This new alternative representation signifies the behavior of the state-to-state transition frequencies over a large time window.Here, seeking to find the time-correlative behavior of the time-series would be of great importance since this would give us an opportunity to see the interaction of zones considered in a collective sense.
Here, we seek to understand the correlative behavior by looking at the two statistical measures, FF sstf and AF sstf , as defined below: The behavior of the two measures, FF sstf and AF sstf , with respect to the optimal timeinterval should shed some light on the correlative behavior of the time-series but also Introduction

Conclusions References
Tables Figures

Back Close
Full on the selective clustering of the certain state-to-state transitions.We consider this knowledge to be useful for forecasting purposes.
In our adaptation of the sum of edge weights for the state-to-state transition frequencies as a new representation of a point-process embedded in the modified Markov chain here, the arguments of Thurner et al. (1997) and Telesca et al. (2001Telesca et al. ( , 2005Telesca et al. ( , 2009Telesca et al. ( , 2011) ) would apply.This means that the FF of the modified Markov chain sequence would follow a power-law with the power-law exponent, α, satisfying 0 < α < 1.
Extending this to FF sstf and AF sstf , as is shown in Fig. 6, we find that the power law exponent calculated, corresponding to the least-squares fit of the data is greater than zero (0.27 and 0.30, respectively).They suggest not only the fractality of the modified Markov chain sequence for optimal time-interval but also the deviation from the Poissonian behavior of earthquake sequencing considered in this present study.
4 Discussion and conclusions Thurner et al. (1997) pointed out that the sequence of counts, generated by recording the number of events in successive counting time-windows of certain length, contained information about the point process depicted by the set of event times.This idea was further tested in understanding the dynamics of earthquake sequencing (Telesca et al., 2009(Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2013), and in particular, the fractal behavior of the sequence of counts.We know that this idea was initially restricted to the sequence of counts for varying windows of interval-times.
In our description of the directed graph of the Markov chain model of any earthquake sequencing, regional or global, we stress the significance of the state-to-state transition probabilities for multiple zones that span the sequence of earthquakes over an optimal time window (Cavers and Vasudevan, 2013;Vasudevan and Cavers, 2013).In other words, the edges of the directed graph carry weights.We conjecture that these weights represent a new definition of the point process.Furthermore, a consideration of the earthquake recurrences within each zone and among zones, following the con-Introduction

Conclusions References
Tables Figures

Back Close
Full cept of recurrences in the record-breaking sense (Davidsen et al., 2008), leads to an empirically-determined distance-dependent weights for the edges.Unlike extending the idea of the sequence of counts where every event occurrence augments the counting value by unity (Thurner et al., 1997;Telesca et al., 2009Telesca et al., , 2011;;Flores-Marquez and Valverde-Esparza, 2013), we consider the summing of the weights for each edge such that the sum represents a "pulse" for each state-to-state transition.We analyse the resulting time-series from the point of view of its Fano factor and Allan factor.There is evidence for fractality of the multi-state modified Markov chain to represent the earthquake sequencing, as is revealed by the power-law scaling behavior present in the Fano and Allan factors with their respective exponents of 0.27 and 0.30.Cavers and Vasudevan (2013) interpreted the Markov chain of 32-states for fivedistinctly different zones to contain the basic combinatoric structure superimposed by the thumb-print of the undulatory structure of the recurrence weights.Since the earthquake sequencing is in general non-linear and non-stationary, we contend that the time-series representing the above Markov chain is also non-linear and non-stationary, and is conducive to an ensemble empirical mode decomposition (EEMD) procedure to understand its intrinsic mode functions (IMFs).The ensemble empirical model decomposition of the time-series leads to nine intrinsic mode functions and a trend.Each one of the IMFs reveals the amplitude fluctuation of the state-to-state transitions.While there is a commonality in the relative dominance of the subduction-style earthquakes, represented by the top right corner grid of the transition frequency matrices, the presence or absence of certain state-to-state transitions in certain IMFs reveals the importance of integral multiples of the optimal time-interval.
A simple observation of the first 6 or 7 IMFs stresses the importance of multiple-zone approach to global seismicity problem in that the earthquake sequencing for the time period we considered has similar oscillatory behavior of the state-to-state transition probabilities from the point of view of the amplitude scaling and the oscillating period.The growth and decay of oscillations in easily identifiable packets in each IMF following

Conclusions References
Tables Figures

Back Close
Full have introduced the Markov chain model to study the earthquake sequencing in a seismogenically active region where the region is partitioned into zones.The functionality of the method is determined by Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | of binary digits b 4 b 3 b 2 b 1 b 0 , where b L = 1 indicates an earthquake occurrence in zone L during the specified time interval ∆t, and b L = 0 indicates the lack of an earthquake Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 .Figure 5 .
Figure 1.A graph representation of earthquake sequencing with arcs (with weights w i j ) representing transitions between states.