Long-Term Forecasting of Strong Earthquakes in North America, South America, Japan, Southern China and Northern India With Machine Learning

Strong earthquakes (magnitude ≥7) occur worldwide affecting different cities and countries while causing great human, ecological and economic losses. The ability to forecast strong earthquakes on the long-term basis is essential to minimize the risks and vulnerabilities of people living in highly active seismic areas. We have studied seismic activities in North America, South America, Japan, Southern China and Northern India in search for patterns in strong earthquakes on each of these active seismic zones between 1900 and 2021 with the powerful mathematical tool of wavelet transform. We found that the primary seismic activity patterns for M ≥ 7 earthquakes are 55, 3.7, 7.7, and 8.6 years, for seismic zones of the southwestern United States and northern Mexico, southwestern Mexico, South American, and Southern China-Northern India, respectively. In the case of Japan, the most important seismic pattern for earthquakes with magnitude 7 ≤ M < 8 is 4.1 years and for strong earthquakes with M ≥ 8, it is 40 years. Every seismic pattern obtained clusters the earthquakes in historical intervals/episodes with and without strong earthquakes in the individually analyzed seismic zones. We want to clarify that the intervals where no strong earthquakes do not imply the total absence of seismic activity because earthquakes can occur with lesser magnitude within this same interval. From the information and pattern we obtained from the wavelet analyses, we created a probabilistic, long-term earthquake prediction model for each seismic zone using the Bayesian Machine Learning method. We propose that the periods of occurrence of earthquakes in each seismic zone analyzed could be interpreted as the period in which the stress builds up on different planes of a fault, until this energy releases through the rupture along faults and fractures near the plate tectonic boundaries. Then a series of earthquakes can occur along the fault until the stress subsides and a new cycle begins. Our machine learning models predict a new period of strong earthquakes between 2040 ± 5 and 2057 ± 5, 2024 ± 1 and 2026 ± 1, 2026 ± 2 and 2031 ± 2, 2024 ± 2 and 2029 ± 2, and 2022 ± 1 and 2028 ± 2 for the five active seismic zones of United States, Mexico, South America, Japan, and Southern China and Northern India, respectively. In additon, our methodology can be applied in areas where moderate earthquakes occur, as for the case of the Parkfield section of the San Andreas fault (California, United States). Our methodology explains why a moderate earthquake could never occur in 1988 ± 5 as proposed and why the long-awaited Parkfield earthquake event occurred in 2004. Furthermore, our model predicts that possible seismic events may occur between 2019 and 2031, with a high probability of earthquake events at Parkfield around 2025 ± 2 years.

Strong earthquakes (magnitude ≥7) occur worldwide affecting different cities and countries while causing great human, ecological and economic losses. The ability to forecast strong earthquakes on the long-term basis is essential to minimize the risks and vulnerabilities of people living in highly active seismic areas. We have studied seismic activities in North America, South America, Japan, Southern China and Northern India in search for patterns in strong earthquakes on each of these active seismic zones between 1900 and 2021 with the powerful mathematical tool of wavelet transform. We found that the primary seismic activity patterns for M ≥ 7 earthquakes are 55, 3.7, 7.7, and 8.6 years, for seismic zones of the southwestern United States and northern Mexico, southwestern Mexico, South American, and Southern China-Northern India, respectively. In the case of Japan, the most important seismic pattern for earthquakes with magnitude 7 ≤ M < 8 is 4.1 years and for strong earthquakes with M ≥ 8, it is 40 years. Every seismic pattern obtained clusters the earthquakes in historical intervals/episodes with and without strong earthquakes in the individually analyzed seismic zones. We want to clarify that the intervals where no strong earthquakes do not imply the total absence of seismic activity because earthquakes can occur with lesser magnitude within this same interval. From the information and pattern we obtained from the wavelet analyses, we created a probabilistic, long-term earthquake prediction model for each seismic zone using the Bayesian Machine Learning method. We propose that the periods of occurrence of earthquakes in each seismic zone analyzed could be interpreted as the period in which the stress builds up on different planes of a fault, until this energy releases through the rupture along faults and fractures near the plate tectonic boundaries. Then a series of earthquakes can occur along the fault until the stress subsides and a new cycle begins. Our machine learning models predict a new period of strong earthquakes between 2040 ± 5 and 2057 ± 5, 2024 ± 1 and 2026 ± 1, 2026 ± 2 and 2031 ± 2, 2024 ± 2 and 2029 ± 2, and 2022 ± 1 and 2028 ± 2 for the five active seismic zones of United States,

INTRODUCTION
Different natural phenomena like the fall of meteorites, tsunamis, volcanic eruptions, droughts, ice ages, the reversal of geomagnetic field, forest fires, droughts, earthquakes, and others can pose a significant danger and threat to human life and humanity's economic developments and resource managements (Murray, 2021).
Earthquakes are caused not only by natural seismic and tectonic processes but often time can also be induced by various anthropogenic activities such as nuclear bomb detonations, large dams, and subsurface exploitation of natural resources. The danger and risk posed by usually low intensity earthquakes induced by anthropogenic activities can be indeed mitigated by reducing or completely stopping the human activities that are responsible by these types of minor earthquakes. In a sharp contrast, especially earthquakes of great intensity that are caused by natural processes cannot be avoided but only forewarned with their often catastrophic and damaging impacts minimized.
Different sources and mechanisms have been suggested as triggers and modulators of earthquakes (see, for example Batakrushna et al., 2022, for a full review). For example, even the Sun's activity has been suggested as a significant agent causing earthquakes (Anagnostopoulos et al., 2021). Other proximate causes discussed in the literature include pole tide (Shen et al., 2005), pole wobble (Lambert and Sottili, 2019), surface ice and snow loading (Heki, 2019), glacial isostatic rebound (Hampel et al., 2007), heavy precipitation (Hainzl et al., 2006), atmospheric pressure (Liu et al., 2009), sediment unloading (Calais et al., 2016), seasonal groundwater change (Tiwari et al., 2021), seasonal hydrological loading (Panda et al., 2020). In addition, the Earth's rotation and tidal spinning have also been suggested as driver of plate tectonic activity.
The present geological paradigm about solid Earth is the plate tectonic theory which describes that the lithosphere is segmented into a series of plates that are in constant motions due to mantle mobility or convection. As a result of their interaction, a series of geological, mainly convergent and divergent, processes take place at their plate margins, ranging from seismicity, orogenic processes, and volcanism. The World Stress Map (WSM) 1 compiles the orientation of maximum horizontal stress (σ Hmax ) where we delimited our study areas in Figure 1 (Heidbach et al., 2016).
The dynamics of the plate tectonics provide a framework to understand the evolutive shape and dynamics of the earth's surface. Plate boundaries involve either divergence, like at oceanic spreading centers and continental rifts, or convergence, such as subduction (ocean to continent or ocean to ocean) and collision zones with different angles of displacement ranging from orthogonal towards subparallel one.
Only minor cases involve transform boundaries that facilitate plate kinematics on the global sphere. These boundaries accommodate plate-parallel relative displacement by strike-slip motion on vertical or steeply dipping faults. Due to these frictional contacts between the different types of plates, seismicity is triggered, producing a succession of earthquakes that progressively decrease in intensity in increasingly distant/ remote areas away from the seismic center/zone.
The sliding between tectonic plates is quite varied. Some plates slide without any consequences on Earth's surface, while catastrophic failures punctuate others. Also, after a few hundred meters some earthquakes stop. Nevertheless, others continue to collapse even after thousands of kilometers (Kanamori and Brodsky, 2004).
The driving mechanisms of plate tectonics remain not well unknown or poorly understood. Are they due to internal factors or external astronomical forces? We are hoping that the analysis of seismic patterns could provide some clues and information about the sources and mechanisms that are responsible for both tectonic movements and earthquakes.
Earthquake forecasting is one of the most difficult areas of research even though it is clear that its early prognosis can save many lives (Jain et al., 2021). Deterministic prediction of the exact coordinates of the epicentre, its depth, magnitude and exact time of one earthquake at the moment remains difficult and possibly impossible (see, for example, Shcherbakov et al., 2019;Beroza et al., 2021). Ogata (1988) suggests that the seismic pattern and temporal variation are usually very complicated. Furthermore, temporal seismic clustering is complex and difficult to discern or anticipate in advance. Different models have been proposed to analyze space-time clusters of seismicity in a region. One example is the Epidemic Type Aftershock Sequence (ETAS) model. This model suggests that the earthquake of a particular magnitude (M) in a region during a period of time can be approximately 1 http://www.world-stress-map.org/casmo/ considered as a Poisson process (Ogata, 1988). In addition, the method of the minimum area of alarm for earthquake magnitude prediction (Gitis and Derendyaev, 2020) and a method for earthquake predictions based on alarms (Zechar and Jordan, 2008) have all been suggested and evaluated.
The studies of earthquake precursors such as observing crustal geochemical fluids and gases, ultra-low frequency magnetic signals, atmospheric effects including ionospheric total electron content measurements, and several recording seismicities in regions experiencing earthquakes in terms of atmospheric, geochemical, and historical information can all help to improve and refine earthquake prediction (see for example, Pulinets and Boyarchuk, 2005;Ouzounov et al., 2018;Pulinets and Ouzounov, 2018). Since 2007, the Collaboratory for the Study of Earthquake Predictability (CSEP) has actively conducted and rigorously evaluated earthquake forecasting experiments as well as the prospective evaluations of earthquake forecast models and prediction algorithms (see, for example, Schorlemmer et al., 2018). CSEP's main targets and focuses are to optimize earthquake forecasting, advance forecast model development, test model hypotheses, and improve seismic hazard assessments.
The medium-term prediction of the strongest earthquakes has been carried out by the M8 algorithm, which is an algorithm for evaluating times of increased probability (TIPs) for strong earthquakes (Keilis-Borok and Kossobokov, 1990) from intensity of an earthquake flow and rate differential on a specific seismic region of earthquake source concentration and clustering. Also, the prediction of extreme events such as earthquakes demonstrates the efficiency and potential of the algorithms based on a pattern recognition approach as example the M8 algorithm (Kossobokov andSoloviev, 2008, 2018). In addition, the M8 algorithm shows that the hypothesis that the largest earthquake events are mere random variations in seismically active regions can be confidently rejected (Kossobokov and Soloviev, 2021). Kossobokov et al. (2015) suggested that "forecasting earthquake information must be reliable, tested, confirmed by evidence, and not necessarily probabilistic". We disagree slightly with this opinion and interpretation of Kossobokov et al. (2015). Probabilistic forecasts in the last century have provided new results to understand natural phenomena (see, for example Wigner, 1967;Landau and Lifshitz, 1988b;Feynman et al., 2011b). In this work, we show the results of a Bayesian model of Machine Learning, which is a probabilistic model. We do agree with Kossobokov et al. (2015) that all forecasts which are either probabilistic or not probabilistic must indeed be confirmed by evidence. We think that only future events will show if our probabilistic Machine Learning predictions are on the right track or not.
In recent years artificial intelligence (AI), deep learning (DL), machine learning (ML) (see, for example, Essama et al., 2021) have been applied to earthquake forecasting. In particular the use of ML in the study of earthquakes has been implemented in the detection, arrival time measurement, phase association, location and characterization (Beroza et al., 2021). In addition, the use of ML has focused on forecasting the exact magnitude of the next strong earthquake in different seismic zones (see, for example, Yousefzadeh et al., 2021).
In this paper, we propose a new method of analysis and algorithm for forecasting strong earthquakes (i.e., magnitude ≥7). We suggest that one promising progress to earthquake forecasting may consist in changing the prediction paradigms from an "exact" approach to probabilistic forecasting of future seismic activity cycles. This work aims to find the temporal seismic patterns of high and low seismicity in four major seismic zones: 1) the United States and Mexico, 2) South America, 3) Japan, and 4) Southern China and Northern India as sketched in Figure 1. We have made a probabilistic long-term earthquake prediction using a Bayesian ML model in these seismic zones based on the seismic patterns deduced from our wavelet analyses.

