Dynamical Features of the Seismicity in Mexico by Means of the Visual Recurrence Analysis

1. Introduction

The recurrence, based in the recurrence Poincare theorem, is a fundamental property of dynamical systems that has been exploited to characterize the behavior of dynamical systems in the phase space. Recurrence is defined when an orbit visits a region of phase space that was previously visited [1]. In this context, the so-called recurrence plot (RP), introduced by Eckmann et al in [2], is a powerful tool for the visualization and analysis of the underlying dynamics of the systems in the phase space like determinism, divergence, periodicity, stationarity among others, for instance the lengths of the diagonal line structures in the RP are related to the positive Lyapunov exponent [3]. Methods based in RP have been successfully applied to wide class of data obtained in physiology, geology, physics, finances and others. RP are especially suitable for the investigation of rather short and non-stationary data [4], and complex systems [5]. For deterministic systems the analysis in the phase space is relatively direct because their solutions are represented as time series, nevertheless, for real natural systems like clime, atmosphere, dimetilsulphure production [6], some of their dynamical variables are gathered as punctual processes. Particularly, the representation of a seismic sequence as a time series is one of the most debated questions in Geophysics, nevertheless natural time analysis, introduced by in [7] by Varotsos et al considers sequences of events and, with this methodology has been possible to identify signals with noises [8] and characterize seismic processes from an order parameter properly defined [9]. Taking in consideration that faults and tectonic plates can be considered as dynamical systems which permit to apply techniques derived from the dynamical systems theory and nonlinear analysis like recurrence plots, clustering, entropy, fractality, correlation, memory among other. In this work we studied properties of the seismicity occurred in Mexico in the period 2006-2014, reported by Servicio Sismológico Nacional, (www.ssn.unam.mx), by considering the occurrence of events as sequences of magnitudes (time series) and determining the recurrence plots. Based on the Visual Recurrence Analysis (VRA), it is possible to get dynamical features of the seismicity. The most important seismic region in Mexico is located along the Mexican Pacific coast. In [10] was proposed a division of 19 regions of the seismo tectonic zone. They took into account seismic characteristics of the existing catalogs for the seismicity occurring in Mexico from 1899 to 2007, for details see [9] and references therein. In order to characterize the Mexican seismicity as a dynamical system driven by the tectonic movements in the Pacific edge, we proposed to investigate the qualitative dynamical features for the dispersion zone, Baja California, and for the subduction zone formed by La Ribera and Cocos plates subducting beneath the Northamerica plate. Taking into account the Geophysical structure, the seismic activity and the b-values in the Gutemberg-Richter law, in this work, six selected regions have been considered. These six regions are named: Baja California (Z1), Jalisco-Colima (Z2), Michoacán (Z3), Guerrero (Z4) and Oaxaca (Z5) and Chiapas (Z6), the first one corresponds with a dispersion zone and the other ones are subduction regions, [11]. The analyzed data set in this work correspond with the Mexican catalogue which is complete for magnitudes greater than 3, within the mentioned period. Since a geophysical point of view, the seismic activity should be individually characterized in each region because the underlying dynamics shows particular features, as is described in the next section. First, seismic data are represented as a temporal sequence of magnitudes, then, the phase space is reconstructed, by estimating the time delay and the embedding dimension. In order to distinguish some features of the underlying dynamics of each Mexican region seismic, the aim of this work is to study the recurrence plot behavior based on the visual recurrence analysis, taking into account the sequence of events (magnitudes) in time and, on the other hand, analyzing the inter-events time series. Our analysis shows important differences in the recurrence maps of each region. Our finding suggest that the patterns obtained could be associated with the local geophysical structures of each subduction and dispersion zones driven by their characteristic nonlinear dynamical features of each region.

2. Methodology

2.1. Phase space reconstruction

To reconstruct the dynamics of the system is necessary to convert the time series in state vectors, being the principal step the phase space reconstruction. Takens [12] showed that it is possible to recreate a topologically equivalent picture of the original multidimensional system behavior by using the time series of a single observable variable. The idea is to estimate a time delay τ, and an embedding dimension m. The parameters, m and τ, must be properly chosen, although there exist some algorithms to do that, appropriated and tested methods are the Mutual Information function to obtain the time delay and, to get the embedding dimension, the False Nearest Neighbors. Once the time delay and the embedding dimension have been approached, the phase space can be reconstructed. For a time series of a scalar variable

$$x\left(t_i\right),i=\mathrm{1,2},3,\cdots ,N$$E1

the next vector is construct in phase space

$$Y\left(t_i\right)=\left[x\left(t_i\right),x\left(t_i+\tau \right),x\left(t_i+2\tau \right),\cdots ,x\left(t_i+\left(m-1\right)\tau \right)\right]$$E2

where i = 1 to N – (m – 1)τ, τ is time delay, m is the embedding dimension and N – (m – 1)τ is number of states in the phase space. According with the embedding theorem [12] dynamics reconstructed using this formula is equivalent to the dynamics in the phase space in the sense that characteristic invariants of the system are conserved.

2.2. Mutual information function

Formally, the Mutual Information is defined, for two stochastic variables X and Y, as I(X;Y) = H(Y,X) – H(X|Y) – H(Y|X) where H(Y,X) is the joint entropy, H(X|Y) and H(Y|X) are the conditional entropies. If X represents the sequence x(t_i) and Y the respective sequence x(t_i + τ) the Mutual Information Function (MIF) is able to calculate the global correlation in a time series taking into account the linear and non linear contributions, being MIF the most frequently used to calculate the time delay and described as follows:

$$\begin{array}{c}x\left(t_i\right)\\ \\ x\left(t_i+\tau \right)\\ P\left(\right)\\ \\ P\\ \frac{P\left(x\left(t_i\right),x\left(t_i+\tau \right)\right)}{}\\ P\left(x\left(t_i\right),x\left(t_i+\tau \right)\right){log}_2\\ I\left(\tau \right)={\sum }_{i=1}^N\end{array}$$E3

The mutual information I(τ) is a measure of the relative entropy between the joint distribution and the product of distributions P(x(t_i)) and P(x(t_i + τ)), where P(x(t_i), x(t_i + τ)) is the joint probability of the signal measured between the times t_i and t_i + τ. The expressions P(x(t_i)) and P(x(t_i + τ)) indicates the and marginal probabilities. MIF is a nonlinear generalization of the linear autocorrelation function. According to [13] the suitable value of τ is attained with the first local minimum of I(τ). When the time series is uncorrelated or random, like white noise, the next equality holds

$$\begin{array}{c}x\left(t_i\right)\\ \\ x\left(t_i+\tau \right)\\ P\end{array}$$E4

and I(τ) = 0 [14]. Typical cases are white noise, periodic and periodic + noise and Rossler time series. In Figure 1 the MIF behavior for this typical cases are depicted.

Figure 1.

MIF behavior for white noise, periodic + noise are uncorrelated signals with a time delay τ = 1. For periodic and Rossler τ = 3 and 13 respectively.

2.3. False nearest neighbors

False Nearest Neighbors method (FNN) is an iterative algorithm which estimates the minimal embedding dimension of the system proposed by [15]. The idea is based in the uniqueness theorem of a trajectory in the phase space. A nearest neighbor P of a point Q in a phase space of d-dimensional embedding is labeled false if these points are close only due to the projection from a higher-dimensional (d+1)-dimensional phase space. Thus, the FNNs will separate if the data is embedded in (d+1)-dimensional space, while the true neighbors will remain close. When all the FNNs have been detected, then the minimal sufficient embedding dimension can be identified as the minima dimension needed to achieve zero fractions of the FNNs being the embedding dimension required. For each vector Y(i), its nearest neighbor Y(j) is looked in a m-dimensional space. In order to do a comparison, the distance d((Yi).Y(j)) is calculated. By iteration of points, the condition:

$$M_i=\frac{‖Y\left(i+1\right)-Y\left(j+1\right)‖}{‖Y\left(i\right)-Y\left(j\right)‖},$$E5

If M_i exceeds a threshold M_th, this vector is marked as a nearest neighbor. When the fraction M_i of vectors that satisfy the condition M_i >M_th, tends to zero the embedding dimension is attained in this case. The FNN is exemplifying with a simple case in Figure 2.

Figure 2.

In 1D, B and B’ seem to be nearest neighbors of A, nevertheless in 2D, B’ is a false nearest neighbor.

For the mentioned examples, the false nearest neighbor fraction is calculated as an embedding function and is depicted in the Figure 3. To get the embedding dimension, first was estimated the time delay for each time series by evaluating the respective MIF. For white noise and periodic with noise time series, the embedding dimension is high and for periodic and Rossler signals, the dimension is short.

Figure 3.

The fraction of false nearest neighbors, as a function of the embedding dimension, is depicted, for white noise (black), periodic (green), periodic with noise (blue) and Rossler (red) time series.

2.4. Recurrence Plot (RP)

In [2] Eckmann et al introduced the so-called Visual Recurrence Analysis (VRA) based on a graphical method designed to locate hidden recurrent patterns, nonstationarity and structural changes observed in the phase space of a dynamical system. The aim of the Recurrence Plot (RP) method is to visualize the recurrences of dynamical systems as a function of time. A brief description is as follows: assuming an orbit of the system in the phase space $\left\{{\overrightarrow{x}}_1,{\overrightarrow{x}}_2,{\overrightarrow{x}}_3,\cdots ,{\overrightarrow{x}}_N\right\}$. Each vector of this trajectory corresponds with a state of the system where their components can represent physical quantities, for example, position and velocity for mechanic systems like a pendulum or, pressure, volume and temperature for thermodynamical states. Then, in order to get the RP, the recurrence matrix must be constructed which is defined as [5]:

$$R_{i,j}=\left\{\begin{array}{c}1if‖{\overrightarrow{x}}_i-{\overrightarrow{x}}_j‖\le \epsilon \\ 0if‖{\overrightarrow{x}}_i-{\overrightarrow{x}}_j‖>\epsilon \end{array}\right.i,j=1,\cdots ,N$$E6

were $\epsilon >0$ is an error. Roughly speaking, the matrix elements represent the distance between two vectors in the times i and j. Once the RP has been obtained from the recurrence matrix, a quantification of the features can be done, for example from periodic patterns to chaotic behavior [5, 16]. As it is well known, in an experiment only a sequence of scalar values can be measured and it is assumed that the information is available on an univariate time series, which is part of a larger n-dimensional (maybe deterministic) model.

Figure 4 shows examples of RP of the four cases described above: white noise, periodic, periodic with noise and chaotic systems.

Figure 4.

In this panel the RP of the cases described above: a) White noise, b) periodic, c) periodic with noise and d) Rossler. Can be observed the different structures of the four examples described above.

3. Seismic regions and data set

We analyzed the whole seismic catalog of the Mexican SSN, (www.ssn.unam.mx) from 2006 to 2014 Due to geophysics features of the Mexican subduction zone, it has been described in [9 and references [20,21,24] therein], sugesting that it can be studied in segments where the six regions are: Baja California (Z1), Jalisco–Colima (Z2), Michoacán (Z3), Guerrero (Z4), Oaxaca (Z5) and Chiapas (Z6), in Figure 5 the six regions are showed.

Figure 5.

The six regions are showed in this map. The region Z1 is a dispersion and the regions from Z2 to Z6 the seismicity is driven by subduction regions.

The following panel of Figure 6 displays the seismicity within 2006-2014 periods for each region. It can be observed the seismicity monitored in Region Z1, Peninsula of Baja California where the Pacific–North America plate boundary in southern California and the north side of Baja California peninsula. The seismicity in Region Z2, Jalisco-Colima is a subduction zone located to the west, where the Rivera plate subducts at a steep angle. In the Region Z3 Michoacán, the dip angle of the Cocos plate decreases gradually toward the southeast. For the region Z4, Guerrero is bounded approximately by the onshore projection of the Orozco and O’Gorman fracture zones, the subducted slab is almost horizontal and moves under upper continental plate. The regions Z5, Oaxaca and Z6, Chiapas are located in the southeastern of Mexico, the dip of the subduction zone gradually increases to a steeper subduction in Central America.

Figure 6.

Seismicity monitored in the six regions: Z1 Peninsula of Baja California. Z2 Jalisco-Colima. Z3 Michoacán. Z4 Guerrero. Z5, Oaxaca. Z6 Chiapas.

As is showed in Figure 6, the number of earthquakes reported by the SSN in each region is different. Firstly, the seismicity in the region Z1, peninsula of Baja California, shows two periods of time, the fist one from 2006 to 2010 and the second one from 2011 to 2014, in the first period little seismic activity compared to the second is observed, however, a sudden change in the seismic activity is observed, this behavior in the seismic activity before and after 2010 can be observed in the regions Z2 and Z3 corresponding with Jalisco-Colima and Michoacán respectively This situation may be due to system upgrades monitoring stations. As mentioned above, the Z1 region evolve with a process in which the peninsula of Baja California is separated from the continental plate, while in the other regions, the dynamics is driven by subduction between continental plate and plates Rivera and Cocos. Specifically, Jalisco-Colima region (Z2) the subduction is given between La Rivera and The north-America plates where the stress and strain fields determine the direction of movement in which La Rivera subducts being different from the case of the Cocos plate. According with the catalogue, Guerrero (Z4), Oaxaca (Z5) and Chiapas (Z6), are the regions with a high seismic activity. It can be observed that in all regions have occurred earthquakes with magnitudes M ≈ 7. By considering this division of six regions and the seismic activity showing seismic clustering, the b-values in the G-R law were recalculated in each case by using the Aki model [17]:

$$b=\frac{log\left(e\right)}{⟨M⟩-M_c},$$E7

where <M> is the average magnitude and M_c is the completeness magnitude of the seismic sequence that represents the minimum magnitude over which the frequency-magnitude distribution behaves as the power-law, N~10^-bM [18]. The b-values calculated for each region are depicted in Figure 7 and resumed in Table 1.

REGION	<M"/>	b-Value
Baja California (Z1)	3.93	1.32
Jalisco- Colima (Z2)	4.20	1.45
Michoacán (Z3)	3.98	2.41
Guerrero (Z4)	3.91	2.07
Oaxaca (Z5)	3.92	1.97
Chiapas (Z6)	4.07	1.61

Table 1.

The b-values calculated for each region.

Figure 7.

b-value calculation from the Gutenberg-Richter law for the six regions. In Table 1, the values are resumed.

Because geophysics shaping the seismic zone as well as the various mechanisms and processes that take place in each of the regions, the statistical properties associated with the seismicity of each of these are reflected as different values to the parameters of the laws scaling as in the case of the Gutenberg-Richter law. This situation indicates that the local stress fields and stress must drive the interactions between different parts of a complex system, such as the Earth’s crust. So, through a study in the context of dynamic systems, where the seismic activity is considered as a response to the underlying dynamics is possible to observe different characteristics of the system and not observed directly from a statistical point of view. Although seismicity is considered, as a sequence of events whose measurable variable is the magnitude, cannot be put aside the temporal component, that is, the distribution of interevent defined as the time between two earthquakes within a region while. It is noteworthy that there are time series studies interevent by multifractal analysis for seismicity in Italy [19].

The Figure 8 shows the interevent time series of the seismicity studied. For this time series the behavior of the MIF and FNN fraction is calculated.

Figure 8.

Interevent time series associated with the seismic activity observed in 2006-2014 period

4. Results

In this work six seismic sequences and their respective interevent time series were analyzed by means of the RP method. Firstly the phase space is reconstructed for each sequence. The MIF and FNN fraction for the magnitudes sequences are depicted in Figures 9 and 10. In order to reconstruct the phase space, the time delay τ, and the embedding dimension m, were estimated by taking the first minima of the mutual information function and the False Nearest Neighbors algorithms respectively. Then the recurrence matrix is obtained whose matrix elements, R_{i j}, are the distances D_ij, between states Y(t_i) and Y(t_j) in the reconstructed phase space, to calculate D_ij, the Euclidean norm was chosen. The τ–values and m-values are contained in the Table 2. Once the recurrence matrix is obtained the distribution of distances between states computed, in Figure 12 the probability distribution function (pdf) for distances is showed. The r_max-value of the pfd for Z2, Z4, Z5 and Z6 are located around r_max ≈1.2 and r_max < 1 for Z1 and Z3. As has been mentioned, Z1 is a region where the peninsula is separating from the continental plate and Z3 is located where the border between La Rivera and Cocos plates are in contact and subducted into the continental plate. Regarding the shape of PDF in all cases an exponential tail is observed.

REGION	τ	m
Baja California (Z1)	5	13
Jalisco- Colima (Z2)	2	12
Michoacán (Z3)	3	7
Guerrero (Z4)	3	9
Oaxaca (Z5)	9	12
Chiapas (Z6)	4	11

Table 2.

The τ–values and m-values for each region for the case of seismicity sequence.

From the behavior of MIF is possible to identify the correlation in seismic sequences. According to this behavior, the region with the lowest correlation is Z2 (Jalisco-Colima) and the highest correlation is Z5 (Oaxaca), as is showed in Figure 9. However, the behavior of the FNN fraction as function of the embedding dimension of the seismic sequences indicates that all of them are in phase spaces of high dimension being the shorter Z3 with m = 7.

To calculate the recurrence matrix, the distances between pairs of vectors were computed in phase space. By definition this matrix is symmetric because D_ij = D_ji, and the principal diagonal D_ii = 0. The Recurrence Plot is the graphical representation of the recurrence matrix. The examples of RP depicted in Figure 4, were drawn in black and white because the recurrence matrix was built according to definition given in Eq. (1). In order to take into account all distances computed and perhaps their distribution (Figure 12), a color code allows characterizing qualitatively some features of the dynamics as are displayed, in the panel of Figure 11, the RP of the six seismic sequences. According [5, 20, 21], RP of some cases of different topologies can be distinguished, for instance: Stationary systems display homogeneous RP like white noise, for periodic or cuasi-periodic systems appears recurrent structures as diagonal lines and checkerboard forms, for non-stationary systems drifts are present, abrupt changes in the systems indicates extreme events, vertical (horizontal) lines represent time intervals where a state remains constant or changes very slowly. In general, the RP of the seismicity occurred from 2006 to 2014 are displayed in Figure 11 where, in a first glance, typical patterns, like periodicity or cuasi periodicity and white noise are not observed. Nevertheless, clusters bordered for vertical and horizontal lines are present suggesting slow changes in the system. The color distribution and the clusters in RP of BC(Z1), JC(Z2), and Ch(Z6) suggest drifts indicating non stationary dynamics associated possibly with their geophysical features because in BC, the Pacific–North America plate boundary in southern California and the north of Baja California peninsula where many faults are connected in a complex geometrical pattern, continuing into a divergent tectonic plate in the Gulf of California. In Jalisco–Colima region, the Rivera plate subducts at a steep angle plate in Central America and for Ch, the Cocos plate subducts beneath the coast but two perpendicular faults, Clipperton and Tehuantepec, contribute wlth their local dinamical evolution. More similar RP structures are dysplayed in M, G and O, which seems to show more stability because the respectives RP are more uniform which could be indicating that the dynamics is driven by the interacion between the Cocos plate which subducts in the same direction beneath the North America plate and the dip angle of the Cocos plate decreases gradually. Our findings are consistent with the results reported in [11] where an analysis of non extensive model of the similar regions were studied indicating that JC region is the most unstable seismic zone in Mexico.

Figure 9.

MIF behavior for the magnitude seismic sequences of the six regions.

Figure 10.

The fraction of false nearest neighbors, as a function of the embedding dimension, is depicted for all regions.

The distribution function associated with the distances between states in phase space allows us to observe the most likely value, and suggest a other possible criterion for determining characteristics distances among all states in the phase space and to determine the ε-value in the definition of the matrix recurrence (Ec. 7). The Figure 12 shows the pdf of the distances for the six sequences of seismic magnitudes. It is observed that the maximum values of the pdf are located in different possitions and possibly this behavior could be associated with the geophysics features of the regions: Z2, Z4, Z5 and Z6 are subduction zones where Z3 is determined by the relative motion between the Rivera plate and the Continental plate, while the other three are determined by the interaction of the Cocos plate subducting under the continental plate. For Z1 the seismic activity is produced by the movement of separation between the peninsula and the continental shelf of North America and the Z3, the seismic activity is determined by the interaction of the plates Rivera in contact with the plate subducting Cocos and both the continental plate.

On the other hand, the results of the interevent time series are showed in Figure 13 for MIF behavior and Figure 14 for FNN fraction. In Table 3 the τ –values and the embeeding dimension m-values are presented. In contrast with the seismic sequences of magnitudes, the interevent time series the correlations are similar, nevertheless the FNN fraction decrese almost monotonically indicanting high embedding dimension. This behavior is similar with the example of periodic with noise time series suggesting that the inerevent time series could be described with a possible detrministic model plus a stochastic process.

Figure 11.

Recurrence Plot of Baja California and Jalisco-Colima (upper). Michoacan and Guerrero (middle), Oaxaca and Chiapas (below)

Figure 12.

Probability distribution function of the distance between states in the phase space for the six regions

Figure 13.

MIF for sequences of interevents.

REGION	Interev	interev
	τ	m
Baja California (Z1)	4	20
Jalisco- Colima (Z2)	5	20
Michoacán (Z3)	4	13
Guerrero (Z4)	5	20
Oaxaca (Z5)	4	20
Chiapas (Z6)	4	20

Table 3.

The τ–values and m-values for each region, in this case for interevent sequence.

Figure 14.

The fraction of false nearest neighbors, as a function of the embedding dimension, is depicted for all regions, now for the interevent sequence.

5. Concluding remarks

It is well known that the Mexican Pacific is an important seismic region where large earthquakes have occurred with devastating consequences producing significant economic downturns and especially many human losses. Due to interactions between subduction zones and the slow separation of the Baja California Peninsula, this region is considered a complex system that evolves as consequence of many processes that occur in the interior of the Earth as well as in the areas of contact between the surfaces involved. In this context, [10] proposed a division of 19 regions of the seismo tectonic zone taking into account seismic characteristics of the existing catalogs for the seismicity occurring in Mexico from 1899 to 2007 and a seismic historic compilation from previous publications and of some catalogs. In order to distinguish some features of the underlying dynamics of each Mexican region seismic, the aim of this work is to study the recurrence plot behavior based on the visual recurrence analysis, taking into account the sequence of events (magnitudes) in time and, on the other hand, analyzing the inter-events time series. Our analysis shows important differences in the recurrence maps of each region. In a similar way by considering the seismicity monitored by SSN within the period 2006-2014, and identifying clusters of earthquakes that can be associated with the geophysical features of the Mexican Pacific, six regions were considered to study (Figure 5). In this paper we studied dynamical features of the six seismic regions located along Mexican Pacific coast. The analyzed data set corresponds with the Mexican seismic catalogue reported by the SSN. First, sequences of magnitude of earthquakes and the interevent time series were studied. Their analysis was performed by means of the phase state reconstruction and the RP of each region. Our findings indicate a possible correlation between the RP calculated and the geophysical features characteristics in each zone (panel of Figure 11). RP displayed of BC, JC and Ch show non periodicities, correlation (not white noise structure), non stationariety. For M, G and O, RP are more similar and stability is observed. The results for the interevents time series, short correlation and large embedding dimension, suggest the possibility to establish a combination between a detremiistic model plus a stochastic noise. Our finding suggest that the patterns obtained could be associated with the local geophysical structures of each subduction and dispersion zones driven by their characteristic nonlinear dynamical features of each region.

Acknowledgments

This work was supported by the Irreversible Processes Physics Research Area of Universidad Autónoma Metropolitana, México. ARR, LRMT and RTPH thanks to Basic Sciences Department of UAM. IRR thanks UPIITA-IPN México. ARR thanks the bilateral Project CNR (Italy) and CONACyT (México).