Estimating the Epicenter of a Future Strong Earthquake in Southern California, Mexico, and Central America by Means of Natural Time Analysis and Earthquake Nowcasting

It has recently been shown in the Eastern Mediterranean that by combining natural time analysis of seismicity with earthquake networks based on similar activity patterns and earthquake nowcasting, an estimate of the epicenter location of a future strong earthquake can be obtained. This is based on the construction of average earthquake potential score maps. Here, we propose a method of obtaining such estimates for a highly seismically active area that includes Southern California, Mexico and part of Central America, i.e., the area N1035W80120. The study includes 28 strong earthquakes of magnitude M ≥7.0 that occurred during the time period from 1989 to 2020. The results indicate that there is a strong correlation between the epicenter of a future strong earthquake and the average earthquake potential score maps. Moreover, the method is also applied to the very recent 7 September 2021 Guerrero, Mexico, M7 earthquake as well as to the 22 September 2021 Jiquilillo, Nicaragua, M6.5 earthquake with successful results. We also show that in 28 out of the 29 strong M ≥7.0 EQs studied, their epicenters lie close to an estimated zone covering only 8.5% of the total area.

Recently, our group has combined [19] NTA and EN together with the properties of the earthquake networks based on similar activity patterns (ENBOSAP) [20,21] for a similar purpose in the Eastern Mediterranean region where seismicity is much less intense than in the region N 35 10 W 120 80 of our present interest. We see that when strong EQs are frequent, as is the present case, we can simplify the approach of Reference [19] just by using NTA and EN.

EQ Data and the Tectonics of the Study Area
The area of interest here is Southern California, Mexico, and Central America, see Figure 1. We used the United States National Earthquake Information Center (NEIC) PDE catalog-these data are available from the United States Geological Survey (USGS), cf. [84]-in the region N 35 10 W 120 80 and considered all EQs recorded during the period of almost 50 years from 1 January 1973 to 21 September 2021 with magnitude M ≥ M thres = 4.0.
The interaction of the continental block with the oceanic provinces that surround Mexico has resulted in the current geographic layout. Mexico is located in a region where five tectonic plates: North America, Pacific, Rivera, Cocos, and Caribbean, are in regular interaction. The major fault zones, the spreading zones, and the subduction zones determine the boundaries between these plates [85]. The North America's displacement is to the southwest, the Eastern Pacific's to the northwest, the Cocos and Rivera's to the northeast, and the Caribbean's to the east; this disposition enhances the likelihood of a large seismic occurrence. Mostly, the interplay that occurs on the Mexican Pacific's southern coast, where the Cocos plate subducts beneath the North American plate, is the primary source of earthquakes in this country.
Oceanic plates are considered to be substantially more rigid than continental plates due to their olivine-rich composition. The Pacific Plate, the largest of the Earth's tectonic plates, is extensively an oceanic plate. The plate border deformation zones in the continental crust are, without a doubt, far larger (tens to hundreds of kilometers) than the normally narrow (10 km) boundaries seen in oceanic plates [86]. In the Pacific region, the Baja California Peninsula is moving northwest, separating from the rest of the continent; the oceanic plate of the Cocos is being assimilated by the continent in the southern Pacific of Mexico, from Cabo Corrientes in the state of Jalisco to Central America; this subduction occurs along an oceanic trench known as the Acapulco or Mesoamerican megashear [87]. Further, in the seismic regions of the Gulf of Mexico and the Caribbean there are geological forces of cortical separation (also referred as tension or distensive) operating on the continental limits, and as a result of the movement of the continental tectonic plates of North America towards the west and the Caribbean towards the east, they advance on the deepest bottoms of the oceanic basins. In the Pacific Plate, in southern California and in Baja California, the plates are migrating northwesterly relative to the North American plate along a series of transformation faults (San Andreas fault) connecting the extension centers, whose activity is gradually separating this territory from the rest of the continent, for which it will become an island in approximately 10 million years. Similarly, oceanic faults allow magma to escape, generating an expansion of the ocean floor [88,89].
The Rivera microplate is located in southern Baja California, right at the gateway to the Sea of Cortez where the magnetic lineaments of the ocean floor indicate how the gap between the Pacific plate and the Rivera plate, positioned between fracture zones, grows. Due to the movement that the Cocos and Rivera plates have towards the northeast of the Mexican Republic, a portion of these plates dips under the North American plate, causing great earthquakes to occur along the coast of Jalisco, Colima, Michoacán, Guerrero, Oaxaca, and Chiapas; yet we cannot determine if these large earthquakes were caused by the Cocos or Rivera plate movement. The Cocos plate is built in the Eastern Pacific mountain range, that goes from the Rivera fault zone to the Galapagos mountain chain. It is located off the coastlines of Michoacán, Guerrero, Oaxaca, and Chiapas, and it dips into the continental crust, leading to a displacement that causes earthquakes throughout the Pacific coast.
Cocos plate subducts beneath the North American plate at a rate of about 12 cm/yr from 20 Ma to 11 Ma and 6 cm/yr from 11 Ma to present [90]. Along the central portion of the Middle American Trench, the subduction interface shows a significant variety in along strike dip angles; also, the middle section of the plate, near Acapulco, has a horizontal slab. The trace/strike of the Neogene volcanic arc, which trends at a 17 degree angle to the trace/strike of the trench, demonstrates the along strike change in the dip. The dip is 50 degrees to the northwest near the Rivera Plate junction, and 30 degrees to the southeast near the Tehuantepec Ridge [91]. The slab in central Mexico has returned to its current location near the southern boundary of the Neogene volcanic arc, as shown by the southward movement of the volcanic arc [92].
A triple point to the southeast of the Tehuantepec Ridge divides the North American plate from the Caribbean plate, and the Cocos plate begins to subduct under it; this poses substantial natural hazards to most of central and southern Mexico. The Yucatan Peninsula, for its part, is rotating clockwise, and the Trans-Mexican Volcanic Belt is still active [93].

