Estimating the Epicenter of an Impending Strong Earthquake by Combining the Seismicity Order Parameter Variability Analysis with Earthquake Networks and Nowcasting: Application in the Eastern Mediterranean

The variance κ1 of the natural time analysis of earthquake catalogs was proposed in 2005 as an order parameter for seismicity, whose fluctuations proved, in 2011, to be minimized a few months before the strongest mainshock when studying the earthquakes in a given area. After the introduction of earthquake networks based on similar activity patterns, in 2012, the study of their higher order cores revealed, in 2019, the selection of appropriate areas in which the precursory minima βmin of the fluctuations β of the seismicity order parameter κ1 could be observed up to six months before all strong earthquakes above a certain threshold. The eastern Mediterranean region was studied in 2019, where all earthquakes of magnitude M≥7.1 were found to be preceded by βmin without any false alarm. Combining these results with the method of nowcasting earthquakes, introduced in 2016, for seismic risk estimation, here, we show that the epicenter of an impending strong earthquake can be estimated. This is achieved by employing—at the time of observing the βmin—nowcasting earthquakes in a square lattice grid in the study area and by averaging, self-consistently, the results obtained for the earthquake potential score. This is understood in the following context: the minimum βmin is ascertained to almost coincide with the onset of Seismic Electric Signals activity, which is accompanied by the development of long range correlations between earthquake magnitudes in the area that is a candidate for a mainshock.

The EQs exhibit power law correlations in time, space, and magnitude (M)  that have been suggested [70,71] to indicate the proximity to a critical point [72][73][74][75][76]. Several proposals suggest [73,74,77] that the EQ dynamics are associated with a second-order phase transition, or with a mean-field spinodal [76] that can be understood as a line of critical points. A first-order phase transition (see, e.g., [78] and references therein) has also been proposed as an alternative. In phase transitions, the new phase corresponds to the mainshock and a quantity-termed order parameter is required to correctly identify the new phase. Landau and Lifshitz [79] state: "To describe quantitatively the change in the structure of the body when it passes through the phase transition point, we can define a quantity η, called the order parameter, in such a way that it takes non-zero (positive or negative) values in the unsymmetrical phase and is zero in the symmetrical phase".
Along these lines, Varotsos et al. [23] proposed the variance κ 1 of the NTA of EQ catalogs as an order parameter for seismicity. As such, it is non-zero before the mainshock, falling abruptly to zero upon its occurrence (see, e.g., the last points in Figure 6 by Márquez et al. [80]) and its scaled distribution exhibits a behavior similar to that of other equilibrium and nonequilibrium critical systems [7,21,23,28,81]. Later, it was shown that, when studying the EQs in a given area, the fluctuations β of this order parameter were minimized a few months before the strongest mainshock [82], and the experimental results revealed that such minima occurred almost simultaneously with the observation of Seismic Electric Signals (SES) activities [83], which are series of low frequency ( f ≤1 Hz) electric signals that precede strong EQs [84][85][86][87][88][89][90]. Subsequent studies have shown that such minima are observed before strong EQs, both regionally [29] and globally [91]. These minima were found to be statistically significant precursors of EQs [92][93][94]. Furthermore, Varotsos et al. [95] showed that regional β min in Japan is observed when long-range correlations develop between EQ magnitudes, which has also been verified in two recent strong EQs in the Americas [40,96].
In recent decades, the theory of networks [97][98][99][100][101] has triggered many applications in seismicity (see, e.g., [47,60,[102][103][104][105][106][107][108]), which showed that the EQ network properties may reveal a seismic hazard [109][110][111][112][113][114]. Tenenbaum et al. [60] introduced the EQ networks based on similar activity patterns (ENBOSAP) which are very useful for the NTA of seismicity [83]. Furthermore, the study of the higher order cores of these networks allowed the selection of appropriate areas in which the precursory minima β min of the fluctuations of the seismicity order parameter κ 1 can be observed up to six months before all strong mainshocks above a certain threshold [115]. In particular, the case studied in 2019 was the eastern Mediterranean region and resulted in five precursory β min before all EQs of magnitude M ≥ 7.1 without false alarms.
It is the scope of the present paper to show that by combining the above two NTA results with the modern method of nowcasting EQs [33], an estimate of the epicenter of the impending EQ with M ≥ 7.1 in the eastern Mediterranean region is achieved. EQ nowcasting employs natural time, which is unique in its characteristics [116], to estimate seismic risk by means of an EQ potential score (EPS) and found many useful applications, both regionally and globally [33,[116][117][118][119][120][121][122][123][124][125][126]. EQ nowcasting [33] has, so far, focused on describing the current state of fault systems.

Study Area and EQ Data
Here, the area of interest was the eastern Mediterranean, see Figure 1. We used the United States National Earthquake Information Center (NEIC) PDE catalog-the data available from the United States Geological Survey (USGS), cf. [127]-in the region N 52 23 E 50 5 , and considered all EQs recorded during the period of almost 45 years from November 1976 to February 2021 with M thres = 3.8. The eastern Mediterranean Sea and its surrounding coastal areas of Greece and western Turkey are one of the most seismically active and rapidly deformed regions on the continents [128]. This area is tectonically dominated by the relative plate motions between Eurasia and African lithosphere plates and their microplates. The African plate moves north to Eurasia from Gibraltar to the east. The plate boundary crosses North Africa and continues into southern Sicily, continuing around the Adriatic through Italy and eastern Balkans to Albania and northwestern Greece. There is a suggestion that the Adriatic is part of the African plate forming a small plate called Apulia [129]. The Apulia plate has a continental structure, e.g., it has a similar mean density to the density of the Eurasian plate, and their interaction leads to a continental collision. This microplate is suggested to rotate counter clockwise and its interaction with Eurasia is responsible for the shallow EQs in the region of Croatia, Albania and northwestern Greece. East of Apulia, the rapid movement of two small plates, the Aegean and the Turkish, takes up the relative movement between the Eurasia and African plates. The Aegean plate consists of the Aegean Sea, part of the Greek peninsula, Crete, and part of the western area of Asia Minor, while the Turkish part includes most of Asia Minor and Cyprus. The collision between Arabia and Eurasia in the Caucasus and eastern Turkey forces the Turkish plate to move westward in relation to Eurasia. Due to this movement, the strain accumulates along the North Anatolian transform fault (e.g., the opposite sides of this fault belong to different plates that move horizontally on top of each other) that accommodates large destructive EQs. The northern front of the African plate, south of the Adriatic Sea, forms the eastern Mediterranean oceanic lithosphere. As this oceanic lithosphere has a higher density than the continental lithosphere of the Aegean plate, it sinks below it, with an average angle of 38 • , and forms the Hellenic Arc (SW of Cefallonia to S. Peloponnese, S. Crete up to S. Rhodes) and is associated with a sharp crustal shortening and an uplift rate of a few mm/yr along the Hellenic Arc due to the accretion of sediments of the African plate below the overriding Aegean plate [130]. Most of the subduction movement is accommodated along the Hellenic arc causing shallow depth EQs at its outer part. The seismically active regions of western Turkey, eastern and northern Greece, and the northern Aegean Sea are dominated by extension forces such as the southern Aegean Sea move in an almost southwesterly direction relative to Eurasia. This extension is thought to be driven principally by the "roll-back" of the subducting slab below the southern Aegean, as it sinks into the mantle [131]. The geophysical properties are different in the Aegean Sea from those in the outer part of the Hellenic arc. Another geotectonic feature causing destructive EQs is the significant right-lateral strike-slip movement along the North Aegean trough, and the Cephalonia-Lefkada transform zone, which is due to the offset between the oceanic-continental convergence to the west and the westward propagation of the Anatolian plate to the east. The combination of all the above factors leads to a significant left-lateral movement along the southeastern front of the Aegean-African interface [130].  [132].  46 5 is depicted with the (red) links and nodes. The red rectangle shows the geographical area (A 6 ) that marginally encloses this 6-core, i.e., N 50 28 E 44 7 . We note that the corresponding 5-core, hereafter labeled A 5 , covers the whole region N 50 25 E 46 5 and is bounded by the green rectangle.

Minima of the Variability of the Order Parameter of Seismicity
In a time series comprising N EQs, the natural time χ k for the occurrence of the k-th EQ of energy Q k is defined as χ k = k/N. This means that the time intervals between consecutive events are ignored, but their order, as well as their energy, Q k are preserved. In NTA, the evolution of the pair (χ k , p k ) is studied, where: is the normalized energy and Q k is estimated by means of the relation [133] Q k ∝ 10 1.5M k , where M k stands for the EQ magnitude. It has been argued [23], as mentioned in the Introduction, that the variance κ 1 = χ 2 − χ 2 of natural time χ weighted for p k , namely: can be considered as an order parameter for seismicity. The study of the fluctuations of this order parameter of seismicity in an EQ catalog was conducted through a fixed-length sliding natural time window containing a number W of consecutive EQs that would occur on average within the crucial scale [82] of a few months, or so, which is the average lead time of SES activities. Following the procedure described by [29,95], we estimated κ 1 from the subexcerpts of consecutive 6 to W EQs within the window of W EQs and, then, calculated the average value µ W (κ 1 ) and the standard deviation σ W (κ 1 ) of the obtained κ 1 values. The quantity: was termed [4,134] variability of κ 1 . The temporal evolution of the β W value ( Figure 2) could then be pursued by sliding the natural time window of W consecutive EQs, event by event, through the EQ catalog and assign to its value, the occurrence time of the EQ which follows the last EQ of the window in the EQ catalog. The corresponding minimum value was labeled β W,min . The quantity β W of Equation (3) is reminiscent [115] of the square root of the Ginzburg criterion introduced in the context of the mean field theories of phase transitions of the Ginzburg-Landau type (see, e.g., p. 175 of Goldenfeld [135]). The latter criterion was also discussed by Holliday et al. [72].
To identify a precursory variability minimum [29,91,95], two window lengths, W and W (W > W), corresponding to the number of EQs that occurred within a few months, were selected and β W and β W (sliding independently event by event) should have exhibited their local minima almost simultaneously. This meant that one required at least 90% of the EQs included in the calculation of the local β W minimum, and were also included in the calculation of the local β W minimum. Once this condition was valid, it was examined whether the ratio r ≡ β W ,min /β W,min laid within the margins defined by the precursory to strong EQ variability minima; see, e.g., [29,91,95]. For example, for the area A 6 of Figure 1, the selection of r ∈ (r 1 = 0.98, r 2 = 1.07) led to 7 variability minima (see Table 2 of [115]) when imposing a threshold β 0 = 0.253 below which the minima of β W for W = 80 must be considered by the aforementioned procedure (cf. W = 80 EQs corresponds to the number of EQs occurring on average in area A 6 of Figure 1 within two months).

Earthquake Networks Based on Similar Activity Patterns
Tenenbaum et al. [60] suggested the construction of the ENBOSAP as follows: The area under study was evenly spaced in geographical coordinates in "squares" of dimensions of about 100 km × 100 km, each marked by two integers (i, j). The center of each "square" is a possible node of the network. The fact that EQs do not occur on a daily basis in every square makes it necessary to coarse-grain the seismic activity over a predefined time interval. There is a trade-off between precision and data richness leading [60] to an optimal coarseness in the 90-day time series. Thus, the seismic activity s t in each square covered 90 days and s t+1 corresponded to the next 90 days of non-intersecting time periods, giving about 4 different values per year. The energy released in each 90-day period inside the "square" with indices (i, j) was given by [60]: where N t (i, j) is the number of EQs displayed in the 90-day t-th time window with epicenter inside the square (i, j) of magnitude M l t (i, j) above a given threshold M thres , i.e., M l t (i, j) ≥ M thres . To define a link between the centers of two squares, say (i, j) and (i , j ), one calculated the Pearson product-moment correlation r x,y between the time series x t ≡ s t (i, j) and y t ≡ s t (i , j ): where . . . indicates the average value of . . ., and σ x and σ y are the standard deviations of x t and y t , respectively. The two grid squares were linked if r x,y was greater than a specified value r c , which is a tunable parameter. Mintzelas and Sarlis [115] found that, for NTA applications of ENBOSAP, r c should be chosen as the one that determines the upper 5% percentile of the positive r x,y among all the different squares; for the eastern Mediterranean region, the value r c = 0.7 was selected. As usual, the degree k of a node was the number of links the node had. Tenenbaum et al. [60] showed that ENBOSAP: (a) displays links without a characteristic length scale with much more links than expected only by chance alone, (b) is much more assortative than random ones, and (c) exhibits significant stability over time. These findings corroborated the existence of high-degree nodes with similar activity patterns over a distance of more than 1000 km away; see, for example, Figure 1.

Earthquake Nowcasting
Rundle et al. [33] proposed EQ nowcasting as a method for estimating the current state of fault systems in the sense that it determines the progress in the EQ cycle. To achieve this determination, one uses an EQ catalog to calculate from the 'small' EQs, which are defined as those with magnitude M < M λ , but above a threshold M σ , i.e., of magnitude M ∈ [M σ , M λ ), the level of hazard for 'large' M ≥ M λ EQs. The lower threshold M σ is typically the completeness threshold of the EQ catalog used [33], while the EQ catalogs employed [33,116,116,117,123,124,136] are global seismic catalogs such as the Advanced National Seismic System Composite Catalog or the NEIC PDE catalog. For such global catalogs, a M σ = 4.0 magnitude threshold was considered [116,117] for applications in EQ-prone areas outside the United States, such as those in Greece, Japan, and India. By employing the natural time concept, one counts the number n of small EQs that occur after a large EQ. In other words, n is the waiting (natural) time or interoccurrence natural time. The current number n(t) of small EQs, since the occurrence of the last strong EQ, is compared to the cumulative distribution function (CDF) of the interoccurrence natural time Prob[n < n(t)]. For such an estimate of Prob[n < n(t)], it should have been ensured [33] that we had enough data to span at least 20 or more strong EQ cycles. In EQ nowcasting, the EQ potential score (EPS) equals the CDF value, and is a measure of the level of current hazard, see Figure 3a. In References [33,116,117,123,136], 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 strong EQ in this circular region, was found. Given that EQs exhibit ergodicity-see, e.g., References [58,137,138]-Rundle et al. [117] 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 (6).
In the present study, the large area considered was N 52 23 E 50 5 which, for M σ = 4.0 and M λ = 6.0, led to the CDF shown in Figure 3a. A study of the empirical CDF comprising 123 EQ cycles in Figure 3a revealed that the fit with the Weibull distribution provided a good approximation with a root mean square of residuals equal to 0.0207; see Chapter 15 of [139]. This property of the Weibull distribution was in accordance with the results found by Pasari and Sharma [123] for Himalayan EQs.  Figure 1. The results came from 123 EQ cycles after using the Weibull model; see Table 1 of [123].

Results
The importance of high-degree nodes of ENBOSAP in the NTA of seismicity was first suggested in Reference [83]. It was observed that the inclusion of high-degree nodes in the geographical area in which we applied the NTA of seismicity was important for the clearer observation of β W,min before strong EQs. In order to determine the important highdegree nodes that should have been included in the geographical areas where we studied seismicity in natural time, Mintzelas and Sarlis [115] searched for k-cores in ENBOSAP. Recall that a k-core in a network is the largest subgraph whose vertices have a degree at least k [99]. Thus, in a k-core, each node has at least k nearest neighbors within that subgraph; see, e.g., the six-core depicted in red color in Figure 1. It is noteworthy that a network is hierarchically organized as a set of successfully enclosed k-cores [99], similar to the Russian nesting doll "Matrioshka". In other words, the five-core covered an area larger than the six-core and, as a network, it included the six-core itself; see Figure 1b,c of Reference [115]. The NTA of seismicity improved when we comparatively examined [115] the results from the areas N 50 25 E 46 5 (labeled A 5 ) and N 50 28 E 44 7 (labeled A 6 ) of the five-core and the six-core of the area N 50 25 E 46 5 , respectively, in a fashion similar to that followed in [95]. More precisely, when focusing on the minima β W,min , which were simultaneously observed in both areas A 5 and A 6 , we found that they preceded within six months, without any false alarm, all strong EQs with M ≥ 7.1 in the smaller region, i.e., A 6 (see the red rectangle in Figure 1). The results for the variabilities β W for both areas A 5 and A 6 for W = 80 and 120 EQs, corresponding to the number of EQs that occurred on average there for the periods of (approximately) two and three months, respectively, are shown in Figure 2 for the almost forty-year period from 1 January 1981 to 9 February 2021. A careful inspection of this figure showed that no additional minima were observed simultaneously since October 2018, which was the end of the period studied by [115]. In the same figure, with (red) open circles we marked the simultaneous β 120,min that preceded the strong EQs with M ≥ 7.1 in area A 6 , as mentioned. The dates corresponding to these minima are given in Table 1. In  the same Table, we also inserted the date of another local minimum in both areas A 5 and A 6 which did not meet the criteria, β 80 < β 0 = 0.253 and r 1 < r < r 2 , mentioned in the Materials and Methods to be considered as precursory, but it could be useful to provide an estimate of the epicenter of the Düzce EQ. Table 1. The dates of the five almost simultaneously (within 3 months) observed minima in the limiting geographical areas A 5 and A 6 of the 5-core and the 6-core of the area N 50 25 E 46 5 , respectively, together with the strong EQs that they had preceded. a They are marked with (red) open circles in Figure 2. In the last line, a local minimum that was also observed before Düzce EQ was also added; see Figure 2. Since the aforementioned minima provided an estimate of the occurrence time (i.e., within the subsequent six months) of the M ≥ 7.1 EQs without any false alarm, we now proceeded to an estimate of the epicenter location of the impending EQ. The NTA already provided an answer to this question [31,32] based on the simultaneous observation of β minima calculated in the large area and in various small areas. The latter was found to cluster around the epicenter location of the impending EQ. Such a method, however, required detailed seismological catalogs for a large geographical region such as those provided by the Japan Meteorological Agency. Unfortunately, in the case of the eastern Mediterranean region studied here, such a catalog does not yet exist and, hence, this method could not be applied at present. To partially overcome this lack, we employed the following procedure: We applied the EQ nowcasting to circular areas of radius R with centers located at the points of a grid of 0.25 • × 0.25 • and estimated for each center an EPS value based on the Weibull approximation of the CDF depicted in Figure 3a. A self-consistent way to obtain an estimate of the epicenter location of the impending EQ would then be to average the aforementioned calculated EPS values within a circular area of radius R around the point P and assign this average EPS, EPS , to it. In this way, we could obtain for R = 250 km the EPS maps shown in Figure 4, which were calculated at the date of the minimum β W,min (see the first column of Table 1) preceding the M ≥ 7.1 EQs whose epicenters were depicted by the (black) open squares. Comparing the EPS maps with these epicenters, we saw a reasonable agreement.
A model to understand the process of self-consistently averaging EPS is shown in Figure 5; to keep this model simple, we assumed that all EQs occurred practically in the center O of the yellow-colored circular region "Y" of radius R. After a long enough period without a large EQ, the EPS took the value of unity for all points inside "Y" and, since we assumed that no other EQ took place, the EPS would be zero outside it. Hence, the value of EPS at a point P lying at a distance d away from O simplified to: where A(d, R) is the (red) line-shaded area of the two intersecting circles shown in Figure 5, and equaled [140]: Hence, the value of EPS in polar coordinates (ρ, θ) with center at O equaled to: for ρ < 2R and zero otherwise. If we now calculate the mean value of EPS : in a circular region of radius R > 2R, we simply obtain that: by virtue of the integrals 2R 0 arccos ρ 2R ρdρ = πR 2 /2 and The validity of Equation (11) was verified by assuming that EQs occurred according to the aforementioned simple model and, using the computer programs used for the calculation of the EPS maps shown in Figure 4, the numerically found m n (R, R ) had the form:  Figure 2) was used. The latter minimum was not as deep as those identified in Reference [115]; see the discussion there about Düzce EQ, which was followed by the deadly 7 September 1999 M6.0 Athens EQ and was discussed later here.
with d f = 1.986 (5), which was very close to two, and R = 1961 (14) km, which was also very close to the radius (1913 km) of an area corresponding to the approximately 11.5 × 10 6 km 2 covered by the area N 53 23 E 50 10 used for the construction of the EPS maps. Additionally, when the quantity m n (R, R ) was studied for the EPS maps of Figure 4, the corresponding quantities were d f = 1.32 (6) and R = 1320(80) km. Interestingly, this value of d f differed only slightly from the value d f ≈ 1.2, which Bak et al. [46] found to describe the fractal dimension of the location of epicenters projected onto the surface of the Earth in a unified scaling law obeyed by the distribution of waiting times between EQs occurring in California and ranging from tens of seconds to tens of years. d P O Y Figure 5. The value of EPS was considered unity inside the yellow-colored circular region of radius R with center at O-this disk was called region "Y" in the text-and zero outside. The EPS was calculated for a point P lying at distance d from the center of the yellow-colored region and equals to the ratio of the (red) line-shaded region over the area πR 2 of the circle having its center at P.

Discussion
An inspection of Figure 4 revealed that the epicenters of the impending EQs correlated favorably with the EPS maps. Of course, M ≥ 7.1 EQs did not occur in areas where EPS was highest, but rather in regions where EPS was rather mediocre. Detailed results showed that the EPS values at the grid points nearest to the epicenters laid in the range from 19% to 50%. Namely, this EPS value equaled 22% for Lesvos EQ, 28% for Vrancea EQ, 21% for Aqaba EQ, 19% for Izmit EQ, 30% for Düzce EQ, and 50% for Van EQ. This result was probably related to the fact that, at the time of β W,min , long-range correlations started to develop in the magnitude time series [40,95,141].
Interestingly, the above results were consistent with respect to the two free parameters of the proposed method (i.e., the grid size and the radius R) were concerned. If we considered the five EQs with M ≥ 7.1 for which simultaneous minima were identified at the same time (i.e., see the first five rows of Table 1) and studied the statistics of EPS of the closest to the epicenters grid points versus R, we found that their mean exhibited a plateau centered at R ≈ 250 km, which was also accompanied by an inflection point for the corresponding standard deviation (STD); cf. Figure 3b. This result was obtained for a 0.1 • × 0.1 • grid and a similar behavior was found for the 0.25 • × 0.25 • grid used for drawing the maps of Figure 4. Indeed, the selection of R = 250 km in the latter figure was determined on this basis. Note that R = 250 km was also used in the NTA of seismicity of Japan [31] for estimating the epicenter of an impending M ≥ 7.6 EQ. At this point, it should be noted that when using the empirical CDF shown in Figure 3a instead of its Weibull approximation, the results did not change significantly and the maps remained almost the same.
In view of the aforementioned ability of the EPS maps calculated at the time of β W,min to provide information on the impending EQ epicenter for M ≥ 7.1, we proceeded to the calculation of such maps upon the occurrence of β W,min that, although were not simultaneous in areas A 5 and A 6 , were followed by an EQ of M ∈ [6.5, 7.0] within six (or at the most seven) months. The dates of such minima are listed in the first column of Table 2, and were used for the construction of the EPS maps shown in Figure 6. Figure 6 was drawn by using a 0.25 • × 0.25 • grid as in Figure 4. In the last column of Table 2, the EPS estimated at the grid points nearest to the impending EQ epicenter are shown (cf. they had an average value 38% and an STD of 20%). A similar calculation revealed an average of 38% and an STD of 18% when a grid 0.1 • × 0.1 • was used, pointing to the stability of the calculation. Table 2. The dates of minima observed in the limiting geographical areas of either the 5-core or the 6-core of the area N 50 25 E 46 5 that preceded the EQs with M ≥ 6.5, but were not simultaneous in both areas. a On the last row, a local minimum that was also observed before the deadly Athens EQ was also added, see Figure 2. In the last column, we also inserted the EPS value at the grid point nearest to the epicenter when a 0.25 • × 0.25 • grid was used as in Figure 6.  Tables 2 and 3 of [115].
Of course, the present method is a long way from becoming a really practical one, and we are still in the progress of improving it, since the cover of areas in Figures 4 and 6 was too large to find the final epicenter of an impending strong EQ. It provides, however, insight to the very difficult question of forecasting strong EQs. At this point, it is worthwhile to add that, in the case of Japan, where better quality seismological data exist, we could obtain a better estimation of a future EQ epicenter location by means of an NTA [31,32].
The two cases of the local minima during September 1999, which corresponded to the last row of Table 1 and the last row of Table 2, concerning Düzce and Athens EQs, respectively, were of particular regional interest. These two minima were not reported before, as they corresponded to values β 80 and β 120 , which were both well above 0.253. Their existence as local minima could be verified in Figure 2, while the corresponding EPS maps shown in Figures 4 and 6 correlate well with the epicenters of the deadly EQs of Düzce and Athens in 1999. Interestingly, the minimum of β 80 on 3 September 1999 roughly coincided with the observation of an SES activity near Lamia, Greece, which was related to the deadly Athens EQ; see Chapter 17 of Lazaridou-Varotsos [142] (see also p. 105 of Reference [90]).
Moreover, a close inspection of Figure 2 revealed that local minima of β 80 in area A 6 coincided with the SES activities recorded on 17 March 2001 [143], 7 November 2007 [144], and 2 October 2018 [136]. These results were remarkable, as they provided additional evidence to support the conclusion [83] that SES activities were accompanied by clear variations of seismicity, giving rise to simultaneous precursory phenomena in independent geophysical observables. This behavior was, of course, more obvious in the case of mega-EQs such as the 2011 M9.0 Tohoku EQ in Japan. In particular, it was found that the deepest minimum β min of the order parameter κ 1 of seismicity during the 27-year period from 1 January 1984 to the occurrence of the 2011 Tohoku EQ occurrence on 11 March 2011 was observed on 5 January 2011 [29], which agreed with the experimental finding that β min was observed [83] when an SES activity was initiated. This was so because the date of 5 January 2011 overlapped with the observation [145][146][147] of the abnormal magnetic field variations on the z-component from 4 to 10 January 2011 at two measuring sites-see Figure 1 of [148]-laying at epicentral distances of around 130 km, a fact which was characteristic of the existence of an SES activity [4].  Figure 6. Same as Figure 4, but here, the EPS was calculated at the time of the β 120 minimum reported in Tables 2 and 3 of [115] for all 'false alarms' preceding the strong M ≥ 6.5 EQs within the area A 5 . The EPS evaluated (at the time of the local minimum of β 80 ) on 3 September 1999 before the deadly M6.0 Athens EQ is also shown; see Figure 2. Note that, although a strong EQ activity was evolving at Izmit at that time, a strong EQ struck the Greek capital.

Conclusions
Here, we showed that upon computing maps of the average earthquake potential score in a self-consistent manner, we could obtain information on the epicenter location of an impending strong EQ. This was achieved by employing earthquake nowcasting at the time of the observation of the minimum β min of the fluctuations of the order parameter of seismicity. Such a minimum almost coincided with the SES activities, which were characterized by the development of long-range correlations, and became more easily detectable when a natural time analysis was employed in the areas covered by the higher order cores of the EQ networks based on similar activity patterns. Data Availability Statement: Earthquake data came from the United States National Earthquake Information Center PDE and were downloaded from the United States Geological Survey, Earthquake Hazards Program [127]. The last date the data were accessed was 9 February 2021. Gnuplot [149] was used for the preparation of the figures. The coast lines were imported from GEODAS Coastline Extractor [150]. 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 the 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: