Perspective Chapter: Testing the Interoccurrence Times Probability Distributions of Earthquakes

Abstract

This chapter aims to provide a methodology to test the probability distributions of earthquakes in terms of the interoccurrence times (ITs), namely, the time between consecutive earthquakes of a specific magnitude. First, we compile a new earthquake catalog for the El Salvador subduction zone within moment magnitude M 5.0–8.12 comprising historical and instrumental data for 1609–2019. Secondly, we explain the fundamentals of the Weibull and Poisson distributions and verify the IT probability fits when considering the clustered catalog. We find that the Weibull distribution fits all ITs, while the Poisson distribution fails to explain the natural seismicity patterns for small magnitude bins. Besides, we test the assumption that the declustering process leads to a Poisson probability distribution when removing foreshocks and aftershocks in the earthquake catalog. Finally, the classical Gutenberg–Richter relationship and conditional magnitude probabilities are calculated as an essential input in any seismic hazard assessment.

1. Introduction

1.1 Earthquake data

The subduction zone comprises earthquakes with focal depths down to 300 km in the Central America Isthmus, where the Cocos Plate submerges beneath the Caribbean Plate. Interface subduction earthquakes occur at shallower depths in the two plates’ boundary near the trench, characterized by thrust focal mechanisms. Intraplate subduction earthquakes are deeper and characterized mainly by normal focal mechanisms. Outer rise earthquakes occur southwest of the trench due to the initial flexure of the oceanic slab, yielding a shallow normal mechanism. The slab dips an average of 45⁰ north-east, with deeper earthquakes away from the subduction trench having offshore and inland epicenters. Correa-Mora et al. [1] suggested employing GPS measurements as a weak coupling between the Cocos and the Caribbean plates with a slip rate of about 7–8 cm/year between the two plates.

The last destructive earthquake in this seismogenic zone occurred on January 13, 2001, with an M 7.7 with its epicenter offshore El Salvador [2]. Large shocks occurred on December 19, 1862 [3] and September 7, 1915 [4] beneath the western El Salvador territory with M 8.12 and M 7.8, respectively. Ambraseys and Adams [4] reported an event dated August 6, 1942, M 8.12 beneath the western Guatemalan territory. On September 2, 1992, an M 7.6 offshore Nicaragua triggered a local tsunami [5]. On July 29, 1773, an earthquake with M 7.6 struck Guatemala [3].

We updated the homogenized moment magnitude catalog of Salazar et al. [6] using online global information as the International Seismological Centre ISC [7], the Preliminary Determination of Epicenters PDE [8], the Centroid Moment Tensor solutions [9], and detail-oriented studies for earthquake locations of large earthquakes [3, 4, 10, 11, 12]. The catalog covers the geographical window between 11 and 16.5⁰ N and 85.5–92⁰. The final subduction catalog contains 2764 events covering 1609–2019 within M 5–8.12 (Figure 1a). Note that earthquakes from other regional seismogenic sources (e.g., volcanic chain, Guatemalan faults, and Honduran–Guatemalan Grabens) were also compiled; however, this study is devoted to only subduction events. Salazar et al. [6] present the criteria to separate upper-crustal and subduction events based on geological and focal mechanism information. We also applied the Gardner and Knopoff [14] method to extract the main events from the original catalog, eliminating the foreshocks and aftershocks. The final declustered catalog contains 889 main events. Note that the method of Gardner and Knopoff [14] calibrated their method to identify aftershocks; however, it has also been applied to eliminate foreshocks with the same temporal and spatial windows of the former (e.g., [15, 16, 17]) and to assure all dependent events are removed from the catalog. About 68% of dependents’ events were eliminated in the decluster process (Figure 1b). Figure 2 shows a hypocenter cross-section for El Salvador and the USGS subduction slab 2.0 model [13]. The depth determination in the Central America Isthmus gives significant uncertainty; however, the USGS slab model generally depicts a reasonable top slab geometry when compared with our hypocentral determinations.

Figure 1.

Earthquake epicenters (circles) for the subduction zone comprising events with moment magnitude M 5.0–8.12, 1609–2019, and focal depth from 1 to 300 km. a) Clustered catalog (all events); the rectangle depicts the area of the cross-section AB-CD in Figure 2 b) Declustered catalog (after removing foreshocks and aftershocks). Blue contour lines depict the top Cocos plate depths every 20 km based on the USGS slab 2.0 model [13].

Figure 2.

El Salvador cross-section illustrates hypocenters for all subduction earthquakes with moment magnitude M 5.0–8.12. Dashed lines depict the top of the Cocos plate based on the USGS slab 2.0 model [13] for the area AB-CD in Figure 1. Earthquakes with a fixed depth of 33 km are not plotted. The red triangle indicates the position of the volcanic chain in Central America, and the letter T indicates the position of the trench between the Cocos and the Caribbean plates.

It is necessary to investigate the year of completeness for the magnitudes, especially because small events are not listed in the catalog for early years. We present in Table 1 the year of completeness of the subduction catalog for several magnitude bins after applying the Tinti and Mulargia [18] method for both the clustered and the declustered catalog after applying the Gardner and Knopoff [14] method; Figure 3 shows an example of completeness analysis for the magnitude bin 5.0–5.5, suggesting that the catalog entirely lists events from the year 1975 to 2019.

M	Year
5–5.5	1975
5.5–6	1955
6.0–6.5	1920
6.5–7	1900
7–7.5	1800
7.5–8	1765
8–8.12	1862

Table 1.

Estimated year of completeness for the subduction zone. We noticed that for all magnitude bins, the year of completeness is the same when employing the clustered and the declustered catalogs (see Figure 3).

Figure 3.

Example of the completeness period for subduction earthquakes within a magnitude bin M 5.0–5.5 employing the a) clustered catalog: Main shocks, foreshocks, and aftershocks; b) declustered catalog: Only main shocks. The catalog is complete for this magnitude interval back to 1975 because the cumulative annual number of earthquakes for this magnitude range is approximately linear back to that date [18]. All other magnitude bins’ periods of completeness are presented in Table 1.

Salazar [16] studied the earthquake interoccurrence times and seismic hazards for earthquakes related to another seismogenic zone, namely, the upper-crustal volcanic chain in El Salvador. This work extends the former article by employing data from the overall subduction zone and studying the interoccurrence times based on the Weibull and Poisson cumulative probability distributions.

2. The Weibull and Poisson cumulative distributions

A straightforward way to investigate if a specific probability distribution fits the time seismicity patterns in a seismogenic zone is by applying the cumulative probability distribution based on the interoccurrence times (ITs), namely, the time between consecutive earthquakes of a specific magnitude bin. Salazar [16] studied the ITs matching the Weibull and Poisson cumulative distributions for the upper-crustal volcanic chain earthquakes in El Salvador. Other authors have also studied the ITs in other parts of the globe [19, 20, 21, 22, 23]. Note that we count the ITs in a magnitude bin (e.g., 5–5.5, 5.5–6.0) since the final objective is to use the probability of occurrence of events in a seismic hazard assessment.

The cumulative F(t) Weibull distribution yields [24]:

$$F\left(t\right)=1-e^{-{\mathit{\alpha t}}^{\beta }}t>0,\alpha ,\beta >0$$E1

where t is the earthquake ITs, and α and β are constants found by the non-linear search algorithm of Bean et al. [25]. The Poisson probability cumulative distribution F(t) yields:

$$F\left(t\right)=1-e^{-\mathit{\lambda t}}t>0,\lambda >0$$E2

where λ is the number of earthquakes per unit of time [26], and t must be calculated based on the completeness analysis explained in the previous section. Table 2 shows the α and β Weibull and the Poisson λ constants for the clustered and declustered earthquake catalogs and their corresponding mean and standard deviations for several magnitude bins. The mean μ for the Weibull distribution is:

	Weibull clustered catalog					Poisson clustered catalog		Poisson declustered catalog
M	α	β	RMS	Mean μ years	σ years	λ	RMS	λ	RMS
5–5.5	0.219	0.767	0.010	0.02	0.03	0.116	0.035	0.028	0.030
5.5–6	0.038	0.822	0.017	0.16	0.20	0.017	0.038	0.008	0.020
6–6.5	0.017	0.733	0.023	0.84	1.17	0.003	0.084	0.002	0.028
6.5–7	0.006	0.762	0.028	2.46	3.27	0.001	0.064	8 x 10^−4	0.050
7–7.5	6 x 10^−4	0.964	0.051	6.19	6.44	3.6 x 10^−4	0.079	9 x 10^−4	0.083
7.5–8.12	4 x 10^−4	0.863	0.085	24.91	28.97	1.0 x 10^−4	0.094	1 x 10^−4	0.094

Table 2.

Weibull and Poisson cumulative probability distribution parameters (Eqs. 1 and 2). σ denotes the standard deviation. M: moment magnitude; λ is expressed above as the number of earthquakes per year; however, for the Poisson cumulative probability calculations, the number of earthquakes per day must be used in Eq. (2).

RMS: root mean square; μ: mean

$$\mu ={\alpha }^{-1/\beta }\Gamma \left(1+\frac{1}{\beta }\right)$$E3

and the variance σ² gives:

$${\sigma }^2={\alpha }^{-2/\beta }\left\{\Gamma \left(1+\frac{2}{\beta }\right)-{\left[\Gamma \left(1+\frac{1}{\beta }\right)\right]}^2\right\}$$E4

The standard deviation σ is the root square of the variance. Γ is the gamma function; we also listed the fit error in terms of the root mean square (RMS):

$$\mathit{RMS}=\sqrt{\frac{{\sum }_{i=1}^N{\left(\mathit{\text{Prob}}\_{\mathit{Obs}}_i-\mathit{\text{Prob}}\_{\mathit{\text{Pred}}}_i\right)}^2}{N}}$$E5

where N is the number of data, Prob_Obs_i is the observed cumulative probability, and the Prob_Pred_i is the predicted cumulative probability at the i event.

We infer by comparing the observed and the predicted ITs using the obtained constants that subduction events pose the Weibull cumulative distribution when considering all the events in the subduction catalog, namely, principal, fore, and aftershocks (clustered catalog) for all magnitude bins with low RMS between 0.01 and 0.085 (Figure 4a,c,e,g,i, and k). Abaimov et al. [19] suggested that the conventional graphs depicting the cumulative probability F(t) vs. the ITs (left panel in Figure 4) are not appropriate for judging, in the first instance, the fit between observed and predicted probabilities. Instead, they propose linearizing the cumulative probability, applying -ln(1-F(t)). The linearized probability Weibull plots confirm the goodness to fit (Figure 4b,d,f,h,j, and l) in all cases. However, the Weibull distribution better predicts the probabilities for longer ITs than short ones in one case after 100 days (e.g., M 6.5–7). The α value yields practically zero for M 7.5–8.12 (Table 2), arguing a Poisson process tendency for the largest subduction shocks. Indeed, the Poisson and Weibull models predict similar probabilities (see the right panel in Figure 4j and l). Although the Poisson probability distribution is thought to be applied to independent events only, we also tested the ITs to the clustered catalog to investigate if such distribution fits some magnitude bins under consideration, especially for big events.

Figure 4.

Interoccurrence times vs. cumulative probabilities for several magnitude bins. (a) M 5–5.5 (c) M 5.5–6 (e) M 6–6.5 (g) M 6.5–7 (i) M 7–7.5 (j) M 7.5–8.12. We compare the observed and the estimated cumulative probability distributions employing the Weibull and Poisson models based on the clustered catalog. The arrow depicts Weibull’s mean. The right side (b, d, f, h, j, and k) depicts the linearized probability plots.

The Poisson cumulative distribution does not fit the ITs for smaller events with magnitudes between 5.5 and 7 when using the clustered catalog (all events in the analysis), yielding RMS from 0.035 to 0.064 (Table 2). The linearized plots on the right panel in Figure 4d,f, and h confirm such a statement for M 5.5–6, 6–6.5, and 6.5–7.0, respectively.

However, when applying the Gardner and Knopoff [14] method to remove the foreshocks and aftershocks and conform to the declustered catalog, the Poisson distribution better fits the ITs for M 5.5–6 and 6–6.5 (Figure 5d and f) but still fail to reproduce the ITs for M 6.5–7.0 (Figure 5h), although there are lower RMS yields for these cases (Table 2). The linearized probability plots confirm the goodness of the Poisson fit in most of the magnitude bins under analysis, including large shocks above M 7.5 (Figure 5j and l). We conclude that the Poisson distributions fit the time seismicity patterns once only when the main shocks are considered in the analysis, except for the magnitude bin between 6.5 and 7.0 when ITs are less than 100 days.

Figure 5.

Interoccurrence times vs. cumulative probabilities for several magnitude bins.

(a) M 5–5.5 (c) M 5.5–6 (e) M 6–6.5 (i) M 7–7.5 (j) M 7.5–8.12. We compare the observed and the estimated cumulative probability distributions employing the Poisson model based on the declustered catalog. The right side (b, d, f, h, j, and k) depicts the linearized probability plots.

3. Seismicity evaluation

Seismicity evaluation is a fundamental part of any seismic hazard assessment where the initial objective is to retrieve the magnitude probability of occurrence. In the case of a time-dependent seismic hazard, once we ensure that the Weibull cumulative distribution fits the subduction seismicity time patterns, we calculate the conditional probability CP(M)t, Δt yielding the probability that an earthquake occurs after an elapsed time Δt once an earthquake has happened at time t. Note that time t is the last year in which an earthquake of a specific magnitude M appears in the catalog [27]:

$$\mathit{CP}{\left(M\right)}_{t,\mathit{\Delta t}}=\frac{R\left(t\right)-R\left(t+\mathit{\Delta t}\right)}{R\left(t\right)}$$E6