Natural Time Analysis Background: An Order Parameter for Seismicity and Its Minima
The natural time χ k for the occurrence of the k-th EQ of the energy Q k in a time series comprising N EQs is defined as χ k = k/N. Hence, the evolution of the pair (χ k , p k ) is studied in the NTA, where is the normalized energy and Q k is estimated by means of the relation [94] Q k ∝ 10 1.5M k , where M k stands for the EQ magnitude. The variance κ 1 = χ 2 − χ 2 of natural time χ weighted for p k , namely can be considered as an order parameter for seismicity [45]. The fluctuations of this order parameter of seismicity in an EQ catalog can be studied by using a fixed-length sliding natural time window containing a number W of consecutive EQs by means of the procedure described in References [9,95]. The window length W is selected to correspond to the average number of EQs that occur within the crucial scale [96] of a few months, or so, which is the average lead time of SES activities. To this end, we estimate all the κ 1 values from the subexcerpts of consecutive 6 to W EQs within the excerpt of the EQ catalog of W EQs and use them for the calculation of their average value µ W (κ 1 ) and standard deviation σ W (κ 1 ). The variability of the order parameter of seismicity κ 1 is given by [7,56] while its temporal evolution can then be pursued by sliding the natural time window of W consecutive EQs, event by event through the EQ catalog, see Figures 2 and 3. In such a procedure, we assign to β W the occurrence time of the EQ which follows the last event of the excerpt of W EQs in the catalog.

Earthquake Nowcasting and Earthquake Potential Score
Rundle et al. [11] proposed EQ nowcasting as a method for estimating the seismic risk through the current state of fault systems in the progress of the EQ cycle identified on the basis of natural time (cf. more recently the construction of time series resembling the EQ cycle has been extensively studied in References [16,18,97]). To estimate the seismic risk, EN uses an EQ catalog to calculate from the number of 'small' EQs, defined as those with magnitude M < M λ but above a threshold M σ , i.e., M ∈ [M σ , M λ ), the level of hazard for 'large' M ≥ M λ EQs. As mentioned in the Introduction, the EQ catalogs adopted [11,12,19,77,78,[78][79][80]98] are global seismic catalogs such as the Advanced National Seismic System Composite Catalog or the NEIC PDE catalog and for M σ the completeness threshold of the EQ catalog is usually selected [11]. Along these lines, the magnitude threshold M σ = 4.0 has been considered [12,78] for applications in areas such as Greece, Japan, and India that lie outside the United States. In EN, one employs the natural time concept and counts the number n of 'small' EQs that occur after a 'large' EQ-n stands for the waiting natural time or interoccurrence natural time. Then, the current number n(t) of the 'small' EQs since the last occurrence of a 'large' one is compared to the cumulative distribution function (CDF) of the interoccurrence natural time Prob[n < n(t)]. To estimate Prob[n < n(t)], it should be ensured [11] that we have enough data to span at least 20 or more 'large' EQ cycles. In EN, the EQ potential score (EPS) equals the CDF value, and measures the level of the current hazard, see Figure 4a. In References [11,12,[77][78][79] the seismic risk for various cities of the world was estimated through the following procedure: After calculating the CDF Prob[n < n(t)] within a large area, the numberñ of the small EQs around a city, i.e., those occurring within a circular region of epicentral distances r < R since the occurrence of the last 'large' EQ in this circular region, is found. Given that EQs exhibit ergodicity-see, e.g., References [99-101]-Rundle et al. [12] suggested that the seismic risk around a city can be estimated by using the EPS corresponding to the current value ofñ by inserting n(t) =ñ in Equation (4).    Table 1).
In the present study, the area considered is N 35 10 W 120 80 , which for M σ = 4.0 and M λ = 6.0 leads to the CDF shown in Figure 4a. When focusing on the period from 1 January 1989 to 1 January 2021, the empirical CDF comprises 218 EQ cycles and is shown in Figure 4a. In this figure, we observe that the fit using the Weibull distribution provides a fair approximation with root mean square of residuals [102] equal to 0.0177. This property of the Weibull distribution is in accordance with the results found by Pasari and Sharma [79] for Himalayan EQs as well as with those later found in Reference [19] for Eastern Mediterranean.  Figure 1 during the period 1989 to 2020. This equals to the EPS according to EQ nowcasting and has been calculated on the basis of 218 EQ cycles. The corresponding Weibull model fit [79] is also shown with the blue curve in lin-lin (a) and lin-log (b) diagram.