SEISMIC STUDY ZONES
We are analyzing the earthquake activity records in four major seismic zones in this work, and in turn we have made the probabilistic prediction for large earthquake (magnitudes M ≥ 7). The probability density function (PDF) of the spatial coordinates of all earthquake zones has been calculated for each seismic zone analysed. The PDFs of the longitudes and latitudes are shown at the top and left panels in Figures 4,9,13,15. • Seismicity related to transform and subduction margins in North America (Southwestern regions of both United States and Mexico) The scope of the tectonics setting related to the seismicity along Mexico's Pacific coast is divided into two regions: Northern and Southern (Figure 4).
In the northern Mexican subduction zone, the Gulf of California spreading center and the triple junction point around the Jalisco and the Michoacán Blocks represent the most active seismogenic belts inducing significant seismic hazard in the Jalisco-Colima-Michoacán region (Dañobeitia et al., 2016). The oblique to sub-parallel motions between the North American and Pacific plates at the latitude of the San Andreas fault produce a broad zone of large-magnitude earthquake activity mainly associated with dextral strike-slip faults extending more than 500 km into the continental interior. This seismic and tectonic activity patterns define the western limits of plate interaction as well as dominate the overall pattern of seismic strain release (Castro et al., 2017). Due to the Rivera Transform Fault, this seismic source corresponds to the shallow seismicity (mean depth value of 16 km), showing strikeslip faulting mechanisms delimiting the southwestern border between the Rivera and the Pacific Plates (Sawires et al., 2021).
Most of the earthquakes in Southwestern Mexico are due to the subduction process between the Cocos and North American tectonic plates (Pardo and Suárez, 1995). The subduction zone extends 1,300 km along the coast of the Pacific Ocean from the Chiapas state to the Jalisco states showing an angle in the range of 12°-45° (Suárez et al., 1990;Singh and Pardo, 1993). The Rivera plate moves with a relative velocity between 2.5 and 5.0 cm/year (Kostoglodov and Bandy, 1995), and the relative velocity of the Cocos plate in the Pacific Coast is in the range of 5.0-7.5 cm/year (Singh and Pardo, 1993). The subduction earthquakes are originated on the Pacific Coast with reverse fault focal mechanism and depths in the range of 10-40 km. The rupture lengths of these earthquakes are between 50 and 250 km, and their widths vary between 75 and 150 km (García et al., 2005). The deeper, in-slab earthquakes are also related to the subduction process, and their epicentres are located inside the continent. The hypocenters of these events have occurred at depths between 40 and 150 km in Mexico's central and western zone, and they are produced by the rupture of the subducted lithosphere (Jara et al., 2015).
• Seismicity in South America (subduction oceanic lithospheric plate versus continental lithospheric plate, Andean case) The South American plate is bounded for the subduction of the oceanic Nazca plate towards the west and the South Atlantic crustal oceanic section of the plate towards the east (Figure 9). This westerly subduction and the easterly spreading stresses due to the opening of the Oceanic Middle Ridge produces a compressional stress pattern on all the continents (Figure 1).
Earthquakes along the Andean cordillera show a progressively deep from the Chilean trench towards east associated with reverse (majority) or strike-slip faulting mechanisms with the principal significant compressional stress (σ 1 ) roughly in E-W direction. Along the Pacific coast, the deadliest earthquakes associated with tsunamis were registered during the last decades.
The seismicity in central Chile observed from 0 to 30 km depth beneath the western Andean thrust is due to the subduction of the Nazca plate. It shows essential seismicity beneath the Principal Cordillera located at a depth of 10 km, and deeper seismicity (~15 km) aligned with the main Andean thrust (Ammirati et al., 2019). Rivas et al. (2019) determine in the foreland Andean region that has 44 seismic locations with focal depths mechanisms showing mainly reverse and in less proportion strike-slip solutions ranging between 10 and 30 km and magnitudes 1.2 ≤ M ≤ 4. The intermediate principal stress (σ 2 ) is also compressional and more significant than the lithostatic pressure (σ v ). In the mid-plate South America, earthquakes seemed related to purely compressional stresses pattern (both σ Hmax and σ hmin larger than σ v ). Along the Atlantic margin, the regional stresses are affected by coastal effects due to transform Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 905792 fault stresses as well as flexural effects from sediment load at the continental margin (Assumpcao et al., 2016). This coastal effect tends to make σ Hmax parallel to the coastline and σ hmin (usually σ 3 ) perpendicular to the coastline. Few available breakout data and in-situ measurements are consistent with the pattern derived from the earthquake focal mechanisms (Heidbach et al., 2016). The Rio de la Plata craton and surrounding areas of Argentina and Uruguay located on the Atlantic margin of the South America plate have always been known as having deficient and shallow earthquakes activity related to preexisting faults (Rossello et al., 2020).
• Seismicity in Japan (subduction oceanic lithospheric plate versus oceanic lithospheric plate) Japan's most densely populated area is subjected to intense crustal stress due to the convergence by subduction of two oceanic lithospheric plates ( Figure 12). This compressional context produces consequently high seismic activity, among other geological processes associated with the subduction, such as volcanism and exhumation (Heidbach and et al., 2016). There are approximately 5,000 minor earthquakes recorded per year, mainly between 3.0 and 3.9 magnitudes and around 160 earthquakes with a magnitude of 5 or higher, which caused significant damages or casualties.
Geological analyses previously performed in the field, such as geomorphological markers, trench surveys and radiometric dating along the Nojima fault, associated with the 6.9 magnitude 1995 Kobe earthquake, revealed significant activity during the Pleistocene-Holocene times. Mainly, at least two significant earthquakes preceding the 1995 earthquake occurred in the last 1800 years (Lin, 2018). According to these authors, there would be no consistency between the recurrence intervals of seismic events proposed in previous contributions.
In this region, among many other areas of interest, the Hinagu-Futagawa fault zone (HFFZ) represents a study case, where the 7.1 magnitude Kumamoto earthquake occurred in 2016, which produced a surface rupture of~40 km in length. Detailed geological studies in the field and radiometric dating (Lin et al., 2017) allowed inferring four events in the last 4,000-5,000 years on this fault, suggesting a mean late Holocene recurrence interval of 1,000 years for the associated morphogenic earthquakes. However, as expressed by the authors mentioned above, these results contradict previous studies estimating recurrence intervals of 3,600-11,000 and 8,000-26,000 years for the target segments of the Hinagu and Futagawa faults, respectively. Different methods have been carried out to predict seismic events in this region (Uyeda, 2013), albeit with numerous difficulties inherent to the complexity of the discipline and the non-linear nature of such complex chains of phenomena.
• Seismicity in Southern China-Northern India (collisional margin continental lithosphere plate versus continental lithosphere plate, Himalayan case) The seismicity of this zone ( Figure 15) corresponds to the sutured obducted margin formed by the collision of the Indian plate against the Asian plate, which is occurring for the last tens of millions of years.
As a result of this collision, the high Himalaya orogen was formed (Tapponnier et al., 1982). As the process is still active even today, intense crustal N-S oriented stress (Heidbach et al., 2016) associated with intense seismic activity and cortical deformation were produced (Figure 1).
Numerous contributions (Bilham, 2019) have been analyzed through different methodologies, rupture zones and rupture propagation directions, regular convergence rates, as well as evaluated the slip potential in different segments of the Himalaya and the occurrence of potential high magnitude earthquakes. There is increasing evidence that Himalayan seismicity may be bimodal: blind earthquakes (up to M~7.8) tend to cluster in the lower part of the seismogenic zone, while infrequent large earthquakes (Mw 8+) propagate up the Himalayan frontal thrust (Dal Zilio et al., 2019).
Recently, Michel et al. (2021), considering many variables and uncertainties, suggest that earthquakes of magnitude greater than 8.7 (such as the one that occurred in 1950) are the most likely candidates to be the largest possible earthquakes in the Himalayas. However, they emphasize that, given the magnitude frequency distribution model used, the probability of a magnitude 8 + earthquake occurring in 100 years ranges from 60-80%. The most likely associated recurrence time for such an event exceeds 1,000 years.

Spatial Clustering
For this work, the orographic basemap comes from the ArcMap map library. The ocean bathymetry layer is obtained from the General Bathymetric Chart of the Oceans (GEBCO) 2 . The seismic records of the period 1902-2021 of the North America, South America, China and Japan regions taken from the U.S. Geological Survey (USGS) 3 , are processed in a geographic information system (GIS) to create a vector layer of geo-referenced points. Continuous surface maps with seismic magnitude information are designed with the Kriging-type interpolation method within a GIS. This probabilistic method is widely used to generate seismic maps (Türker and Bayrak, 2018;Teves-Costa et al., 2019;Moradia et al., 2020) to its efficiency in predicting information from one variable to the other. Through the spatial structure of its discrete values. Eight classes with equal intervals (0.5) are defined to better represent and interpret the spatial results. The processing and plotting of seismic magnitude data on a map is done using GIS. Plate Boundary and Movement Information from USGS was added to the interpolated map to find out its geographic location and plate type, as well as its displacement. A spatial data filter is applied to the seismic record layer to extract values ≥7 for inclusion in the interpolated map to locate the largest earthquakes. Using spatial analysis tools, density zones from the seismic records are defined for the filtered data.

Seismic Events
In order to analyze and search for any coherent seismic patterns, it is necessary to have an excellent catalogue of seismic activity. Seismic patterns should show periods when there is seismic activity as well as when there is no seismic activity at a certain magnitude if it exists. In this work, we have analyzed the seismic activity in four different and important seismic zones: 1) Seismicity in North America (Figure 4) 2) Seismicity in South America (Figure 9) 3) Seismicity in Japan ( Figure 12) 4) Seismicity in Southern China-Northern India ( Figure 15) The public data on seismic activity have been obtained from USGS. These data contain the list of all registered seismic activity. The table also contains information on the date of each seismic event, its magnitude, depth, type of magnitude, longitude, and latitude. We will analyze seismic activity for strong earthquakes magnitudes ≥ 7 available from 1900 to 2021. Also, we only analyze the earthquakes that are delimited in the areas/zones shown and discussed in Figures 4, 9, 12, 15.

Time Series With Data Gaps
One of the biggest problems of analyzing incomplete time series (such as the seismic activity catalogues) is to extract the information of the phenomena. A solution used regularly to analyze this type series is to apply different interpolation methods Sturges (1983). However, any interpolation can lead to an underestimation of spectral power at both higher and lower frequencies. Another technique sometimes used is to extract the mean value of the data (Carroll et al., 1997). Though, this technique can fail for data records that have a trend.
In particular, gaps in geological, geophysical, geographic, and seismic databases exist because those records contain errors, or were not complete nor homogeneous, and were created with data from different epochs (Jopek and Kaňuchová, 2017;Soon et al., 2019). In this case, a Bayesian block in the geosciences record can be applied to suppress the inevitable corrupting observational errors (see, for example, Scargle et al., 2013). The statistical inference using a Bayesian approach is used for analyzing an incomplete database (Gelman and Meng, 2005). In addition, spectral analysis has also been used to study time series with missing data (Maoz et al., 1997;Ding et al., 1998;Ramírez-Rojas et al., 2019). For example, the Fourier Transform (FT) is applied for analyzing these data. Nevertheless, this method may not often be suitable for non-stationary and irregularly spaced time series (Velasco Herrera V. M. et al., 2022). Therefore a new and reliable method is required to study a time series with gaps such as the seismic events.
The classical wavelet technique (see, e.g., Torrence and Compo (1998); Velasco Herrera et al. (2017)) is used to analyzed nonstationary times series, but this technique can only be used for regularly spaced time series. This is why a modified wavelet technique has been proposed and used to analyse incomplete time series. This technique has been named gapped wavelet (Frick et al., 1997;Soon et al., 2019). In our study, we have used the gapped wavelet spectral algorithm to analyze seismic events. Another variant for the analysis of irregularly spaced time series with wavelet has been proposed in Salcedo et al. (2012) recently.

Gapped Wavelet
In order to find seismic patterns in an area, it is necessary to analyze the periodicities that the seismic time series could have. Since earthquakes do not occur continuously, this implies frequent inactive or zero-activity "gaps" in the seismic catalogues. That is why we propose to analyze the seismic catalogues with the gapped wavelet transform. The gapped wavelet transform (W g ) of a time series with data gaps f g (t) is a matrix (Ω) and is defined by Frick et al. (1997) as: where t is the time index and a is the wavelet scale, the superscript (*) indicates the complex conjugate, and ψ is the mother function. We applied the Morlet's mother function (ψ) to analyze the power spectral density (PSD) of seismic activity since this mother function does not only provide a higher periodicity resolution but also is a complex function that allows calculating the inverse wavelet transform (Torrence and Compo, 1998;Velasco Herrera et al., 2017;Soon et al., 2011Soon et al., , 2019. The meaningful wavelet periodicities with a confidence level greater than 95% must be inside the cone of influence (COI), and thin black contours mark the interval of 95% confidence (Torrence and Compo, 1998). The global spectra show the power contribution of each periodicity inside the COI. Also, we established the significance levels in the global wavelet spectra with a simple red noise model that increases power with decreasing frequency (Gilman et al., 1963). The uncertainties of every peak position are obtained from the peak full width at half maximum (Mendoza et al., 2006).

Inverse Wavelet Spectral Analysis
The decomposition of f g in channel (y n ) can be obtained from the inverse wavelet (Torrence and Compo, 1998) as: where j 1 and j 2 define the scale range of the specified spectral bands, ψ o (0) is an energy normalization factor, C δ is a reconstruction factor, and δ j is a factor for scale averaging. For the Morlet wavelet, δ j = 0.6, C δ = 0.776, and ψ o (0) = π −1/4 . The input data in the W g are the seismic catalogues between 1900 and 2021. The W g has two main outputs as shown in Figures 2, 5, 7, 10, 13; the global spectrum, which shows the periodicities existing in the seismic record with the 95% confidence level above the red noise spectrum drawn as red dashed line (left panel) and the wavelet power spectral density (PSD) that shows the evolution over time of these periodicities (central panel).

Non-Linear Autoregressive eXogenous Model
The system's state can describe the dynamics of seismic activity from its input (V)-output (Y) behavior, which describes its evolution over time. Different models can approximate the state of the system. In particular, we use the Non-linear Autoregressive eXogenous (NARX) (Suykens et al., 2005) model in order to create forecasting models of seismic activity variation (Ŷ) that is defined as:Ŷ where Ξ is a transfer function of the state of the system at the moment "k" to the moment "k + 1", that depends intrinsically on the input and output data (V and Y, respectively), p and Q are the delay times, andŶ is the estimated seismic activity at a time "k + 1". We used the Least-Squares Support-Vector Machines (LS-SVM) algorithms to estimate the transfer function (Ξ), a nonlinear function (Suykens et al., 2005) as: where D is training data, in our case the seismic records. Also D k denotes the input data, i.e., the seismic records at time "k" (discrete time index from k = 1, . . . , n), ω k is the weighting factor which in turn has functional dependence on V k and B is the bias term. We use the Bayesian inference ML model (Suykens et al., 2005) obtained from the seismic records to provide a probabilistic earthquake prediction of the variation in the seismic activity. Bayes's theorem (Bayes, 1763) is the basis of our ML model and can be expressed as follows: where Ξ is the Least-Squares Support-Vector Machines (Eq. 9) regression model. Furthermore, Bayes's theorem is used to deduce the optimal parameters in the LS-SVM model (Eq. 9). In addition, we use the radial basis function (RBF) kernel in this LS-SVM method. In this work, we have applied and modified the LS-SVM algorithms and toolbox by Suykens et al. (2005).
3.6.2 Algorithms for the Estimation of the Next High or Active Phase of Seismic Activity We apply the following iterative steps to forecast the next high seismic activity season: 1) Use wavelet transform (Eq. 1) to find the periodicities (seismic patterns) in each seismic zone analyzed. The results are shown in Figures 2, 5, 7, 10, 13.
2) The decomposition of the seismic record in time series called "channels" with the periodicities obtained in step (1) can next be obtained using the inverse wavelet (Eq. 7). 3) Selection of the model lags P and Q for each Bayesian inference model that has been analyzed in each seismic zone. 4) Use the Radial Basis Function (RBF) kernel. 5) For training, validation, testing and deduction of the hyperparameters of the model. Use the K-fold cross-validation. 6) Set aside 1/K of data. Train the model with the remaining (K − 1)/K data. Measure the accuracy obtained on the 1/K data that we had set aside. K independent training is therefore acquired. The final accuracy will be the average of the previous K accuracies. Note that we are hiding or withholding a 1/K part of the training set during each iteration. This is applied at the time of training. After these K iterations, we obtain K accuracies that should be similar to each other; this would be an indicator whether the model is working well or not. In this work, K = 10 is adopted, but, it is possible to vary K between 5 and 10. 7) Determination of the weight and bias. 8) Estimation of next high cycle of earthquake activity using Eq. 9. Before forecasting the following period of strong earthquake activity, it is necessary to quantify the ability of the Machine Learning to "predict" the recent clustered earthquakes. We use 80% of the Bayesian clustered model (that is, data from 1900 to 1996) as input data to "forecast" the remaining 20% of the Bayesian clustered model (i.e., 1997-2021) in each seismic zone analyzed. The Bayesian clustered model of the historical earthquakes shows that all the historically strong earthquakes were manifested during the positive phase of the Bayesian clustered model; this fact indicates that our model has no overtraining nor undertraining. Furthermore, the wavelet analysis shows that the high and low seasons of strong earthquakes have a multiannual and multidecadal variations, so the Bayesian model we deduced is not overly complex, which implies that the validation is simple. We do not show the validation figures but instead choose to concentrate on the forecasting result. 9) Computation of a cost function. 10) Test of the accuracy of the estimate next high cycle of earthquake activity. 11) Test of the cost function: if this function was small enough, we stopped. Otherwise, we change one of the parameters and repeat from step (2) onwards.
We have used and modified the LS-SVM algorithms and toolbox by Suykens et al. (2005) for this goal.
We want to add that the accuracy of any forecast of seismic activity is limited by an uncertainty principle (Velasco Herrera et al., 2015). Great precision in the spatial location forecast implies a significant uncertainty in the temporal forecast. This is why in this work, we focus on the problem of temporal forecasting by proposing a new Bayesian Machine Learning composite method. In our case, the possible zones where the following strong earthquakes in each seismic zone analyzed could occur have been essentially clustered and pre-determined according to the methodology described in Section 3.1.

RESULTS
In this section, we show the time-frequency seismic patterns of strong earthquakes (M ≥ 7) from 1900 to 2021 in the following four major seismic zones: 1) the United States and Mexico, 2) South America, 3) Japan, and 4) Southern China and Northern India, using the wavelet transform. After finding seismic patterns in each of these seismic zones, the oscillation with a given periodicity that groups the historical earthquakes (M ≥ 7) into high and low seismicity will be obtained using the inverse wavelet transform. Each of these oscillations is used to create a probabilistic long-term earthquake prediction model for each seismic zone analyzed using the Bayesian Machine Learning method. Figure 2 shows the wavelet analysis of the strong earthquake records (M ≥ 7) for the southwestern United States and northern Mexico (see Figure 4) between 1900 and 2021. The top panel shows that only seven strong earthquakes have been recorded in this century-long interval and they are heterogeneously distributed between 1900 and 2021. The global wavelet spectrum (left panel) shows periodicities (seismic patterns) at 1.2 ± 0.5, 2.4 ± 0.9, 9.2 ± 2, 15.7 ± 4, and 55 ± 10 years. The time evolution of the power spectral density (PSD) for these periodicities is illustrated in the central panel.

Southwestern United States and Mexico
Each of the periodicities shown in the global wavelet spectrum (seismic patterns) obtained with the wavelet transform implicitly provides information about the intrinsic properties of the tectonic plates of the southwestern United States, the interaction between these plates, the sources both internal and external that modulate the tectonic movement as well as the dynamics of strong earthquakes. In particular, we are selecting seismic pattern of 55 years that indicates the period of recurrence of strong earthquakes. The other periodicities shown in the global wavelet spectrum will be discussed in other future analysis because in this work the objective is to make a forecast for strong earthquakes. Figure 3A shows the probabilistic earthquake prediction model (blue line/shade). This is a model with a 55-year periodicity, and it can be seen that the historical seismic events of the southwestern United States and northern Mexico (vertical blue bars) could be grouped into three groups (I-III). The fact that historical earthquakes can be grouped would indicate that the activity and event are not mere random process but that they are the result of a complex interaction between the tectonic plates and the internal and external factors that modulate the tectonic movement and trigger a series of earthquakes that occur only in the positive phase of the 55-year periodicity. This pattern holds for our all new results for other seismic zones discussed below.
We want to note that there are historical records of strong earthquakes before 1900 and that these earthquakes also occur in the positive phase of the 55-years oscillation. However, in this paper, we do not show these older historical earthquakes. According to the Bayesian ML model, the next active seismic period (M ≥ 7) could probably start in 2040 ± 5 and end in 2057 ± 5 (cluster IV), and no earthquake would be expected in this seismic zone from 2058 ± 5 to 2077 ± 5. Then another new active period of M ≥ 7 earthquakes would begin again around 2078 (cluster V).
A cursory study of Figure 4A shows that earthquakes could apparently occur anywhere. The PDF of the longitude and latitude shows that the seismic zone has a bimodal and trimodal distribution with maxima at −117°and −93°; and 15°, 25°, and 37°, respectively. However, Figure 4B shows the seismic activity's grouping into different classes (see Methodology). In Figure 4B earthquakes with magnitudes between 5 and 6 are shown in shades of green. Earthquakes with magnitudes between 6 and 7 are shown in all yellow-orange shades. Earthquakes greater than seven are shown in shades of red-brown. If strong earthquakes (red triangles) have a random distribution, then no more than one should occur in the same area. In addition, it can be seen the grouping of strong earthquakes that are in the areas delimited by a black curve, which are practically around the geological faults. This result can be interpreted in terms of two scenarios for the successive strong earthquakes.
The first scenario is that the next high seismic season, including the "Big One" could occur according to our model between 2040 ± 5 and 2057 ± 5 (see Figure 3), around any of the areas outlined with a black line ( Figure 4B). Also, it is possible that occur around the yellow-orange areas, which are the areas FIGURE 4 | Earthquakess in the Southwestern United States and Mexico. (A) Seismic activity between 1900 and 2021 is shown as yellow triangles for 5 ≤ earthquakes magnitudes < 7. Earthquakes magnitudes ≥7 are marked with red triangles. Plate motion vectors (rotated based on the direction of plate movement), plate limits and earthquakes data were taken from the USGS. The PDF (probability density distribution) of the longitude and latitude shows that the seismic zones have a bimodal and trimodal distribution, respectively (see main text for more details and discussion). (B) Distribution of seismic hazards in the Southwestern United States and Mexico. The probabilistic spatial clustering occurrence for 5 ≤ earthquakes magnitudes < 6 is shown in shades of green; for 6 ≤ earthquakes magnitudes < 7 is shown in shades of orange; and for earthquakes magnitudes ≥ 7 is shown in shades of red.
Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 905792 where historically, earthquakes of categories 6 and 7 have occurred.
The second scenario is that successive strong earthquakes can occur arbitrarily and sporadically in any zone. However, the fact that strong earthquakes are clustered temporally and spatially may indicate that there are both temporal and spatial seismic patterns that had not previously been considered for their study and forecast. So from the ML point of view, this second scenario is less likely.

Forecasting Earthquakes at Parkfield, California
The Parkfield section of the San Andreas fault (California, United States) was officially recognized by the United States government as a seismic physics laboratory for developing earthquake forecasts due to its apparent regularity in six moderate earthquakes (magnitudes between 5 and 7) since 1857 (Bakun and Lindh, 1985). The interval between these seismic events at Parkfield on 9 January 1857; 2 February 1881; 3 March 1901; 10 March 1922; 8 June 1934; and 28 June 1966 is, on average, 21-22 years (Bakun and McEvilly, 1984;Bakun and Lindh, 1985;Bakun et al., 2005). Based on the average recurrence time of 22 years, the Parkfield recurrence model by Bakun and Lindh (1985) forecasted that the next following characteristic Parkfield earthquake could have occurred on 1988.0 ± 5.2 years, i.e., the possible next Parkfield earthquake should occur between 1983 and 1993. However, this expected Parkfield earthquake never occurred within the anticipated interval until 2004 (Bakun et al., 2005).
Our methodology proposed for analyzing strong earthquakes can also be applied to study moderate earthquake activity and events at Parkfield. We seek to analyze and explain why the seventh earthquake could never occur in 1988.0 ± 5.2 years as predicted, but instead occurred between 1993 and 2007. In addition, we made the forecast for the eighth earthquake in Parkfield. One of the differences between the methodology proposed by Bakun and Lindh (1985) and the methodology in this work is that while Bakun and Lindh (1985) had chosen an average value between the events, we use the periodicities (seismic patterns) of the Parkfield earthquakes. A mean value is a global characteristic of a phenomenon or event, while a periodicity is an intrinsic property of the phenomenon or events (see, for example, Velasco Herrera et al., 2021; Velasco Herrera et al., 2022a; Velasco Herrera et al., 2022b). Furthermore, a periodicity is intrinsically related to spectral and temporal power, which is a fundamental concept of physics (Landau and Lifshitz, 1988a;Feynman et al., 2011a). In another deeper sense, a periodicity is also related to the underlying symmetry of the physical phenomenon (Wigner, 1967).
The mean value has the characteristic that all events associated with this value oscillate around it. In fact, the Parkfield seismic events vary between 12 and 38 years, so the mean value of these seismic events is 21 years. Therefore, using the mean value of seismic events to forecast earthquakes from the point of view of signal theory, signal processing, and machine learning is not the most appropriate. While, the periodicities in the wavelet spectra of earthquakes allow us to identify the intrinsic properties of the earthquakes, determine the characteristics of the earthquake's source, and ultimately determine the interaction between the earthquake and its source and/or the interaction between the earthquake and the faults involved (see for example, Ramírez-Rojas et al., 2019;Soon et al., 2019). In addition, the periodicities cluster the events in high and null seasons (Velasco Herrera V. M. et al., 2022), which allows for constructing models for the forecasts of events with Machine Learning (Velasco Herrera et al., 2021;Velasco Herrera et al., 2022a;Velasco Herrera et al., 2022b). Figure 3B shows the grouping of earthquakes in Parkfield with the period of 22 years proposed by Bakun and Lindh (1985) for the first six seismic events between 1857 and 1966. For the first five events, the group of historical earthquakes are in the positive phase of the 22-years oscillation grouped by the clusters I to V. Also, these five events occur practically during the five maxima of this periodicity. However, in cluster VI, there was no Parkfield earthquake, but the sixth characteristic earthquake in Parkfield occurred in the positive phase of Cluster VII, which was the 1966 event.
With these six seismic events in Parkfield, we can offer a prognosis for the next characteristic earthquakes of Parkfield with the 22-year period proposed by Bakun and Lindh (1985). The result of this forecast is clustering labeled VII to X. According to the periodicity of 22 years, the seventh characteristic earthquake must have occurred in the positive phase of cluster VII which was between 1979 and 1990. But, the seventh event occurred in 2004 (this seventh characteristic event seismic is not shown in or directly indicated on Figure 3B because it was never used for the forecast) i.e., during the positive phase of cluster IX. Therefore, earthquakes did not occur in Parkfield in two clusters (VI and VIII). Based on these results, the empirical evidence suggests that the periodicity of 22 is not a characteristic seismic pattern of earthquakes around the Parkfield seismic zone. This is because, according to the grouping of the 22-years oscillation, the forecast either can be fulfilled or cannot be fulfilled. The fact that an earthquake did not occurred (as predicted) could mean that there were no human and economic losses for the population that lives near these active seismic areas. However, for the development of the science of earthquake forecasting, the fact that a predicted earthquake or earthquake event did not occurred means that the proposed model does not have all the necessary information. In turn, such failure simply means that it is necessary for a more serious and careful re-analysis of the data and a deeper understanding of the particular seismic zone.
Almost 2 decades after the 2004 Parkfield earthquake and nearly 4 decades after the pioneering Bakun and Lindh (1985) forecast, it is necessary to explain why physically and geologically, an earthquake could not occur in 1988.0 ± 5.2 years as proposed. In addition, it is necessary to correct the pioneering/original model proposed by Bakun and Lindh (1985). Wavelet analysis (figure not shown) of the earthquake activity around the Parkfield seismic section show periodicities of 6.2, 11.1, 20.8 and 35.1 years. We are not going to focus on the periodicities of 6.2 and 11.1 years in this paper. We note that the periodicity of 20.8 ± 5 years, with its associated uncertainty, is practically the 22-years oscillation proposed by Bakun and Lindh (1985), which was presented in Figure 3B.
Next we will focus on the result of the 35.1± 7-year periodicity from our wavelet analysis. Figure 3C shows the groupings of historical earthquakes with a periodicity of 35.1 ± 7 years. It can be observed that all, absolutely all, Parkfield earthquakes are occuring during the positive phase of the six clusters identified. This is the first contrast of our results with the model proposed by Bakun and Lindh (1985). Furthermore, it is recognized that more than one event can occur within a single cluster. Such is the case of cluster IV, where two seismic events occurred.
According to the Bayesian Machine Learning model with the chosen period of 35 years, the characteristic earthquakes in Parkfield can never occur in the negative phase, especially between 1978 and 1992, as proposed by Bakun and Lindh (1985). Our model suggests that during cluster VI (which is between 1994 and 2007), the seventh seismic event would occur. Again, in this forecast, we do not use the 2004 earthquake to make the actual forecast; which is the reason why we did not plotted the 2004 event in Figure 3C. Indeed, we wish to highlight in this paper that the seventh event at Parkfield is correctly "forecasted" within the cluster VI. In addition, we want to note that the next season in which at least one characteristic Parkfield earthquake can occur would be around cluster VII starting in 2019 and ending in 2032. Except for cluster IV, Parkfield seismic events occur around the maximum of each cluster. So with a high probability, the eighth earthquake in Parkfield seismic zone/ section can be reasonably expected to occur at around 2025 ± 2 years. This date will be coming soon, so it will be possible to follow in great detail the seismic events in this well-recognized recognized seismic laboratory at Parkfield, California. We are cautiously hopeful that this is indeed a significant area for the development of the science of earthquake forecasts in order to proffer effective and reliable early warnings with the worthy objective of minimizing human and economic losses. Figure 5 shows the wavelet analysis of the earthquake records (M ≥ 7) for southwestern Mexico (see Figure 4) between the years 1900 and 2021. A more significant number of earthquake is observed here when compared to the southwestern United States and northern Mexico. This shows a completely different seismic activity and pattern between Northern and Southern Mexico. Also, this could indicate that in southwestern Mexico there is less viscosity between the tectonic plates when compared to the southwestern United States and northern Mexico. The global wavelet spectrum (left panel) shows periodicities at 1.3 ± 0.4, 2.2 ± 0.5, 3.8 ± 0.6, 7.7 ± 1, 16.4 ± 3, 30.1 ± 4, and 56 ± 5 years. The time evolution of the power spectral density (PSD) for these periodicities is shown in the central panel. We have selected the periodicity of 3.7 years to group the strong earthquakes in southwestern Mexico. Figure 6 shows the probabilistic earthquake forecast (M ≥ 7) for southwestern Mexico. This model has adopted a nominal recurrent periodicity of 3.7 years and it can be seen that this model groups the historical earthquakes (black bars) into nineteen clusters. Furthermore, it can be seen that these seismic events also occur in the positive phase of the 3.7-years oscillation. This characteristic shows that strong earthquakes in southwestern Mexico are not temporally random. According to this model, the next period of greater earthquakes (M ≥ 7) would start in 2024 ± 1 and last until 2026 ± 1. This is clearly a testable scientific proposition.

Mexican Earthquake
In Figure 4B, the cluster of strong earthquakes in southwestern Mexico can be studied. It can be seen that these earthquakes occur preferentially in the subduction zone. In addition, the following strong earthquakes may occur in the cluster zones and in the yellow-orange areas where magnitude 6 to 7 earthquakes have historically occurred. In addition to the interplate earthquakes in southwestern Mexico, there are also intraplate earthquakes in central Mexico. If successive earthquakes occur in this area, it could cause significant damages in the cities of central Mexico with serious human and economic losses (e.g., Novelo-Casanova et al., 2013). Figure 7 shows the wavelet analysis for earthquakes (M ≥ 7) in the South American seismic zone. The global wavelet spectrum shows the periodicities of 1.10.3, 2, 2 ± 0.5, 3.6 ± 0.7, 4.5 ± 0.7, 7.7 ± 1.7, 12.1 ± 2.5, 24.6 ± 5.5, and 46.8 ± 8.4 years.

South America
In order to group strong earthquakes, we use the periodicity of 7.7 years. This nominal choice of longer recurrent periodicity may indicate or imply that there is greater viscosity between the tectonic plates in South America than in southwestern Mexico. Figure 8 shows the probabilistic earthquake forecast (blue line/ shade) model of 7.7-years that grouped historical earthquakes into fifteen clusters. According to this model, the next seismically active period would begin in 2026 ± 2 and end in 2031 ± 2 (cluster XVI).
In Figure 9B, the spatial clustering of strong earthquakes in South America is shown. It can be seen that these earthquakes occur preferentially in the subduction zone. In addition, the following strong earthquakes may occur in the cluster zones and yellow-orange areas where magnitude 6 to 7 earthquakes have historically occurred.

Japan
Owing to the significant and relatively more frequent seismic activity in Japan, we have divided earthquakes in the Japanese zone into two groups. The first group consists of earthquakes 7 ≤ M < 8. After the 11 March 2011s "Great East Japan Earthquake", offshore of the Tohoku region (see, for example, Davis et al., 2012, for a full discussion about the missed opportunity for disaster preparedness), there is a new fundamental question in developing earthquake forecasts in Japan: When could another similar or equal magnitude earthquake occur? To offer a possible answer, we analyze earthquakes in Japan for magnitudes equal to or greater than 8. Therefore our study of the second group is focused on the strongest earthquakes M ≥ 8.  7). Bayesian inference of LS-SVM model (blue line amd shade) compared with the historical earthquakes clustered in fifteen groups (I-XV) and shown with vertical blue bars from 1900 to 2021. In addition, the probabilistic earthquake prediction is shown for the following two periods (cluster XVI and XVII). The blue shaded area represents the 95% confidence intervals of the Bayesian model. Figure 10A shows the wavelet analysis for the Japanese earthquakes of the first group. The global wavelet spectrum shows the periodicities of 1.2 ± 0.3, 2.4 ± 0.5, 4.1 ± 0.7, 10.9 ± 1.5, and 36.8 ± 5 years. Figure 10B shows the wavelet analysis for the second group of Japanese earthquakes M ≥ 8. The global wavelet spectrum shows the periodicities of 1.1 ± 0.3, 1.7 ± 0.5, 5.1 ± 0.7, 23.1 ± 4.5, and 40 ± 7 years. Figure 11A shows the probabilistic earthquake forecast (7 ≤ M < 8) model (blue line/shade) of 4.1-year that grouped historical earthquakes into twenty-nine clusters. We want to highlight a high seismic activity between clusters IX and X, XIV and XV, and XVI and XVII. Furthermore, there was no seismic activity in cluster XXIV. Seismic cluster number XXIX began in 2019 ± 1 and will end in 2022 ± 1. Therefore, it is possible that strong earthquakes may occur during 2022. During the preparation of this work, on 16 March 2022 an earthquake of magnitude 7.3 (37.702 o N, 141.587 o E) was recorded in Japan. So this strong earthquake verifies the accuracy of the Bayesian machine learning event classification and forecast for Japan in real time. According to this model, the next active seismic period would begin in 2024 ± 2 and end in 2029 ± 2 (cluster XXX). Here is another imminently testable forecast contributed by our Bayesian ML algorithm and analyses. Figure 11B shows the probabilistic forecast of earthquakes M ≥ 8. This model has a 40-years oscillation and groups historical earthquakes into three clusters. It is observed that earthquakes occur in the positive phase and that the next period of high seismic activity would begin in 2035 ± 5 and end in 2050 ± 5 (cluster IV).
In Figure 12A shown the seismic activity in Japan between 1900 and 2021. The PDF of the longitude and latitude shows that the seismic zone has a trimodal and bimodal distribution with maxima at 131°, 142°, and 147°; and 37°and 43°, respectively. In Figure 12B, the spatial clustering of strong earthquakes in Japan is shown. It can be seen that these earthquakes occur preferentially in the subduction zone. In addition, the following strong earthquakes may occur in the cluster zones and yelloworange areas where magnitude 6 to 7 earthquakes have historically occurred.
We would like to highlight that the earthquake of 16 March 2022 that was recorded in Japan with a magnitude of 7.3 (with epicenter at 37.702°N, 141.587°E) which was recorded off the coast of Japan's Fukushima prefecture. It is within the areas where, according to our spatial model shown in Figures 12A,B strong earthquake would be expected. Figure 13 shows the result of wavelet analysis for strong earthquakes with M ≥ 7 around the Southern China-Northern India zone. The global wavelet spectrum shows the periodicities of 1 ± 0.3, 1.8 ± 0.5, 3.4 ± 0.7, 6.9 ± 1.2, 8.6 ± 1.3, 13 ± 2, and 20.7 ± 5 years. Figure 14 shows the probabilistic forecast of the Southern China-Northern India earthquakes with M ≥ 7. This model has a 8.6-years oscillation and groups historical earthquakes into twelve clusters. It is observed that earthquakes occur in the positive phase and that the next period of high seismic activity would begin in 2022 ± 1 and end in 2028 ± 2 (cluster XIV).

Southern China-Northern India
In Figure 15A shown the seismic activity in Southern China-Northern India between 1900 and 2021. The PDF of the seismic activity longitude and latitude shows that the seismic zone has bimodal distributions with maxima at 77°and 94°; and 30°and 40°, respectively. In Figure 15B, the spatial clustering of strong earthquakes in Southern China-Northern India is shown. It can be seen that these are intraplate earthquakes. In addition, the following strong earthquakes may occur in the cluster zones and yellow-orange areas where magnitude 6 to 7 earthquakes have historically occurred and we would like to highlight that strong seismic activities are preferentially distributed around fault lines.
FIGURE 11 | (A) Probabilistic earthquake prediction for Japan with magnitudes 7 ≤ M < 8. Bayesian inference of LS-SVM model (blue line and shade) compared with the historical Japanese earthquakes clustered in twenty nine groups (I-XXIX) and shown with vertical blue bars from 1900 to 2021. In addition, the probabilistic earthquake prediction is shown for the following three periods (cluster XXX-XXXII). (B) Probabilistic earthquake prediction for Japan with magnitude greater than or equal to 8. Bayesian inference of LS-SVM model (blue line and shade) compared with the historical Japanese earthquakes clustered in three groups (I-III) and shown with vertical blue bars from 1900 to 2021. In addition, the probabilistic earthquake prediction is shown for the following period (cluster IV). The blue shaded area represents the 95% confidence intervals of the Bayesian model.

DISCUSSION AND CONCLUDING REMARKS
Mechanisms conducting the plate motions and the Earth's geodynamics related to the triggering, the persistency and the potency of earthquakes are still not entirely clarified (Kanamori and Brodsky, 2004;Doglioni and Panza, 2015;Senapati et al., 2022, and several references cited in these papers) Although it is not the subject of diagnosis in the present work, the episodic elastic accumulation and release of energy produced by earthquakes depend on the multiple prevailing characteristics when the shear stress exceeds the internal friction between the sliding plates.
The relationships between stress fields and frictional responses (Lockner and Beeler, 2002) result from the rheological behavior of the materials that influence the viscosity of the affected plate contacts. The pressure by crust overloading and temperature conditions varies considerably from depth, affecting in consequence the origin, intensity and frequency of earthquakes. The velocity and the geometry of spatial orientations of the tectonic convergences also interact both in planar sense (magnitudes of the transcurrent components) and in depth (subduction dip angles).
Other complementary surficial factors such as tidal friction of the Earth could also induce the earthquake's triggering (Zschau, 1986;Wilcock, 2001) due to the accumulation of energy for centuries (Doglioni and Panza, 2015). Based on these multiple geophysical factors, we consider the main causes of the greater recurrency and intensity of earthquakes that are expected in longlived tectonic contexts of near to orthogonal convergence and deeper subduction inclinations. Kossobokov (2004) suggested that an apparent irregularity and a certain infrequency of earthquake occurrences may ultimately hinted that earthquakes are ultimately unpredictable phenomena rooted in stochastic processes beyond any deterministic rules or probabilities. But the appearance of stochastic nature or process may be due to the fact that abrupt or sporadic events such as earthquakes, and explosions, among others, must be analyzed in a different way than gradual processes such as temperature, atmospheric pressure, and other physical phenomena. Velasco Herrera et al. (2017) suggested that there is another type of natural phenomena that occurs only in a specific phase of a very particular oscillation such as strong earthquakes. The apparent irregularity of the strong or moderate earthquakes analyzed in this work could be observed in the heterogeneous distribution shown by the seismic events in the positive phase of each cluster obtained with Machine Learning.
The grouping of historical earthquakes makes it possible to find the periods of high and null seismic activity, but at the moment, it cannot answer the deeper questions posed by Kossobokov (2004): Why, where and when do earthquakes occur? In addition, Kossobokov (2004) raises a question that, FIGURE 12 | Earthquake zones in Japan. (A) Seismic activity between 1900 and 2021 is shown as yellow triangles for 5 ≤ earthquakes magnitudes < 7, for 7 ≤ earthquakes magnitudes are marked with red triangles. Plate motion vectors (rotated based on the direction of plate movement), plate limits and earthquakes data were taken from the USGS. The PDF (probability density distribution) of the longitude and latitude shows that the seismic zone has a trimodal and bimodal distribution, respectively. (B) Distribution of seismic hazards in Japan. The spatial probability of occurrence for 5 ≤ earthquakes magnitudes < 6 is shown in shades of green; for 6 ≤ earthquakes magnitudes < 7 is shown in shades of orange; and for earthquakes magnitudes ≥ 7 is shown in shades of red.  to date, still does not have an answer. Are earthquakes predictable? Although a negative answer, according to Kossobokov (2004), is merely a guess. In addition, Kossobokov (2004) asks if there are intrinsic temporal and physical characteristics of earthquakes that can be used in other forecasting types to can be used to reduce human and economic losses? Different algorithms and models, 4,5 have been used to forecast earthquakes in different seismic zones (see, for example, Michael and Werner, 2018;Schorlemmer et al., 2018, for a full review). The medium-term prediction has been carried with the M8 algorithm for strong earthquakes (Keilis-Borok and Kossobokov, 1990). In addition, the M8 algorithm demonstrated that the largest earthquakes are not a random process (Kossobokov and Soloviev, 2021). Also, the strong earthquakes are forecasted with a limited precision in their ranges of time, space, and magnitude (Kossobokov and Soloviev, 2021). We propose a probabilistic algorithm for forecasting earthquakes (strong or moderate) that can be applied and implemented in different seismic zones. Different algorithms and methodologies for forecasting earthquakes require different tests. For example, we may note that the tests for the M8 algorithm by V. Kossobokov and colleagues(see, Kossobokov and Soloviev, 2021) and for our Bayesian Machine Learning models are very different (see the Methodology Section 3.6.2).
However, the proposed deterministic methodologies have not yield effective results nor robust conclusions. The study of seismic precursors has increased in recent years, but there is no reliable method to predict earthquakes over long time horizons/windows of multi-years or even decade to multidecades (Pulinets and Boyarchuk, 2005;Ouzounov et al., 2018;Pulinets and Ouzounov, 2018).
Predicting earthquakes must be one of the most significant challenges in modern science and technology. Earthquake prediction is necessary to minimize the enormous earthquake hazard risks in a seismic area. In particular the prediction of earthquakes is essential to minimize human tragedies and economic losses. The occurrences of an earthquake are complex and largely non-linear in character, which is why there is no deterministic model that can predict the exact location, magnitude, and time of an earthquake of any significant magnitudes.
There is currently a great debate in the scientific community about the origin of earthquakes. While some consider that it may not possible to predict earthquakes (Geller et al., 1997), others like us suggest that it could be a predictable phenomenon. Our point of view is informed by the seismic patterns found in the seismic zones analyzed adopting the methodology applied in this work.
Furthermore, we suggest that temporal forecasting of strong earthquakes should consist of forecasting magnitude widths/ FIGURE 15 | Earthquakes in Southern China-Northern India. (A) Seismic activity between 1900 and 2021 is shown as yellow triangles for 5 ≤ earthquakes magnitudes < 7. Earthquakes magnitudes ≥7 are marked with red triangles. The PDF (probability density distribution) of the longitude and latitude shows that the seismic zone has both a bimodal distribution, respectively. Plate motion vectors (rotated based on the direction of plate movement), plate limits and earthquakes data were taken from the USGS. (B) Distribution of seismic hazards in Southern China-Northern India. The spatial probability of occurrence for 5 ≤ earthquakes magnitudes < 6 is shown in shades of green; for 6 ≤ earthquakes magnitudes < 7 is shown in shades of orange; and for earthquakes magnitudes ≥ 7 is shown in shades of red. ranges rather than forecasting any "exact" magnitude. Moreover, from the geological point of view, the exact spatial forecast is not a correct concept since the geological faults are not point-like entities. In addition, the release of accumulated energy does not occur at a single point, nor is it instantaneous. Therefore, the energy released from an earthquake is an average of the interrelated cascade processes of rupture of the geological fault in a volumetric geographical area that occurs within a time interval. In other words, an earthquake does not occur at one point, nor is it a momentary event. We suggest that the problem of earthquake predictions should be viewed as a question of probability. One possible solution to earthquake forecasting is to calculate the predicted probabilities of future seismic cycles.
In our point of view, the key challenge in prediction of seismic activity today is to shift paradigms from deterministic forecasts to a probabilistic approach to predicting earthquakes with reliable estimates or quantifications of uncertainties. So earthquake prediction should be a multidisciplinary task that accounts for the recent advances in artificial intelligence and should be widely applied for earthquake forecasting (Beroza et al., 2021).
We have studied seismicity in North America, South America, Japan and Southern China and Northern India with Machine Learning by analyzing seismic patterns of variations for earthquakes with magnitude 7 or greater between 1900 and 2021. We then created a probabilistic earthquake prediction model for each seismic zone analyzed using the Bayesian Machine Learning method. Each model obtained groups the seismic of magnitude greater than 7 in clusters. Our result also partially explains the periods of earthquakes of magnitude 7 and the apparent seismic tranquillity for earthquakes of magnitude 7 or greater. We want to clarify that the periods where no earthquakes of magnitude 7 do not imply the total absence of seismic activity, indeed earthquakes of lesser magnitude must still occur during this period.
We suggest that if the dynamics of the tectonic plates have not changed in the last thousands of years, then no abrupt change in the tectonic movement would be expected, and therefore, the seismic pattern found should remain stationary, at least temporarily and spatially stable and unchanging, for the following next few decades.
We speculate that the seismic patterns found, that is, the periods of occurrence of earthquakes in each seismic zone analyzed could be interpreted as the period in which the energy accumulates and this energy releases through the rupture along faults and fractures near the plate tectonic boundaries.
As it can be seen in the results we obtained, each one of the analyzed areas has different characteristic frequencies for earthquakes of the analyzed magnitudes. This is probably due to the type of plate boundary, the lithological composition of the related plates and their consequent different friction coefficients, the geometry of the margin, the convergence angle of the plates, as well as the velocity of the approach. We leave this briefly expressed hypothesis awaiting further future studies.
As a highlight of our results, we note that the larger number of clusters are found and defined for the Himalayan collisional margin. We could speculate that this phenomenon is related to the diverse lithological nature of the constituent materials of the two continental crustal plates that come into contact in a frontal collision with very different friction coefficients. The cases of the subduction convergence zones of Southwestern Mexico and South America seem to have similar frequencies, sensu lato, a fact that is not surprising because of the type of crust of the plates involved in the boundaries (oceanic crust versus continental crust).
To the contrary, convergence by subduction of two oceanic crust plates, as in the case of Japan, presents a notably higher frequency. This would be related to the physical properties of oceanic basalts and their response to imposed stresses.
In the Southwestern North American transform margin, we observed a very low frequency oscillation and modulation. In this case, the low frequency of earthquakes of magnitude greater than 7 could be related to the rectilinear geometry of the plate contact (plus its transcurrence) and the relation of this with the accumulation-release of elastic energy.
In addition, we would like to highlight the forecast for the earthquakes in the Southwestern United States of America since one of the biggest concerns is any major earthquake in the densely populated region of San Andres fault, such as the catastrophic event on 18 April 1906. This event is known as the San Francisco great earthquake. So after more than a hundred years, the "Big One" could occur according to our model between 2040 ± 5 and 2057 ± 5. Although this great earthquake may occur "tomorrow", we still have a little time to refine our forecasts of such strong earthquakes.
In summary, we have analyzed seismicity in North America, South America, Japan, and Southern China-Northern India for M ≥ 7 earthquakes from 1900 to 2021. The primary seismic patterns found for M ≥ 7 earthquakes in the analyzed seismic zones are 55, 3.7, 7.7, and 8.6 years, respectively, for the southwestern United States and northern Mexico, southwestern Mexico, South American, and Southern China-Northern India. In the Japanese zones, the primary seismic pattern for 7 ≤ M < 8 earthquake is 4.1 and 40 years for M ≥ 8 earthquakes.
Our Machine Learning models show that there are periods where there are earthquakes magnitude ≥7 and periods without earthquakes with magnitude ≥7 in the analyzed seismic zones. In addition, our Machine Learning models predict a new seismically active phase for earthquakes magnitude ≥7 between 2040± 5and 2057 ± 5, 2024 ± 1 and 2026 ± 1, 2026 ± 2 and 2031 ± 2, 2024 ± 2 and 2029 ± 2, and 2022 ± 1 and 2028 ± 2 for the five seismic zones in United States, Mexico, South America, Japan, and Southern China-Northern India, respectively. Finally, we note that our algorithms can be further applied to perform probabilistic forecasts in any seismic zone.
Our algorithm for analyzing strong earthquakes in extensive seismic areas can also be applied to smaller or specific seismic zones where moderate historical earthquakes with magnitudes between 5 and 7 occur, as is the case of the Parkfield section of the San Andreas fault (California, United States). Our analysis shows why a moderate earthquake could never occur in 1988 ± 5 as proposed by Bakun and Lindh (1985) and why the long-awaited characteristic Parkfield earthquake occurred in 2004. Furthermore, our Bayesian model of Machine Learning adopting a periodicity of 35 years predicts that possible seismic events may occur between 2019 and 2031, with a high probability of event(s) around 2025 ± 2. The Parkfield section of the San Andreas fault is an excellent seismic laboratory for developing, testing, and demonstrating earthquake forecasts. In a few years, it will be possible to demonstrate whether our algorithm effectively forecasts strong and moderate earthquakes. We may note in anticipation that if some of our forecasts are not fulfilled in some of the analyzed seismic zones, it may not be the sole fault of Machine Learning algorithms. Instead the complex issues and problems may lie with our conjectures and the proposed models for each seismic zones, so we will have to carefully re-analyze again the seismic zone or zones where the forecast was not fulfilled. If all our forecasts for the next high season of earthquakes are fulfilled, then we must incorporate more elements such as the seismic precursors (see, for example, Pulinets and Boyarchuk, 2005;Ouzounov et al., 2018;Pulinets and Ouzounov, 2018, for a full review) of each zone analyzed in our models to give more accurate earthquake forecasts in order to provide earlier warnings and greater security to people living in these earthquake zones.
Spatial forecasting models Zechar and Jordan (2008) have been suggested, involving variables and metrics such as 1) Relative Intensity (RI) alarm function, 2) Pattern Informatics (PI), and 3) the United States Geological Survey National Seismic Hazard Map (NSHM). These models represent different assumptions about the spatial distribution of earthquakes. Concerning RI: the hypothesis is that future earthquakes are more likely to occur where seismicity is historically higher. On PI: the hypothesis is based on the fact that anomalous seismic activities indicate the locations of future earthquakes. Finally, on NSHM: the hypothesis suggests that future earthquakes will occur where previous earthquakes have occurred, and that some earthquakes may occur anywhere. In our case, the possible zones where the following strong earthquakes in each seismic zone analyzed could occur have been essentially clustered and pre-determined according to the methodology described in Section 3.1. According to this spatial grouping, the fact that strong earthquakes occur probabilistically close to where strong earthquakes have historically occurred could indicate certain information about the tectonic plates. For example, precisely in those areas, there is greater fracturing compared to areas where strong earthquakes have not occurred. In addition, the probabilistic spatial clustering of Figures 4, 9, 12, 15 could show the areas with the highest probability of strong earthquakes. Therefore, it would be essential to analyze these highly probable areas of strong earthquakes with all the models and precursors that are currently known in order to minimize economic losses and human losses.
In conclusion, we believe that our results demonstrate that our methodology is a good alternative to traditional deterministic earthquake prediction. Thus, the problem of earthquake predictions should be considered as a question of probability. We believe the challenge in the study of seismic activity is to modify the forecast paradigm to a probabilistic earthquake prediction.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.