where the reliability function is $R\left(t\right)=1-F\left(t\right)$ . Figure 6 shows the conditional probabilities for each magnitude bin under consideration until 2120. For M 5–5.5 and 5.5–6.0, the conditional probability is practically unity after 2019, while for the rest of the magnitudes, the probabilities increase at a lower rate as the size of the earthquakes increases and time passes. For example, for large shocks with M 7.5–8.12, the probability of occurrence after 2001 (the last destructive earthquake in the region happened that year) yields 0.28, 0.73, and 0.95 for the years 2025, 2050, and 2100, respectively. Such values are suitable for time-dependent seismic hazard assessment [16]. M 6–6.5 and 6.5–7 yield similar conditional probabilities as the elapsed time increases.

Figure 6.

Weibull conditional probabilities vs. elapsed times (Eq. 6). Note that the probabilities for M 5.5–6 completely overlap the ones from M 5–5.5 since the yield is unity.

We finally test our subduction catalog deriving the classical Gutenberg–Richter (G-R) relationship, yielding:

$$logN=A-\mathit{BM}$$E7

where N is the number of earthquakes per year with a magnitude equal to or above M. A and B are constants obtained by regression analysis. Note that the values of N must be taken after the year of completeness analysis (see Section 1 and Table 2) to avoid underestimation of seismicity levels in the magnitude bins. The relationships are obtained for the clustered and the declustered catalog (see Figure 7a). The G–R relationships yield log N = 6.43–1.01 M, σ = ±0.045 and log N = 7.73–1.17 M, σ = ±0.097 for the declustered and clustered catalogs, respectively. Note that we retrieved a B value of 1.0 with a very low standard deviation σ for the subduction earthquakes when using only the main events in the analysis, which is a characteristic value of tectonic earthquakes worldwide used in several seismic hazard analyses. Note that such G–R relationship must be truncated to the maximum possible magnitude when employing it in a seismic hazard assessment. Since we have yet to compute design spectra, we opt to present our results using the classical relationships based on Eq. (7).

Figure 7.

a) Classical Gutenberg–Richter relationships: The observed number of earthquakes per year N above a specific magnitude M (Eq. 7). The circles represent observed data; the lines represent N’s estimation of the G–R relationships after a regression analysis. The G–R yields log N = 6.43–1.01 M, σ = ±0.045 and log N = 7.73–1.17 M, σ = ±0.097 for the declustered and clustered catalogs, respectively; b) discrete probability function for 20 magnitudes between M 5 and 8.12 for the declustered catalog (Eq. 8).

For time-independent seismic hazard analysis, the probability that a magnitude will fall within a magnitude m1 and m2 is given by:

$$P\left(M\right)={\int }_{m1}^{m2}f\left(M\right)\mathit{dM}$$E8

where f(M) is the magnitude probability density function based on the constants A and B in Eq. (7) [28]. We present in Figure 7b the magnitude probability for the overall subduction zone suitable for time-independent seismic hazard analysis. The probability of the occurrence of small shocks is higher than the ones of large shocks. The sum of all probabilities is 1.0.

4. Conclusions

We have compiled a new earthquake catalog for historical and instrumental subduction earthquakes in El Salvador and surrounding areas that represent the input data for the ITs analysis and the seismicity evaluation. Subduction seismicity patterns fit the Weibull probability distribution well for all events (including main, foreshocks, and aftershocks). The Poisson probability distribution is not suitable for the natural seismicity for small events in the subduction zone. However, such distribution fits when dependents shocks are removed from the catalog for most of the magnitude range under scrutiny but fails to explain the interoccurrence times less than 100 days for the magnitude bin 6.5–7.0. Such failure might be explained because the Gardner & Knopoff [14] method was calibrated for Southern California crustal events rather than subduction zones. Indeed, each region should have its declustering dynamic window on time and space that belongs to the tectonics of the region under study. We conclude that the catalog declustering process is necessary when employing the classical Cornell seismic hazard assessments on time-independent schemes as previous works have applied for the region [6, 29, 30, 31]. In other words, despite the Poisson probability distribution not fitting the natural seismicity of subduction earthquakes, the declustering process produces a Poisson distribution with only independent events. However, the correspondent earthquake loads might be revised when employing the whole catalog in a seismic hazard analysis.

Since all the IT analyses belong to the overall subduction zone, we expect that future research must separate interface and intraplate shocks to consider the change in the seismic activity in depth; it is clear that below 120 km, there is a decrease in the seismicity in the Cocos Plate (see Figure 2). The depth determinations contain a significant error in the Central America Isthmus that can make such differentiation cumbersome, especially when there is no reported focal mechanism for some events.

A B value of 1.0 retrieved from the Gutenberg–Richter recurrence relationship after removing dependent events validates the quality of the subduction catalog compiled in this study. Conditional probabilities that depend on the last event listed in the catalog for a certain magnitude are also given that are useful in time-dependent seismic hazard assessments. The methodology provided in this chapter would improve the quality of the following seismic hazard assessment in the country.

Acknowledgments

The Catholic University of El Salvador financed this research under the Academic Researcher grants program devoted to full-time professors.