Construction of Average Earthquake Potential Score Maps
Upon the selection of a date, one can estimate by means of the EQ catalog the number of 'small' EQsñ that occurred since the last 'large' EQ inside a circle of radius R for each point in the study area. By insertingñ in Equation (5) we can obtain the corresponding EPS value. The EPS maps are constructed by evaluating first the value of EPS, EPS ij , at the points (x ij , y ij ) of a 'square' lattice that covers the study area and then by spatially averaging these EPS values within disks of the same radius R: where the summation is restricted to the lattice points whose distance d(x, y; x ij , y ij ) from the observation point (x, y) is smaller than or equal to R; N stands for the number of lattice points included in the sum.
In Section 3 of Reference [19], it has been shown that the EPS mean value m n (R, R ) where the summation is made over all the N ij lattice points (x ij , y ij ), estimated numerically in an EPS map drawn with self-consistency radius R and covering an area A of average radius R (≈ √ A/π) may take the form where d f is related to [19] the fractal dimension [103] of EQ epicenters.

Results
A self-consistent method of producing EPS maps using a radius R has been suggested and applied to Eastern Mediterranean in Reference [19], as already mentioned. The date on which β min has been observed -on the basis of the analysis of the properties of ENBOSAP [21]-has been selected for drawing these EPS maps. The latter are produced by first estimating EPS for disks of radius R centered at each point of a square 0.25 • × 0.25 • lattice covering the whole region of study and then by averaging these EPS values within the same radius R (see Section 2.4). A clear relation between the such made EPS maps and the epicenter of an impending strong EQ has been observed. Moreover, the stability of the method when changing the lattice 'constant' from 0.25 • to 0.10 • has been secured.
Here, we focus on the area N 35 10 W 120 80 of intense seismicity depicted in Figure 1. For example, during the 32 year period from 1 January 1989 to 1 January 2021 we have in total 19,184 EQs with M ≥ M σ = 4.0, i.e., approximately 50EQs/month, 28 of which have magnitude M ≥ 7.0, i.e., approximately 0.875 EQs/year (see Table 1). Due to this frequent occurrence of EQs of magnitude M ≥ 7.0, the application of the properties of ENBOSAP k-cores in order to identify the minima preceding these strong EQs as made in Reference [19] is not considered necessary. In other words, since on average every year we have approximately one strong EQ we just need to identify minima of the variability β W for a reasonable value of W as suggested by Varotsos et al. [96]. We selected W = 200 corresponding to the number of EQs that occur on average every four months.
This way, the dates of the minima β W,min identified six months before each strong EQ are shown by the vertical cyan lines in Figure 2 and are also inserted in the second column of Table 1. In the case of some strong EQs, i.e., those numbered 2, 17, 18, and 27 in Table 1, a minimum of β 200 cannot be identified clearly before the strong EQ, thus minima of β 150 and/or β 100 have been used. The detailed behavior of the variability of the order parameter of seismicity before each strong EQ can be seen in the excerpts shown in Figure 3. An inspection of the latter figure together with Table 1 reveals that some minima may correspond to more than one strong EQ. This does not impose any actual problem, because the EPS maps produced for such a date can be easily compared with the epicenters of the strong EQs that followed the minimum due to the large dimensions of the study area. Table 1. The dates of the variability minima β W,min for W = 200 observed within 6 months before all the strong M ≥ 7.0 EQs in the study area N 35 10 W 120 80 that have been selected for drawing the EPS maps shown in Figure 5. The lines corresponding to EQs of magnitude M ≥ 7.5 are indicated by typing the magnitude in bold. In the last line, the β 200,min , which was observed before the very recent 2021 Guerrero, Mexico EQ, is inserted (almost two weeks later the 2021 Jiquilillo, Nicaragua, M6.5 earthquake took place with an epicenter at 12.16 • N 87.85 • W). In Figure 5, using a 0.2 • × 0.2 • grid, we depict the EPS maps determined for R = 100 km for each β W,min date of Table 1 together with the location of the epicenters of the strong EQs that followed within the next six months. We observe that the EQ epicenters compare favorably with the contours of EPS . This relation will be further elaborated in the next section.

