Seismicity characterization of oceanic earthquakes in the Mexican territory

We analyzed the seismicity of oceanic earthquakes in the Pacific oceanic regime of Mexico. We used data from the earthquake catalogues of the Mexican National Service (SSN) and the International Seismological Centre (ISC) from 1967 to 2017. Events were classified into two different categories: intraplate oceanic (INT) and transform fault zone and mid-ocean ridges (TF-MOR) events, respectively. For each category, we determined statistical characteristics such as magnitude frequency distributions, the aftershocks decay rate, the nonextensivity parameters, and the regional stress field. We obtained b values of 1.17 and 0.82 for the INT and TF-MOR events, respectively. TF-MOR events also exhibit local b-value variations in the range of 0.72–1.30. TF-MOR events follow a tapered Gutenberg–Richter distribution. We also obtained a p value of 0.67 for the 1 May 1997 (M_w=6.9) earthquake. By analyzing the nonextensivity parameters, we obtained similar q values in the range of 1.39–1.60 for both types of earthquakes. On the other hand, the parameter a showed a clear differentiation, being higher for TF-MOR events than for INT events. An important implication is that more energy is released for TF-MOR events than for INT events. Stress orientations are in agreement with geodynamical models for transform fault zone and mid-ocean ridge zones. In the case of intraplate seismicity, stresses are mostly related to a normal fault regime.

1 Introduction

Mid-ocean ridges and transform fault zones are two of the main morphological features of oceanic environments. Most of the oceanic earthquakes take place in areas close to the active spreading ridges where the seismogenic zone is narrow. For this reason, large aspect ratios are often required to generate moderate-size strike-slip oceanic earthquakes. Nevertheless, the rupture process of oceanic events is still poorly understood. Previous studies showed that these types of events have peculiar characteristics. For example, estimates of seismic coupling for oceanic transform faults indicate that about three-fourths of the accumulated moment is released aseismically (Abercrombie and Ekström, 2003; Boettcher and Jordan, 2004), and some oceanic events exhibit slow slip ruptures (Kanamori and Stewart, 1976; Okal and Stewart, 1992; McGuire et al., 1996). Earthquakes that have longer durations than those predicted by scaling relationships are considered as slow (Abercrombie and Ekström, 2003). These “slow” ruptures are mainly interpreted as having low rupture velocities. On the other hand, others proposed that the slow ruptures may be explained as numerical artifacts generated by the inversion procedures (e.g., Abercrombie and Ekström, 2001, 2003). Several oceanic strike-slip events were reported as being energy deficient at high frequencies (Beroza and Jordan, 1990; Stein and Pelayo, 1991; Ihmlé and Jordan, 1994) or having high apparent stresses (Choy and Boatwright, 1995; Choy and McGarr, 2002). On another front, oceanic earthquakes also occur as intraplate events but to a lesser extent. The reason is that the oceanic plate interiors do not experience significant strain over long periods of time (Bergman and Solomon, 1980; Bergman, 1986). Oceanic intraplate earthquakes originate from the following processes: stresses of the oceanic crust, in regions that concentrate significant deformation; reactivation of faults; or thermoelastic stresses (Bergman, 1986).

From the statistical perspective, previous studies showed that the magnitudes of the major events in the mid-oceanic ridges and transform fault zones are relatively smaller ($\mathrm{6.0}\le M_{\mathrm{w}}\le \mathrm{7.2}$) compared to continental events. The b value in oceanic environments showed significant variability. For example, Tolstoy et al. (2001) reported high b values (b∼1.5) in the Gakkel Ridge associated with volcanic activity. In the Southwest Indian Ridge, Läderach (2011) found b values of about 1.28. Bohnenstiehl et al. (2008) quantified the b value in the East Pacific Rise, obtaining estimations in the range of 1.10<b<2.50. Global studies have also shown that the mid-ocean ridge transform seismicity follows a tapered frequency-moment distribution (Kagan and Jackson, 2000; Boettcher and McGuire, 2009). Cowie et al. (1993) studied the seismic coupling of mid-ocean ridges. They found that fast-spreading ridges (≥9.0 cm yr⁻¹) are weakly coupled. On the contrary, slow-spreading ridges (≤4.0 cm yr⁻¹) are strongly coupled (Cowie et al., 1993). In Mexico, oceanic earthquakes have been poorly studied. There are no systematic studies of their statistical characteristics. In this article, we characterized the oceanic seismicity in Mexico. We determined the orientation of the principal stresses, the b and p values, and the nonextensivity parameters. The results may help to understand the ocean tectonics, particularly in Mexico.

2 Tectonic setting

The Pacific oceanic regime of Mexico is an active area exhibiting ongoing tectonic plate interactions. These interactions involve the Cocos (CO), Pacific (PA), Rivera (RI), and North American (NA) plates. The Gulf of California and the Middle America Trench (MAT) are separated by the Tamayo Fracture Zone (TFZ). The convergence rate between the RI and NA plates decreases northward along the MAT (averaging about 2–3 cm yr⁻¹ in the RI Plate, which is slower than the adjacent CO Plate, about 5–7 cm yr⁻¹) (NUVEL-1A model, DeMets et al., 1994). Sea-floor spreading takes place along the northernmost segment of the East Pacific Rise in the Cocos and Rivera segments (EPR-CS and EPR-RS, respectively). In the EPR-RS, the spreading rates range from 5.3 cm yr⁻¹ at the northern to 7.3 cm yr⁻¹ at the southern end of the rise (Bandy, 1992). The spreading rates at the EPR-CS are as follows: 7.0 cm yr⁻¹ near the Rivera Fracture Zone (RFZ), 8.2 cm yr⁻¹ near the Orozco Fracture Zone (OFZ), 10.1 cm yr⁻¹ near the Clipperton Fracture Zone (CFZ), and 10.7 cm yr⁻¹ near the Siqueiros Fracture Zone (SFZ) based on the NUVEL-1A model (DeMets et al., 1994; Pockalny et al., 1997). The Rivera Transform (RT) is a left transform fault with fast slipping (∼7.0 cm yr⁻¹) (Bandy et al., 2011) (Fig. 1). Due to these differences in subduction and spreading rates and convergence direction of the RI and CO plates, complex seismicity patterns are generated in this region. In the last century, some intermediate-size earthquakes (6.8<M<7.1) have taken place in the Pacific oceanic regime of Mexico (Table 1 and Fig. 2).

Figure 1 Main tectonic features in the oceanic environment off the Pacific coast of Mexico discussed in the text. CO is the Cocos Plate; NA is the North American Plate; PA is the Pacific Plate; RI is the Rivera microplate; TMVB is the Trans-Mexican Volcanic Belt; TFZ is the Tamayo Fracture Zone; EPR-RS is the East Pacific Rise Rivera segment; EPS-CS is the East Pacific Rise Cocos segment; and RT is the Rivera Transform. Blue triangles are volcanoes. Dashed lines show contour lines of the subducted slab. Arrows indicate the motion of the PA, CO, and RI plates. R1 to R6 are the regions in which the study area was divided for analyzing stress and seismicity characteristics. Red number indicates the slipping rates. Pink numbers indicate convergence rates, and black numbers indicate spreading rates.

Table 1 Major oceanic earthquakes in Mexico (M>6.8).

(1) Data from the Decade of North American Geology Project (DNAG) of the National Geophysical Data Centre (NGDC) and the Geological Society of America. (2) Pacheco and Sykes (1992). (3) ISC earthquake catalogue. (4) Abe (1981). (5) Global CMT catalogue.

3 Data and methods

3.1 Data

We used earthquake catalogues of the Mexican National Service (SSN) and the International Seismological Centre (ISC) from 1967 to 2017. Events with no reported magnitude were excluded from our analysis. Reported magnitudes (based on surface, M_s; body, m_b; and coda, M_c, waves) were converted to moment magnitude (M_w). The SSN reports M_w for events in Mexico. For the case of the ISC events, M_s and m_b were converted to M_w using the scaling relationships of Scordilis (2006). We classified the seismic events into two different categories: (1) intraplate oceanic events (INT, red dots in Fig. 2) and (2) transform fault zone and mid-ocean ridge events (TF-MOR, green dots in Fig. 2). The INT catalogue consists of 177 events with magnitudes in the range of 2.9–6.0. The TF-MOR catalogue includes 2074 earthquakes with magnitudes between 2.7 and 6.9. We also used the Global CMT focal mechanism catalogue (Dziewonski et al., 1981; Ekström et al., 2012) with solutions from 1976 to 2017. For the stress analysis, the focal mechanism catalog was divided into six sub-catalogues shown in Fig. 8 (R1 to R6).

Figure 2 Seismicity in the oceanic environment off the Pacific coast of Mexico from 1899 to 2017. The size of the circles represents magnitude. Brown circles are relevant historical earthquakes shown in Table 1 with M>6.8. Red circles are intraplate oceanic events, and green circles are transform fault zone, and mid-ocean ridges earthquakes. Epicenters were compiled from the Mexican National Service (SSN) and the International Seismological Centre (ISC) catalogues.

3.2 Methods

3.2.1 Moment/magnitude earthquake distributions

The Gutenberg–Richter law describes the earthquake magnitude distribution (Ishimoto and Iida, 1939; Gutenberg and Richter, 1944). Mathematically, this law is expressed by the following equation: log${}_{\mathrm{10}}N\left(M\right)=a-bM$, where N(M) is the cumulative number of earthquakes with a magnitude larger than a given magnitude limit (M), the constant b (or b value) describes the slope of the magnitude distribution, and the constant a is proportional to the seismic productivity. The b value describes the distribution of small to large earthquakes in a sample, and it is considered to be specific for a given tectonic environment (e.g., Scholz, 1968; Wyss, 1973; Smith, 1981; Wiemer and Benoit, 1996; Wiemer and Wyss, 2002). In several tectonic environments, b is close to 1 (Utsu, 1961), with deviations affected by many factors. Among them, high thermal gradients and rock heterogeneity (Mogi, 1962; Warren and Latham, 1970) increase the b values. On the contrary, increments in effective and shear stresses (Scholz, 1968; Wyss, 1973; Urbancic et al., 1992) reduce the b value. The b value differs not only between unrelated fault zones (Wesnousky, 1994; Schorlemmer et al., 2005) but also for specific space and time periods (Nuannin et al., 2012). Schorlemmer et al. (2005) found a global dependence of the b value on focal mechanism, which was corroborated at a regional level by Rodríguez-Pérez and Zúñiga (2018). According to those authors, the highest b values correspond to normal-faulting events, followed by strike-slip and thrust earthquakes, respectively. To characterize the b value of oceanic earthquakes and compare the results with other tectonic environments, we calculated the b value with a robust method which has proven its validity in many studies. We estimated the b value by means of the maximum likelihood formulation of Aki (1965) and the completeness magnitude (M_c) employing the maximum curvature method (Wiemer and Wyss, 2000) with the aid of the ZMAP software package (Wiemer, 2001).

As reported by previous authors, seismicity on the mid-ocean transform faults is better represented by a tapered frequency-moment distribution (e.g., Boettcher and McGuire, 2009). This distribution has the following form (Kagan, 1997, 1999; Kagan and Jackson, 2000; Kagan and Schoenberg, 2001; Vere-Jones et al., 2001):

$$\begin{array}{}\text{(1)}& N\left(M\right){=\mathrm{N}}_{\mathrm{o}}{\left({\displaystyle \frac{M_{\mathrm{o}}}{M}}\right)}^{\mathit{\beta }}\exp \left({\displaystyle \frac{M_{\mathrm{o}}-M}{M_{\mathrm{m}}}}\right),\end{array}$$

where β is one of the parameters to determine, $\mathit{\beta }=(\mathrm{2}/\mathrm{3})b$, where b is the b value, N₀ is the cumulative earthquake number over a completeness threshold seismic moment (M₀), and M_m is the maximum expected moment. We analyzed if this frequency distribution is suitable for describing the seismicity of oceanic events in Mexico. In order to calculate the tapered Gutenberg–Richter distribution, we used the MATLAB function Get_GR_parameters.m developed by Olive (2016). The tapered Gutenberg–Richter moment distribution is fitted by mens of a least-squares inversion following Frohlich (2007).

3.2.2 Temporal distribution of aftershocks

The frequency distribution of the decrement of earthquake aftershocks is described by the modified Omori's law (Utsu, 1961; Utsu et al., 1995) as

$$\begin{array}{}\text{(2)}& R\left(t\right)={\displaystyle \frac{k}{{\left(\mathrm{t}+\mathrm{c}\right)}^p}},\end{array}$$

where R(t) is the rate of occurrence of aftershocks within a given magnitude range, t is the time interval from the mainshock, k is the productivity of the aftershock sequence, p is the power-law exponent (p value), and c is the time delay before the onset of the power-law aftershock decay rate. Variations in p values exist for different tectonic regimes and each aftershock sequence. Many authors have related the p value with crustal temperature, heat-flow, or rock heterogeneity in the fault zone. Thus, relevant information can be extracted from these aftershock parameters in order to have a better understating of the rupture process of oceanic earthquakes. As before, we used the ZMAP software package (Wiemer, 2001) for estimating the p value of the aftershock sequence of the 1 May 1997 earthquake (M_w=6.9).

3.2.3 Fragment-asperity model

Alternative statistical models that relate the earthquake magnitude distribution with the rheology of the fault have been proposed. Among them, we have the fragment-asperity model. This model was introduced by Sotolongo-Costa and Posadas (2004) to describe the earthquake dynamics in a Tsallis entropy nonextensive framework (Tsallis, 1988). This model takes into consideration the irregular surfaces of two fault planes in contact and the rock fragments of different shape and sizes that fill the space between them. According to this model, earthquakes are triggered by the interaction along the fault planes of these rock fragments. Considering that large fragments are more difficult to release than small ones, the resulting energy is assumed to be proportional to the volume of the fragment (Telesca, 2010). Silva et al. (2006) improved the model and found a scaling law between the released energy (ε) and the size of asperity fragments (r) by the following proportional factor: ε∝r³. The nonextensive statistics are used to describe the volumetric distribution function of the fragments. A parameter that represents the proportion between ε and r is introduced. This parameter is known as the a value or parameter a (Silva et al., 2006; Telesca, 2010). The parameter a is defined using a volumetric distribution function of the fragments applying the maximum entropy principle for the Tsallis entropy (for details in the mathematical expressions see Silva et al., 2006; Telesca, 2010). The magnitude cumulative distribution function becomes

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\log }_{\mathrm{10}}\left(N>M\right)& ={\log }_{\mathrm{10}}\left(N\right)+\left({\displaystyle \frac{\mathrm{2}-q}{\mathrm{1}-q}}\right){\log }_{\mathrm{10}}\\ & \left[\mathrm{1}-\left({\displaystyle \frac{\mathrm{1}-q}{\mathrm{2}-q}}\right)\left({\displaystyle \frac{{\mathrm{10}}^K}{a^{\mathrm{2}/\mathrm{3}}}}\right)\right],\end{array}\end{array}$$

where N is the total number of earthquakes; N (>M) represents the number of events with magnitude larger than M; a is a proportionality parameter between ε and r, and q is the nonextensivity parameter. K is defined as K=2M (Silva et al., 2006) or K=M (Telesca, 2011). The magnitude (M) is related to ε by the following relation: $M=\mathrm{1}/\mathrm{3}$ log(ε) (Silva et al., 2006). Telesca (2011) considered that the relation between ε and M is given by $M=\mathrm{2}/\mathrm{3}$ log(ε). Neither of the two models is preferred over the other. We used both models in order to quantify the variability in the nonextensive parameters. According to Telesca (2010), the physical meaning of the q parameter consists in that it provides information about the scale of interactions. It means that if q is close to 1, the physical state is close to the equilibrium. As a result, few earthquakes are expected. On the other hand, as q rises, the physical state goes away from the equilibrium state, this implies that the fault planes are able to generate more earthquakes, thus resulting in an increment in the seismic activity (Telesca, 2009, 2011). The physical meaning of the a value lies in the fact that it provides a measure of the energy density. It means that the a value is large if the energy released is large (Telesca, 2011). For example, high a values are expected when the events with the highest magnitude take place. Previous studies have shown that the q value ranges mainly from 1.50 to 1.70 (Vilar et al., 2007; Vallianatos, 2009; Rodríguez-Pérez and Zúñiga, 2017, among others). We obtained the a and q parameters by minimizing the root mean square error (RMS) with the Nelder–Mead method (Nelder and Mead, 1965).

3.2.4 Stress inversion

Focal mechanisms are reliable indicators of the state of stress in a tectonic region. In order to study the regional stress field for oceanic earthquakes, we performed stress tensor inversion from focal mechanisms reported in the Global CMT catalogue (Dziewonski et al., 1981; Ekström et al., 2012) with the iterative joint inversion developed by Vavryčuk (2014). From the stress inversion, we obtained the orientation of the principal stress axes σ₁, σ₂, and σ₃ (where ${\mathit{\sigma }}_{\mathrm{1}}\ge {\mathit{\sigma }}_{\mathrm{2}}\ge {\mathit{\sigma }}_{\mathrm{3}}$) and the stress ratio R. We now briefly explain each method. The first method (the iterative joint inversion) provides an accurate estimation of R and stress orientations (Vavryčuk, 2014). In this method, the ratio is defined as $R=({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}})/({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{3}})$ (Gephart and Forsyth, 1984). A fault instability constraint is applied, and the fault is identified with that nodal plane, which is more unstable and thus more susceptible to faulting (Vavryčuk, 2014). By incorporating a fault instability constraint into the inversion, an iterative procedure is imposed. The uncertainties are determined as the differences between the inverted results considering noisy data (Vavryčuk, 2014). The stress inversion was carried out with the STRESSINVERSE software developed by Vavryčuk (2014). The maximum horizontal stress (SH_max) was calculated using the formulation of Lund and Townend (2007). The stress inversion was performed for each of the six different regions shown in Fig. 7.

Figure 3 Main statistical characteristics for intraplate oceanic events (INT). (a) Magnitude earthquake histogram; (b) frequency–magnitude distributions with M_c and b values. (c, d) The normalized cumulative number of events as function of magnitude for intraplate oceanic events (INT). Color curves show the best fit for the nonextensivity parameters q and a for the Telesca (2010) (red lines) and the Silva (2006) (green lines) models, respectively.

Figure 4 (a) Main statistical characteristics for the transform fault zone and mid-ocean ridge (TF-MOR) events (regions R1 to R5). (b) Magnitude earthquake histograms and frequency–magnitude distributions with M_c and b values for each of the different subregions shown in Fig. 8.

4 Results

There is a large span of b values (Table 2) which nevertheless sheds light on the seismicity characteristics of oceanic earthquakes in Mexico. INT events exhibit higher b values and M_c than TF-MOR events (Figs. 3, 4a and Table 2). In particular, TF-MOR events also show local b-value variations in the range of 0.72–1.30 (Fig. 4b) for each of the subregions R1 to R5 (Table 2). Previous studies had shown large fluctuations in b values of oceanic events. For example, Tolstoy et al. (2001) reported b values of about 1.5 associated with volcanic activity in the Gakkel Ridge. Läderach (2011) reported b values of 1.28 in the Southwest Indian Ridge. In a global study, Molchan et al. (1997) estimated the b value for mid-ocean and transform zones, obtaining values of the following interval 0.97–1.47. In general, our b-value estimates agree with reported b values in previous studies. On the other hand, our results showed that M_c for oceanic events is higher than reported M_c for the subduction zone and continental regions of Mexico, which reflects the capability of the global and regional networks to appropriately register events in that region. The magnitude completeness for oceanic earthquakes differs for different parts of the world, but in most cases, it is in the range of 4.0–5.0 on average considering most of the global catalogues.

Table 2 Statistical parameters.

INT are intraplate oceanic events; TF-MOR are transform fault zone and mid-ocean ridge events; M_c is the completeness magnitude; b is the slope of the Gutenberg–Richter distribution; q_S and a_S and q_T and a_T are the nonextensive parameters based on Silva et al. (2006) and Telesca (2011), respectively.

Our results also showed that transform fault zone and mid-ocean ridge events follow a tapered Gutenberg–Richter distribution, as suggested in previous studies (Boettcher and McGuire, 2009). The tapered Gutenberg–Richter distribution was fitted with the following parameters: β=0.64 and an estimated corner magnitude of M_m=6.7 (Fig. 5a). These results are in agreement with previous studies such as that of Bird et al. (2002), which studied the tapered Gutenberg–Richter distribution for spreading ridges and oceanic transform faults based on global data, obtaining a β value of about 0.67 for both types of events. They reported that M_m varies from 5.8 to 6.6–7.1 for mid-ocean ridge and transform faults, respectively. The results for the nonextensive parameters are shown in Table 2. We found higher q values for TF-MOR events than for INT events (Fig. 5), meaning that TF-MOR events are farther from the equilibrium than INT events. The results showed a better fit for cumulative distribution functions using the Telesca model for TF-MOR and each of the regions (Fig. 6). In regions R1-R5, our results showed that q varies from 1.31 to 1.52 and from 1.57 to 1.63 using Telesca's and Silva's models, respectively. In the case of subduction zones, the q value can vary from 1.35 to 1.70. For example, in the Hellenic Subduction Zone, q is in the range of 1.35–1.55 (Papadakis et al., 2013); in the Mexican Subduction Zone, Valverde-Esparza et al. (2012) found that q varies from 1.63 to 1.70. Thus, our results conform to values obtained in regional studies.

Figure 5 (a, b) The cumulative annual seismic moment frequency distribution for the transform fault zone and mid-ocean ridge (TF-MOR) events (regions R1 to R5). The blue lines are the moment, tapered Gutenberg–Richter distributions. The red lines represent the ordinary moment Gutenberg–Richter distributions. (c, d) The subregions that do not follow an ordinary moment Gutenberg–Richter distribution are subregions R1 and R2.

Figure 6 The normalized cumulative number of events as function of magnitude for the transform fault zone and mid-ocean ridge (TF-MOR) events. Blue triangles show the completeness magnitude (M_c). Red curves show the best fit for the nonextensivity parameters q and a for Telesca's (2010) model (red lines). Green curves show the best fit for the nonextensivity parameters q and a for Silva's (2010) model (green lines).

The analysis of the aftershock sequence of the 1 May 1997 earthquake (M_w=6.9), yielded a p value of 0.67±0.33 (Table 3). The magnitude of the largest aftershock of the 1997 event was M_w=5.3 (Table 3). Oceanic strike-slip events seem to have lower p values than mid-ocean ridge events. For example, Bohnenstiehl et al. (2004) found a p value of 0.95 for the 15 July 2003 (M_w=7.6) central Indian Ridge strike-slip event. For the Siqueiros, Discovery, and western Blanco transforms, the p value varies from 0.94 to 1.29 (Bohnenstiehl et al., 2002). Davis and Frohlich (1991) determined a p value of 0.928±0.024 for the combined ridge and transform environments. Our results fall within the range of global studies that showed that the p value varies from 0.6 to 2.5 (Utsu et al., 1995). We also reported a c close to 0 for the aftershock sequence of the 1 May 1997 (M_w=6.9) (Table 3). Shcherbakov et al. (2004) found that the parameter c of the Omori's law decreases as the magnitude of events considered increases. According to the study, this observation is due to the effect of an undercount of small aftershocks in short time periods. This provides an explanation for our result of c∼0 because of the limited magnitude detection reported in the regional and global catalogues used.

Table 3 Aftershocks characteristics of 1 May 1997 event.

M_m is the magnitude of the mainshock; M_a is the magnitude of the largest aftershock; D is the difference in magnitudes of the mainshock and its largest aftershock; p, c, and k are the coefficients of Omori's law.

We classified the focal mechanisms used in the stress inversion into seven categories: (1) reverse, R; (2) reverse with lateral component, R-SS; (3) strike-slip with reverse component, SS-R; (4) strike-slip, SS; (5) strike-slip with normal component, SS-N; (6) normal with lateral component, N-SS; and (7) normal, N (Fig. 7). This classification was performed to identify the dominant type of faulting for each subregion. Region R1 is composed of strike-slip (70.3 %), strike-slip with normal and reverse components (21.6 % and 5.4 %, respectively), and normal-faulting (2.7 %) focal mechanisms (Fig. 7b). Region R2 exhibits the following focal mechanism distribution: strike-slip (82.4 %) and strike-slip with normal and reverse components (9.5 % and 8.1 %, respectively) (Fig. 7b). In region R3, the focal mechanism classification shows the following distribution: strike-slip (62.5 %), strike-slip with normal component (25 %), normal-faulting with strike-slip component (6.3 %), and reverse (6.3 %) (Fig. 7b). Region R3 consists of strike-slip (70.8 %), strike-slip with normal and reverse components (8.3 % and 16.7 %, respectively), and reverse earthquakes (4.2 %) (Fig. 7b). Region R5 exhibits the following focal mechanism distribution: strike-slip (53 %), strike-slip with normal and reverse components (23.5 % and 17.6 %, respectively), and reverse (5.9 %). For the case of earthquakes in R6, the classification shows the following distribution: normal (83.3 %) and normal-faulting with strike-slip component (16.7 %) (Fig. 7b).

Table 4 Stress inversion results.

Stress ratio is defined by $R=({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{2}})/({\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{3}})$; ^* stress inversion based on Vavryčuk (2014) and Lund and Townend (2007). Location of the regions are shown in Fig. 1.

Table 4 summarizes the results from the stress inversion. Based on the orientation of stress axes, a dynamical description of the tectonics of the oceanic earthquakes in Mexico can be carried out. A quantitative comparison with other oceanic regions is discussed in what follows. Region R6 is only dominated by N and N-SS earthquakes (Fig. 8). In regions R4 and R5, stress results showed moderate similarities. The differences in these regions may also be related to the variability in focal mechanisms (here we have SS, SS-N, SS-R, and to a lesser extent R events) (Fig. 8). Variations are very significant in regions R1 to R3 (particularly in σ₂) (Table 4). These regions also showed different types of events: SS, SS-N, and SS-R for R1; SS, SS-N, and SS-R for R2; and SS, SS-N, N-SS, and R for R3 (Fig. 8). In these regions, strike-slip earthquakes are the dominant type of faulting. Events with unusual mechanisms have also been reported in other oceanic regions. According to Wolfe et al. (1993), most of the anomalous seismic activity is associated with mislocations, complex fault geometry, or large structural features with an influence on the slip of the fault. DeMets and Stein (1990) showed that the strike direction and earthquake slip vectors in the Rivera Transform are rotated clockwise from the expected direction of the Pacific–Rivera Euler vector. This deviation can be the result of morphologic features resulting in unusual patterns of epicenters and focal mechanisms.

Figure 7 Focal mechanism solutions of oceanic earthquakes in Mexico reported by the Global CMT catalogue from 1976 to 2017. (a) Focal mechanisms are divided into six regions (R1 to R6) for the stress inversion analysis. (b) Focal mechanism classification based on the Kaverina et al. (1996) projection technique implemented by Álvarez-Gómez (2015): reverse, reverse with lateral component, strike-slip with reverse component, strike-slip, strike-slip with normal component, normal with lateral component, and normal (R, R-SS, SS-R, SS, SS-N, N-SS, and N, respectively).

In the case of the East Pacific Rise Rivera segment (region R1), σ₂ is almost vertical, and SH_max is ∼170^∘, suggesting a strike-slip regime (Table 4). The main orientations of the P axes are in the N–S, NW–SE, and E–W directions. The orientations of the P axes are NW–SE and, to a lesser extent, E–W directions (Fig. 9). For the case of the Rivera Transform (region R2), σ₂ is quasi vertical, and the SH_max is 157^∘, suggesting a strike-slip regime. The orientation of the P axes is in the NW–SE direction, and the orientation is in the NE–SW direction for the T axis. In region R3, σ₂ is almost vertical, and the SH_max is also 157^∘, suggesting a strike-slip regime. The orientation of the P axis is in the NW–SE direction. The main orientation of the T axes is NE–SW, but E–W directions occur as well. For the region R4, σ₂ is 76, and the SH_max is 22^∘, suggesting a strike-slip regime. The predominant orientations of the P and T axes are NE–SW and NW–SE, respectively. In R5, σ₂ is from 69^∘, and the SH_max is 120^∘, suggesting a strike-slip regime. The main orientation of the P axes is NW–SE, while that of the T axis is NE–SW. In R6, the principal axes are related to a normal fault regime. σ₁ is almost vertical, and the SH_max is ∼45^∘. The orientation of the T axes is in the NW–SE direction. The Mohr circle diagram showed that most of the studied events are clustered along the outer Mohr's circle in the area of validity of the Mohr–Coulomb failure criterion (Fig. 9). Reported focal mechanisms confirm Sykes's model for mid-ocean ridges (Sykes, 1967), where events in transform zones tend to have strike-slip mechanisms, while ridge crest events have mainly normal faults. The obtained orientation of the principal axes supports this model.

Figure 8 Orientation of horizontal axes. (a) Maximum horizontal stresses (SH); (b) Minimum horizontal stresses (Sh).

5 Discussion

One of the main problems for studying oceanic seismicity is that the epicenters are located far from most of the recording stations in mainland Mexico. This has a direct effect on the earthquake magnitude distributions (M_c and b value). We first discuss the magnitude completeness of oceanic earthquakes. Global studies showed that the magnitude completeness for oceanic earthquakes is in the range of 4.0–5.0. Our results are in agreement with these global studies. However, as expected, several microseismic surveys which have been conducted in different oceanic environments (e.g., Smith et al., 2003; Simão et al., 2010; McGuire et al., 2012, among others) can yield lower-magnitude thresholds. As a result of these studies, precise hypocenter locations and earthquake distributions with a broader magnitude range were obtained. Thus, lower M_c has been reported for studies based on microseismic surveys. For example, in the Mid-Atlantic Ridge, M_c∼3.0 with several smaller events (M_w<2.5) were reported (Bohnenstiehl et al., 2002; Smith et al., 2002, 2003).

Figure 9 Mohr circle diagrams for all the regions (left column). P and T axes distributions for all the regions (right column). Red circles represent pressure, while blue crosses represent tension.

Another factor that has to be discussed is the accuracy in the location of the epicenters. The location uncertainty plays an important role when earthquakes are assigned to an intraplate or a mid-ocean ridge/transform fault environment. For example, some studies reported that for faults located at 4S on the EPR, teleseismic locations could be off by as much as 50 km (McGuire, 2008; Wolfson-Schwehr, 2014). As a consequence, some TF-MOR events are probably classified as INT events and vice-versa (for example, epicenters in color in Fig. 2). Some events located in the Tamayo Fracture Zone close to the Rivera Subduction Zone may also be misidentified. This mislocation effect introduces uncertainties in the estimation of the statistical parameters useful for understanding the tectonics of the region. In order to have precise locations and avoid mislocation, ocean-bottom seismometers off the Mexican coast would be needed. Being aware of this, one should avoid over-interpretation of the results. Local monitoring of oceanic events represents an improvement of more than an order of magnitude relative to the regional and teleseismic detection levels.

Previous studies also showed that the seismicity near oceanic transform faults that connect mid-ocean ridges may be thermally controlled (Abercrombie and Ekström, 2001; Boettcher et al., 2007). The thermal effect is most evident in the seismogenic zone. It is essential to mention that faults along the middle and southern segments of the EPR are shorter and faster slipping. The faster slip rates and shorter fault lengths result in narrower seismogenic zones because the thermal structure is shallow. On the other hand, the Rivera Transform is longer and has a slower slip rate, resulting in a wider seismogenic zone. However, heat is not the only factor that regulates seismicity, because the largest events break a small part of the rupture areas predicted by thermal models (Boettcher and Jordan, 2004; Roland et al., 2010). Thus, most slip occurs without producing large earthquakes (Goslin, 1999; Boettcher and Jordan, 2004; Roland et al., 2010). This can explain the occurrence of a few events with M>6.5 in the Rivera Transform. According to McGuire et al. (2012), the apparent lack of large events on mid-ocean ridge transform faults may also be related to the heterogeneity of materials on the fault plane. The maximum magnitude for transform fault events on the East Pacific Rise (in the latitude interval of 3^∘ < Lat < 5^∘) is about 6.5 (McGuire et al., 2005). On the other hand, earthquakes in the Rivera Transform and on the northern segment of the East Pacific Rise (in Mexico) have relative larger magnitudes (M>6.8) based on reported seismicity in different catalogues (Fig. 1). This highlights a differentiation between the middle and southern and northern segments of the East Pacific Rise.

A further aspect of the analysis of oceanic earthquakes is their capacity to generate aftershocks as well as the their occurrence, duration, and magnitude. Earthquake statistical studies showed that large oceanic events in transform faults, fracture zones, and intraplate regions release low energy levels in their aftershock sequences (Houston et al., 1993; Boettcher and Jordan, 2001; Antolik et al., 2006). Boettcher et al. (2012) found that earthquakes on transform faults have an order of magnitude fewer aftershocks than intraplate events. According to some authors, a low aftershock-to-mainshock energy ratio indicates an efficient rupture or complete stress drop in the mainshock presupposing a weak fault (Hwang and Kanamori, 1992; Velasco et al., 2000). Many factors can affect the aftershock productivity. For example, the age of the lithosphere and the heat flux have a direct influence on the rock strength (Antolik et al., 2006), thus, explaining the low energy release in the aftershock sequence of oceanic events. The observed low aftershock energy seems to be a common feature of oceanic earthquakes (Antolik et al., 2006). In this regard, we studied the 1 May 1997 (M_w=6.9) strike-slip event in the Rivera Transform and its largest aftershock (M_w=5.3). By considering the energy magnitude as log $E=\mathrm{1.5}{M}_{\mathrm{w}}+\mathrm{11.8}$, we obtain that the energy of the mainshock is 1.41×10¹⁵ Nm, and the energy of the largest aftershock is 5.62×10¹⁵ Nm, resulting in an aftershock-to-mainshock energy ratio of 0.003. This value is considered as low and representative of strike-slip events, as shown by the comparison with the results reported by Velasco et al. (2000).

A similar analysis comes from Båth's law by considering the magnitude difference between the mainshock and the largest aftershock. We determined that the magnitude difference for the 1997 event is 1.6, which is higher than the theoretical value of 1.2. Both magnitude difference and the aftershock-to-mainshock energy ratio showed large scatter (e.g., Velasco et al., 2000; Utsu, 2002), and results ought to be taken with caution. The aftershock decay rate is the product of the strain relaxation around the rupture plane. Aftershock studies have shown that oceanic ridges are prone to having larger p values than those of subduction zone regimes due to the high temperature of the oceanic crust, which results in rapid strain release (Kisslinger, 1996; Rabinowitz and Steinberg, 1998; Klein et al., 2006). According to previous studies, extremely high p values (p>2) and short aftershock durations are related to high temperatures (Bohnenstiehl et al., 2002; Simão et al., 2010) and/or migration of hydrothermal fluids (Goslin et al., 2005). We found a p value of 0.67±0.33 for the 1 May 1997 (M_w=6.9) strike-slip event in the Rivera Transform. This p value is consistent with other oceanic regions, but it does not seem to conform to a high-temperature regime.

Regarding the magnitude distribution of oceanic events, our b-value estimates are in agreement with global oceanic studies but differ from local studies. For example, along the East Pacific Rise (in the latitude interval of 5^∘ N < Lat < 9.90^∘ N), b-value estimations fluctuate from 1.10 to 2.50 (Bohnenstiehl et al., 2008). Bohnenstiehl et al. (2008) determined the b value of 9000 microearthquakes with magnitudes in the range of −1.5–1.0 located in the southern part of our study zone. Due to this overlap, we compare their results with our results for region R5. For this region, we obtained a b value of 0.94 with a M_c of 4.2, while they found that the b value approaches 2.5 at very shallow depths (< 0.3 km) (with $M_{\mathrm{c}}=-\mathrm{1.3}$). At depths of 0.5–1.5 km, the b values drop to a value of 1.10 (with $M_{\mathrm{c}}=-\mathrm{0.4}$). According to Bohnenstiehl et al. (2008) at very shallow depths, the uppermost oceanic crust is structurally heterogeneous because of the extrusion of lava and the repeated emplacement of sheeted dikes. As a consequence, there is a large proportion of small versus large earthquakes resulting in high b values. The b values decrease with depth due to the decreasing heterogeneity and/or changes in ambient stress levels. Considering that events in our catalogue for R5 occur at a different depth interval, and assuming the decreasing heterogeneity, fewer magnitude events would be expected (reducing the b value). Another explanation for the differences between our results and the results of Bohnenstiehl et al. (2008) is that the magnitude ranges of the earthquake catalogues are extremely different. This highlights how the b value is affected by magnitude completeness.

Statistical studies suggested that the β value mainly takes values between 0.60 and 0.70 for a global range (Kagan, 2002). Our estimates of β agree with global oceanic studies. It is essential to discuss the tectonic implications of this parameter. Bird et al. (2002) also found a dependence of β value on the relative plate velocity. According to them, the β value is higher (with M_m=7.1) when the velocity is < 36 mm yr⁻¹ than when the velocity is > 67 mm yr⁻¹ (with M_m=6.6) for spreading ridges and oceanic transform faults, respectively. These observations are in agreement with our estimate of β=0.64 and M_m of 6.6 for oceanic earthquakes in Mexico (Fig. 5). For intraplate events, we obtained a β>0.70. According to Kagan (2010), β values > 0.70 may be related to the mix of earthquake populations with different maximum magnitudes (M_m). In the case of intraplate events, we associated the somewhat high β values with the mix of some intraplate and mid-ocean–transform events. This could be related to incorrect hypocenter locations due to the difficulty of precisely locating oceanic events by the land-based networks.

The seismicity models based on nonextensivity consider the interaction of two irregular fault surfaces (asperities) and rock fragments filling them. However, these models differ in their assumption of how energy is stored in the fragments and the asperities. This difference is expressed through the constant a, which represents the proportionality between the released energy E and the fragment size r. This explains the difference in a parameter between Telesca's and Silva's models (Fig. 5). Both models showed that a for TF-MOR is higher than a for the INT events (Fig. 5). This implies that more energy is released for TF-MOR earthquakes. On the other hand, the q value indicates if the physical state of a seismic area moves away from equilibrium. The physical state is at equilibrium when q is equal to 1, and as q increases, the system is in an instability state in which a more significant amount of seismic energy is released.

Finally, we discuss the focal mechanisms and the calculated state of stress for oceanic earthquakes in Mexico. Focal mechanisms provide useful information about the structure and settings of faults and can describe the crustal stress field in which earthquakes take place. Our analysis is limited because we only used focal mechanisms based on teleseismic data. The teleseismic detection threshold for oceanic events in the East Pacific Rise is dependent on the region of the EPR. For example, Riedesel et al. (1982) report a magnitude detection threshold in the range of 4.0–5.0. For the Quebrada, Discovery, and Gofar faults, the CMT catalogue is only complete to M_W=5.4. (McGuire, 2008; Wolfson-Schwehr et al., 2014). Another limitation of our study is that we combine different types of earthquakes into a single region, resulting in inaccurate estimations of the stress state for that specific region. Under these circumstances, our study provides information on the stress field of major structures or the stress associated with the dominant types of earthquakes.

In oceanic environments, the largest magnitude events along transform fault or intraplate earthquakes usually show strike-slip mechanisms (Wiens and Stein, 1984; Kawasaki et al., 1985). In the areas adjacent to the oceanic ridges where the oceanic lithosphere is young, Wiens and Stein (1984) reported a large variety of focal mechanisms and stress orientations. For example, in the East Pacific Rise, in the Mexican territory, Wiens and Stein (1984) reported thrust and normal mechanism solutions for near ridge intraplate seismicity. This explains the strike-slip with normal components, as well as thrust events in regions R3, R4, and R5 (Fig. 7). In R3 and R4 (Fig. 7), the maximum horizontal axes (compression) of thrust events show a preferred orientation perpendicular to the spreading direction. On the other hand, in region R5 (Fig. 7), the compression axes, showed a weak preferred alignment with respect to the spreading direction. In the Rivera Transform, focal mechanisms showed right lateral strike-slip motion implying oblique horizontal stresses (Fig. 7). Although most of the events in the Rivera Transform (R2 in Fig. 7) are strike-slip events, some events with unusual mechanisms have been reported (normal faulting events) (Wolfe et al., 1993). Normal faulting events may be related to extensional offsets or internal deformation of the Rivera Plate (Wolfe et al., 1993).

6 Conclusions

We analyzed the seismicity of oceanic events in the Pacific oceanic regime of Mexico. Oceanic earthquakes were classified into two different categories: intraplate oceanic (INT) and transform fault zone and mid-ocean ridge (TF-MOR) events, respectively. We conducted a stress state estimation for the different regions. Because of the combination of different types of earthquakes into the regions, our results only provide information on the stress field of major structures or the stress associated with the dominant types of earthquakes. It is important to be aware of this limitation in order to avoid an over-interpretation of the results. TF-MOR events have strike-slip, strike-slip with normal and reverse components, normal and normal-faulting with the strike-slip component, and reverse focal mechanisms. On the other hand, INT events have only normal and normal-faulting with strike-slip component focal mechanisms. The stress field from INT and TF-MOR events agree with global studies. Regarding the aftershock productivity, we found that the aftershock decay rate of the 1 May 1997 (M_w=6.9) strike-slip event in the Rivera Transform is also consistent with oceanic p-value estimations. Despite the limitation of the catalogues used, our results provided a comprehensive insight into the seismicity of oceanic environments. The main problem is the location uncertainty and mislabeling of the earthquakes. The b value for INT events (1.17) is higher than that for TF-MOR events (0.82). Our b-value estimations are in agreement with other regional studies but differ from b-value estimates based on microseismicity studies. Our b-value estimates for mid-ocean ridge/transform fault environments are lower (0.72<b<1.30) than those derived from microseismicity studies (1.1<b<2.5). Our results also showed that TF-MOR events mostly follow a tapered Gutenberg–Richter distribution.

From the nonextensivity analysis, we observed that TF-MOR events are farther from the equilibrium than INT events. Thus high q values take place in mid-ocean ridges and transform fault zones. This means that mid-ocean ridge and transform faults are able to produce more seismicity. Low q values are also reported during relatively quiet periods, characterized mainly by the occurrence of small-magnitude events. This can be an explanation for the low q values of regions R1 and R5. Our results also showed that a values are higher for TF-MOR events than for INT events using both models. This implies that more earthquakes with larger magnitude occur (or more energy is released) in mid-ocean ridge/transform fault environments than in an oceanic continental environment. Telesca's model fits better with the cumulative magnitude distribution functions making a better option to study the oceanic seismicity in Mexico.