1. Introduction
Analysis of the distribution of seismicity that falls within an area is important in determining the area where a strong earthquake is likely to occur in the future. In fact, some models, based on the analysis of seismicity [1] [2] [3], can identify areas of higher hazard, where strong earthquakes are more likely to occur.
Instrumental seismological network catalogs are the main tool from which to take data, which are analyzed by analysis algorithms, to derive the location of a probable future earthquake. The lack of data, causes large errors in predicting epicenters.
Our study starts from an assumption of proximity between the seismogenic source and the strong earthquake associated with it. This assumption can be verified by observing some seismic sequences (Figure 1).
Seismogenic sources are represented by faults that, subjected to strong tectonic stresses, can reach the rupture point, causing a strong earthquake, where the epicenter is located at the beginning of their rupture. Some studies [4] [5] [6] [7], have shown that asperities are areas of a fault plane with different rheological properties, where there is an increase in resistance to the rupture and fault movement is blocked compared to the rest of the surface that would continue to move.
The surface traces of seismogenic sources can be represented by a curve, characterized by concavities that, in our hypothesis, are ascribable to areas of asperity or earthquake clusters, where a future strong earthquake may most likely nucleate.
The model we propose is to analyze the distribution of seismicity related to earthquake sequences preceding strong earthquakes, through a polynomial function (Figure 2). The curve, which is derived, is a line characterized by the succession of Up and Down concavities and inflection points. It allows obtaining a band, where within it falls the strong earthquake of the analyzed sequence. Some cases studied show how the polynomial curve approximates the surface trace of the seismogenic source.
Figure 1. Location of the epicenters of some strong earthquakes. The red-colored star indicates the epicenter of the EQ, the red-colored line the trace of the seismogenic source, the black-colored arrow-vector indicates the dipping direction of the source surface.
Figure 2. Schematics of the geometry of the curve trace obtained from seismicity analysis with the polynomial function, with Up and Down concavity. The red-colored circles indicate the point of maximum and minimum concavity, the blue-colored circle indicates the point of bending, and the red-colored arrow-vector indicates the type of concavity Up or Down.
2. Data and Methods
2.1. Function Model
The proposed model is derived from the mathematical analysis of seismic sequences preceding strong earthquakes, through the use of a polynomial function. The analysis yields a curve showing a band where the strong earthquake of the analyzed sequence falls.
The model was calculated based on retrospective analysis of 180 earthquakes of magnitude M ≥ 7.0 and 24 between 6 - 6.9 M, occurring in various areas of the mode and using the catalogs of the National Institute of Geophysics and Volcanology (INGV) [8] and the U.S. Geological Survey (USGS) [9] and Japan (NIED) [10].
The number of earthquakes analyzed before the big shock varied from a minimum of 40 to a maximum of more than 1000 events (most frequently 80 - 120 events analyzed). The minimum magnitude and maximum depth used were 2.0 M and 50 km, respectively.
Analysis of the majority of seismic sequences associated with strong earthquakes shows that the mainshock epicenter is most frequently located within the concavity. In other words, the depth and width of the concavity define the area in which the strong earthquake is located. Other sequences, however, show how the strong earthquake is located in areas slightly outside the concavities.
Overall, the concavities and inflection points of a curve can be regarded as unstable zones, where most tectonic deformation accumulates that could contribute to nuclear a strong earthquake in the future.
Polynomial curves composed of a single Up or Down concavity are typical of areas where normal or reverse faults are present, while those composed of multiple Up and Down concavities are typical of transcurrent faults or transcurrent faults associated with reverse or normal faults (Figure 3).
Figure 3. Polynomial curve type and corresponding fault.
Seismicity data (Latitude and Longitude of earthquakes preceding a strong event) of a studied area, are analyzed by a polynomial function (1), which defines a curve where each of its points is described by the coordinate pair Xi and Yi (longitude and latitude of the i-th point).
In general, given the pairs (Xi, Yi) i = 0, …n, the solution of the polynomial interpolation is posed in the following terms: look for that polynomial
of degree n,
(1)
such that
(2)
where n is the degree of the highest term of the polynomial, at most equal to 10,
are the longitude and latitude of the analyzed earthquake series, and
is the coefficients of the calculation procedure.
By writing this relationship for each i, the following linear system is derived:
(3)
Then the coefficients αi i = 0, …, n of the polynomial P(x), which satisfies the conditions P(xi) = yi, i = 0, …n, are the solutions of the previous linear system.
Table 1 shows an example of polynomial curve calculation, using the last 20 earthquakes preceding the main event of Mw 7.8 in Türkiye on 06/02/2023.
Figure 4(a) shows the curve calculated with the polynomial function of degree 5. A similar result can be obtained using other methods to represent the interpolating polynomial (Figure 4(b)).
The epicenter of the mainshock was located near the inflection point of the polynomial curve. The inflection point and the two concavities represent the most dangerous areas and where a strong earthquake may occur.
Table 1. Polynomial curve calculation procedure. The epicenter data of the 20 earthquakes preceding the mainshock with Lat. (yi) and Long. (xi) preceding the strong earthquake in Türkiye on 06/02/2023, are given in ascending order according to Long. (Ascending data). The Calculate function columns show the longitude and latitude values of the polynomial curve points calculated with the Function model.
| Earthquake Data (source: USGS) |
Ascending data |
Calculate function |
| time |
Lat (yi) |
Long (xi) |
depth |
mag |
magType |
Long (xi) |
Lat (yi) |
Long (Xi) |
Lat (Yi) |
| 2023-0118T10:08:14.166Z |
38.4334 |
44.9458 |
19.102 |
5.7 |
mww |
25.8264 |
39.488 |
25.8234 |
39.47584 |
| 2023-01-18T10:42:39.023Z |
38.472 |
44.9169 |
16.66 |
4.3 |
mb |
28.4965 |
35.7086 |
28.4965 |
35.77241 |
| 2023-01-18T13:30:47.860Z |
36.9911 |
36.1694 |
10 |
4.1 |
mb |
31.2583 |
35.1548 |
31.2583 |
35.02887 |
| 2023-01-18T14:55:50.392Z |
38.3997 |
44.9477 |
10 |
4.4 |
mb |
35.7431 |
35.7692 |
35.7431 |
36.47892 |
| 2023-01-19T02:22:27.227Z |
39.488 |
25.8234 |
10 |
4.3 |
mwr |
36.1694 |
36.9911 |
36.1694 |
36.67677 |
| 2023-01-19T05:55:22.616Z |
38.5878 |
44.7612 |
10 |
4.3 |
mb |
36.4998 |
37.2003 |
36.4998 |
36.83167 |
| 2023-01-24T05:59:29.165Z |
35.1548 |
31.2583 |
10 |
4.1 |
mb |
40.3446 |
38.4311 |
40.3446 |
38.48665 |
| 2023-01-25T12:37:05.134Z |
35.7086 |
28.4965 |
29.59 |
5.9 |
mww |
44.717 |
38.5594 |
44.717 |
38.57678 |
| 2023-01-28T18:14:45.855Z |
38.4199 |
44.9097 |
16 |
5.9 |
mww |
44.7507 |
38.4819 |
44.7507 |
38.56077 |
| 2023-01-28T18:16:56.220Z |
38.4811 |
44.9282 |
10 |
5.1 |
mb |
44.7612 |
38.5878 |
44.7612 |
38.5557 |
| 2023-01-28T18:16:56.220Z |
38.4811 |
44.9282 |
10 |
5.1 |
mb |
44.8139 |
38.5153 |
44.8139 |
38.52973 |
| 2023-01-28T20:09:59.211Z |
38.6005 |
44.8287 |
10 |
4.5 |
mb |
44.8223 |
38.6392 |
44.8223 |
38.5255 |
| 2023-01-29T03:48:09.650Z |
38.6392 |
44.8223 |
10 |
4.1 |
mb |
44.8287 |
38.6005 |
44.8287 |
38.52227 |
| 2023-01-29T14:41:14.127Z |
38.4906 |
44.8693 |
10 |
4.7 |
mb |
44.8693 |
38.4906 |
44.8693 |
38.50142 |
| 2023-01-29T15:13:25.532Z |
38.4819 |
44.7507 |
10 |
4.5 |
mb |
44.8784 |
38.5428 |
44.8784 |
38.49667 |
| 2023-01-29T16:12:38.398Z |
35.7692 |
35.7431 |
10 |
4.3 |
mb |
44.9097 |
38.4199 |
44.9097 |
38.48012 |
| 2023-01-29T16:38:47.878Z |
38.5594 |
44.717 |
10 |
4.2 |
mb |
44.9169 |
38.472 |
44.9169 |
38.47626 |
| 2023-01-29T22:00:19.418Z |
38.5153 |
44.8139 |
10 |
4.2 |
mb |
44.9282 |
38.4811 |
44.9282 |
38.47018 |
| 2023-01-31T04:13:04.953Z |
38.5428 |
44.8784 |
10 |
4 |
mb |
44.9458 |
38.4334 |
44.9458 |
38.46061 |
| 2023-02-03T11:05:09.361Z |
37.2003 |
36.4998 |
10 |
4.2 |
mb |
44.9477 |
38.3997 |
44.9477 |
38.45957 |
| Mainshock |
| 2023-02-06T01:17:34.342Z |
37.2256 |
37.0143 |
10 |
7.8 |
mww |
|
|
|
|
![]()
Figure 4. Polynomial curve calculated using the last 20 earthquakes (Table 1) preceding the Mw 7.8 mainshock of Türkiye on 06/02/2024.
Figure 5 shows the results of the main calculation steps of the model using the last 1108 earthquakes preceding the main event of Mw 7.8 in Türkiye on 06/02/2023.
The calculation steps are summarized as follows:
Step 1. Calculate the mean epicenter of all earthquakes falling in the analyzed area (MEn);
Step 2. Divide the study area into four quadrants (I, II, III, IV) taking the mean epicenter MEn (blue color lines) as reference;
Step 3. Calculate the mean epicenter of each quadrant ME1,2,3,4 (purple colored circles);
Step 4. Locate the last five epicenters that occurred in the map (light blue, yellow, green, brown and black colored circles);
Step 5. Calculate the dynamic epicenters DEM11 and DEM22 [3] of the analyzed area (magenta and green color circles);
Step 6. Calculate the mean ME12 epicenter of the twelve epicenters identified above falling within the analyzed area (red filled circle);
Step 7. Order the longitude and latitude data of all epicenters in ascending order according to longitude (Ascending data in Table 1);
Step 8. Calculate and plot the polynomial curve of degree 5 (red dashed line);
Step 9. Draw the curves parallel to the polynomial curve passing through the four middle epicenters ME1,2,3,4, the midpoint MEn (purple color curves), the last five epicenters (blue, yellow, green, brown and black color curves). For points DEM11 and DEM22 (havana color curves) and the midpoint ME12 (red color dashed curve);
Step 10. Calculate the absolute and relative maxima and minima and the inflection points of the function and from them draw the normal lines to the function (Z1, Z2 and Z3).
Mathematical analysis with the polynomial function of the earthquake sequence leading up to the Mw 7.8 earthquake in Türkiye on 06/02/2024 yielded the interpolating curve and the bundle of red and orange curves passing through points DEM11, DEM22 and the last epicenter of the analyzed series (Figure 5).
Where these curves are closest, they represent the range within which there is the highest probability of occurrence of a strong earthquake.
In addition, it can be seen that the UP (Z2) and Down (Z1 and Z3) concavities identify unstable areas, where higher seismicity is observed and which according to our model may be the site of a strong earthquake.
The areas near the thirteen characteristic points (Steps 1, 3, 4, 5 and 6) calculated above are also those with a higher hazard of occurrence of a strong earthquake (red color box).
The most dangerous quadrant is generally where the dynamic epicenters DEM11 and DEM22 (quadrant I) and/or the middle epicenter ME12 (quadrant II) fall. The least dangerous is the opposite quadrant (quadrant III). The epicenter of the main event of Mw 7.8 falls within the areas obtained from the analysis with the polynomial function.
In Figure 6, straight lines perpendicular to the polynomial curve and passing through the characteristic points have been drawn. As can be seen, four danger zones (F1, F2, F3 and F4) can be identified. The strong earthquake of Mw 7.8 occurred in the vicinity of zone F3.
Analyses performed on many seismic sequences have shown that 97% of the epicenters of the future EQ lie near the concavities or mean inflection points of the calculated polynomial curve. Only 3% fall at the edges of the area, where there is little data to plot the curve of the polynomial function.
Figure 7 shows the results obtained by applying the model to some seismic sequences of earthquakes that occurred in various parts of the world.
Figure 5. Epicenter map of the 06/02/2023 Türkiye sequence region. Red and orange curved lines indicate the higher hazard zone. The purple-colored solid circle (MEn) indicates the midpoint. The quadrants are indicated by the letters I, II, III and IV in red color. The red-colored star indicates the epicenter of the strong earthquake on 06/02/2023. The solid circles of pink and green color DEM11 and DEM22. indicate the dynamic epicenters. The solid circles of blue, green, yellow, brown and light blue color indicate the last five epicenters of the occurred earthquakes, the solid circle of red color indicates the midpoint of the previous twelve points.
Figure 6. Epicenter map of the Türkiye sequence region of 06/02/2023. The red and black straight lines indicate the higher hazard bands (F1, F2, F3 and F4).
Figure 7. Function procedure applied to the Tajikistan 23/02/2023 (Figure 7(a)), Indonesia 26/12/2004 (Figure 7(b)), Italy 24/08/2016 (Figure 7(c)), Japan 16/03/2022 (Figure 7(d)) and Taiwan Region 03/04/2024 (Figure 7(f)) earthquakes. The pink and green filled circles indicate the dynamic epicenters calculated by DEM11 and DEM22 procedure, the red colored star indicates the mainshock epicenter.
2.2. Application Filters
The filters are designed to select concavities or inflection points identified by the function model that are most likely to be affected by a strong earthquake.
Through the filters, which can be applied alone or in combination, it is possible to control not only the data to be analyzed but also the data to be excluded in the analysis of seismic sequences.
Some application filters that can be used during the data analysis performed with the Function model are described below.
2.3. Circles Filter
This filter makes it possible to identify three circular areas (Figure 8) of higher hazard within the analyzed area and reduce the number of hazardous areas identified with the Function model.
The procedure to identify the circular areas can be summarized in the following steps:
Step 1. Draw from the MEn midpoint a circle (yellow circle) using the greatest distance between the MEn midpoint and the four quadrant midpoints (purple filled circles) as the radius;
Step 2. Draw from the DEM22 point the red-colored circle using the radius defined above;
Step 3. Draw from the DEM22 point (origin) the orange-colored circle using the distance between the last epicenter of the analyzed series (blue-colored filled circle) and the DEM22 dynamic epicenter as the radius.
Analyses of the seismicity within an area, show that the areas of higher hazard are the ones within the red circle. In addition, it was noted that the strong earthquake (7.8 Mw) most likely falls within the area between the red and yellow circles.
In Figure 9, F1, F2, F3 and F4 and the area between the red and yellow colored circles were drawn.
Figure 8. Procedure applied to the Türkiye earthquake seismic sequence of 06/02/2023. The map shows the circular areas of higher hazard (red and orange circles).
Figure 9. Epicenter map of the Türkiye region. Procedure applied to the earthquake sequence of the Türkiye earthquake of 06/02/2023. Red-colored areas indicate the higher hazard bands (A1, A2, A3 and A4), while yellow-colored areas are the lower hazard bands.
The areas of higher hazard (red colored areas A1, A2, A3 and A4), are those between the red and orange curved lines and the one passing through the last epicenter (blue colored solid circle).
Areas A3 and A4 falling within the red-colored circle, are those, where a strong earthquake may most likely occur, while area A3 between the two circles is the most dangerous, as evidenced by its proximity to the main earthquake of 7.8 Mw.
2.4. Reduction Filter
The proposed polynomial function model can be improved in terms of the thickness of the band, which defines the area with the highest probability of occurrence of a strong earthquake.
This filter starts with an initial estimate of the thickness of the band of curves, which is subsequently narrowed by eliminating those earthquakes that fall outside the identified area. The thickness of the band is refined iteratively, generating a new band of progressively thinner curves, which minimizes the error.
The optimal band is obtained as a result of following steps, until the number of interactions exceeds a threshold that does not change the thickness of the band of curves.
The calculation procedure is as follows:
Step 1. Apply the polynomial function to all epicenters in the analyzed series (Figure 10);
Step 2. Select the epicenters that fall within curves 1 and 2 (upper magenta and lower light blue lines) and calculate the polynomial function again;
Step 3. Repeat the procedure given in step 2 above, one or more times until a narrower range than the initial one is obtained.
Figure 11 shows the result obtained after five interactions. This procedure gives good results, when the epicenter of the mainshock falls near the trace of the polynomial function (dashed red color line) or in the middle part of the band of curves.
Figure 10. Epicenter map of Türkiye region. Function procedure applied to the earthquake sequence of Türkiye earthquake of 06/02/2023. The red dashed curve indicates the function trace, while number 1 and 2 indicate the upper and lower curves parallel to the polynomial curve.
Figure 11. Epicenter map of Türkiye region. Reduction procedure applied to the earthquake sequence of Türkiye earthquake of 06/02/2023. The red dashed curve indicates the trace polynomial function, while number 1 and 2 indicate the upper and lower curves parallel to the polynomial function obtained after five interactions.
2.5. Filter More Functions
The number of concavities and points of inflection increase as the degree of the polynomial function increases.
This filter allows the location of all the maxima, minima and inflection points of some sequences that are difficult to interpret with only the polynomial function of degree 5 to be identified more accurately by calculating and plotting the curves of degree 5, 6, 7, 8, 9, 10.
Looking at Figure 12, in which the curves calculated by varying the degree of the function are shown, five concavities can be seen that can be interpreted as more dangerous areas.
Applying the filters described above, we obtain that concavities 4 and 5 are those where a strong earthquake can most likely fall, as evidenced by the occurrence of the 7.8 Mw mainshock.
Figure 13 shows the curves calculated by varying the power of the function (5, 10, 15, 20, 25 and 30); six concavities are noted that can be interpreted as the most dangerous areas.
Using the filters described in the previous sections, we obtain that concavities 4, 5 and 6 are those where the epicenter of a strong earthquake may be located.
Figure 12. Procedure applied to the earthquake sequence of the 06/02/2023 Türkiye earthquake. The map shows the highest hazard band (red lines) and the maxima, minima and inflection points (indicated by progressive numbers).
Figure 13. Procedure applied to the earthquake sequence of the 06/02/2023 Türkiye earthquake. The map shows the highest hazard band (red lines) and the maxima, minima and inflection points (indicated by progressive numbers).
2.6. Inclined Seismic Band Filter
This filter makes it possible to identify certain areas of an area where an earthquake may occur using a simple procedure based on the alignments between the thirteen characteristic points calculated with the Function model.
The procedure for identifying areas can be summarized in the following steps:
Step 1. Divide the study area into four quadrants (I, II, III, IV) taking as reference the last epicenter of the analyzed series (blue color lines);
Step 2. Sort according to longitude the thirteen characteristic points calculated with the Function model;
Step 3. Identify the alignments formed by the thirteen characteristic points, comparing each point with the next ones, and calculate the angle αi of the best fit between the three points (Figure 14);
Step 4. Calculate the sum of the angles of the alignments falling in the first and third quadrants, the total number of angles and the average angle α1,3;
Stop 5. Calculate the sum of the angles of the alignments falling in the second and fourth quadrants, the total number of angles and the average angle α2,4;
Step 6. Draw the straight lines of mean angle α1,3 (Figure 15) passing through the thirteen characteristic points(red lines);
Step 7. Draw the lines of mean angle α2,4 passing through the last thirteen characteristic points (green color lines).
The epicenter of the expected earthquake is located within the outermost perimeter of the thickest red-colored box in Figure 15, obtained from the intersection of the outermost red- and green-colored straight lines drawn previously (steps 6 and 7), while the inner red-colored box may be home to the epicenter of the expected earthquake with a probability of 80%. Within this area falls the Mw 7.8 earthquake.
The joint application of this filter and the Circles filter makes it possible to reduce the amplitude of the area where a strong earthquake may occur, on which subsequent verifications can focus.
Figure 16 shows the projections of the characteristic points on the narrowest band, and the most dangerous areas are indicated with boxes. Boxes B, C, and E falling in the concavities obtained with the Function model indicate the most dangerous areas, as evidenced by the strong earthquake of Mw 7.8.
The analysis of seismicity preceding the Greece earthquake of 08/06/2008 (Figure 17), shows that the most dangerous area is represented by the innermost red box, as evidenced by the location of the epicenter of the strong Mw 6.4 earthquake.
Figure 14. Sloping band filter. Red-colored lines indicate alignments with angle αi falling in the first and third quadrants, and green-colored lines indicate alignments with angle αi falling in the second and fourth quadrants. Solid circles of different colors indicate the thirteen characteristic points.
Figure 15. Slant band filter The red and green colored lines indicate the seismic bands. The red-colored areas indicate the regions where a strong earthquake may occur, the middle box indicates the most dangerous area, the red and orange dashed curved lines indicate the circle traces obtained with the Circles filter.
Figure 16. Sloping bands filter. The red and green color lines indicate the seismic bands. Red-colored boxes indicate areas where a strong earthquake may occur.
Figure 17. Epicenter map of the region of Greece. Procedure applied to the earthquake sequence of the Greece earthquake of 08/06/2008. Blue color circles indicate the location of the thirteen characteristic points and red and green color lines the seismic bands. The thicker red color lines indicate the region where a strong earthquake may occur, respectively, while the middle box indicates the most dangerous area.
2.7. Fault Filter
This filter can be useful when you want to analyze the seismicity that develops along a main fault before the mainshock.
It is possible to select the epicenters closest to the fault using a selection box or band (Figure 18(a)) and then apply the polynomial function to process the curve and the filters above to select the concavity or inflection point.
Below are the results obtained by applying the “faults” filter to the magnitude M 7.3 earthquake of 5/28/2009 (Figure 18(b)), which struck the Atlantic coast of Honduras and whose epicenter was located near a left-trending fault [11] [12].
In the analyzed territory, there are two main faults: the Mid America Trench on the southern sector and the left transcurrent fault on the northern sector
Figure 18. Epicenter map of the Honduras region and main tectonic elements. Procedure applied to the seismic sequence of the Honduras earthquake of 05/28/2009. The red color indicates the main faults. Blue and green-colored lines indicate the data selection boxes and bands (Figure 18(a)), while blue-colored curves (Figure 18(b)) indicate the trace of the polynomial functions.
(Swan fault).
We divided the territory into two boxes (a and b) and analyzed the seismicity that preceded the earthquake with the Function model.
The analysis produced similar polynomial curves composed of two concavities (UP and Down) and an inflection point.
The epicenter of the earthquake was located near the Down concavity closest to the trace of the left transcurrent fault.
3. Conclusions
Understanding the spatial distribution of the epicenters of earthquakes that occurred in a region before a strong earthquake is a good procedure for identifying the area where a strong earthquake is most likely to occur.
The analysis of the seismicity of an area, by means of the polynomial function, makes it possible to derive a curve characterized by one or more Up and Down concavities that connect at the bending points. Strong earthquakes from the sequences analyzed show that the epicenter lies within the range defined by the polynomial curve and to a good approximation near the concavities and flexures. In addition, the application of filters makes it possible to identify more accurately the areas where the strong earthquake of the analyzed sequence falls.
The polynomial curve of our model, in many cases analyzed tends to approximate well the surface trace of the seismogenic source, which represents the fault that generated the strong earthquake.
In our model, the polynomial curve, which approximates the seismogenic source, is represented by a series of concavities that are ascribable to areas of asperity or earthquake clusters, where a future strong earthquake may most likely nucleate.
In some cases, such as the Central Italy earthquake, it is observed that the polynomial curve does not approximate the seismogenic source because of a strong dispersion of epicenters with respect to the mainshock location, due to the presence of multiple fault segments in the analyzed area.
In such cases, in order to obtain a good approximation, an initial filtering of epicenters should be performed, selecting those that fall near the trace of the seismogenic source.
The proposed model can be used to identify areas where a strong earthquake may occur.
In those cases, where the polynomial curve approximates the surface trace of the seismogenic source, the UP and Down concavities indicate the areas of the fault surface where more deformation accumulates due to tectonic stress, and which act as blocked areas compared to the surrounding areas.
These areas could experience increased seismicity with the possible nucleation of a strong earthquake.
By analyzing seismic sequences with relatively little data, the graphs produced by the model contain enough information to identify the most likely areas of occurrence of a strong earthquake.