Discussion
We first studied the statistics of the EPS values closest to the epicenters of the strong EQs with M ≥ 7.0 during the period 1 January 1989 to 1 January 2021, i.e., the first 28 EQs of Table 1. The results for the mean value µ(R) and the standard deviation (STD) σ(R) are plotted versus R in Figure 6a, which indicates a characteristic behavior of a monotonically increasing µ(R). A further inspection of Figure 6a reveals that it is compatible with the functional form µ(R) = 1 2π log R 6.68 (6) − R 3580(50) (9) which includes the characteristic logarithmic singular part of the Green's function in two dimensions, see, e.g., Equation (9.96c) in paragraph 3 of Section 9.2.2.3 of Bronshtein et al. [104], together with a linear (non-singular) correction term. The presence of this singular logarithmic part when analyzing EPS closest to the strong EQ epicenters should be considered as a strong indication of the interrelation between the location of the epicenters and the corresponding EPS maps. It signifies that the future epicenters seem to act similar to sources in the two dimensional EPS field. The above point is further strengthened by the fact that when considering the EPS mean value EPS over all the lattice points of the grid (square lattice) for all the EPS maps of Figure 5 that correspond to the first 28 EQs of Table 1 the singular behavior disappears, see Figure 6b. The functional form observed in Figure 6b indicates a value of d f = 0.95 (2) which is compatible with the almost one-dimensional fault structure of Figure 1 (cf. the latter is dominated by the faults along the Pacific coast). Moreover, a comparison of Figure 6b with Equation (8) reveals that R = 2160(80) km, which compares favorably with the quantity √ A/π ≈ 2000 km when considering the study area A of 25 • × 40 • . Furthermore, we studied the statistics of EPS values closest to the epicenters of the 8 EQs with M ≥ 7.5 (shown with bold letters in Table 1 and Figure 5). It was found that for R = 250 km, a particularly interesting behavior is observed: Namely a bimodal distribution appears with mean value 0.55 and two lobes ±0.15 from the mean. This, in conjunction with the fact that at R = 250 km, Figure 6a reveals an average value of 0.50, encouraged us to investigate the maps of | EPS − 0.5| for values of R around 250 km (cf. in previous studies [19,58] R = 250 km has also been found of particular importance). Such a study revealed that | EPS − 0.5| maps may be useful and the optimum results (shown in Figure 7        The results shown in Figure 7 do not, of course, solve the very difficult problem of finding the epicenter of a future strong EQ since the areas covered by the tones of cyan color are spatially distributed in a rather large area, but clearly indicate that EPS maps definitely include information concerning the epicenter of a future strong EQ. Our efforts to improve these results are in progress by applying this method to other seismically active areas around the globe (cf. the operation of SES measuring stations [106] is an essential factor for such an improvement since additional information on the future epicentral area can be thus obtained).

Conclusions
Here, we studied the variability β W of the order parameter of seismicity κ 1 (introduced by natural time analysis) together with earthquake nowcasting within the highly active seismic region N 35 10 W 120 80 that covers Southern California, Mexico, and part of Central America. We suggest a self-consistent method of constructing EPS maps to obtain an estimation of the epicenter location of a future strong EQ of magnitude M ≥ 7.0 in this region. The study of EPS values closest to the strong EQ epicenters showed a logarithmic behavior, which is reminiscent of the Green's function in two dimensions. This is compatible with the view that the future epicenters act like 'sources' in these two dimensional maps. Using NTA and EN the epicenter of a future strong M ≥ 7.0 EQ was estimated to lie in the vicinity of a region covering on average only 8.5% of the total study area with a hit rate 28/29(≈96.5%). Data Availability Statement: Earthquake data come from the United States National Earthquake Information Center PDE and were downloaded from the United States Geological Survey, Earthquake Hazards Program [84]. The last date the data were accessed was 21 September 2021. Gnuplot [107] was used for the preparation of the figures. The coast lines were imported from GEODAS Coastline Extractor [108]. The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. Acknowledgments: Jennifer Perez-Oregon deeply appreciates support from the National Council of Science and Technology (CONACYT) of Mexico (CVU No. 376516).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: