Electromagnetic and Radon Earthquake Precursors

Abstract

Earthquake forecasting is arguably one of the most challenging tasks in Earth sciences owing to the high complexity of the earthquake process. Over the past 40 years, there has been a plethora of work on finding credible, consistent and accurate earthquake precursors. This paper is a cumulative survey on earthquake precursor research, arranged into two broad categories: electromagnetic precursors and radon precursors. In the first category, methods related to measuring electromagnetic radiation in a wide frequency range, i.e., from a few Hz to several MHz, are presented. Precursors based on optical and radar imaging acquired by spaceborne sensors are also considered, in the broad sense, as electromagnetic. In the second category, concentration measurements of radon gas found in soil and air, or even in ground water after being dissolved, form the basis of radon activity precursors. Well-established mathematical techniques for analysing data derived from electromagnetic radiation and radon concentration measurements are also described with an emphasis on fractal methods. Finally, physical models of earthquake generation and propagation aiming at interpreting the foundation of the aforementioned seismic precursors, are investigated.

1. Introduction

Earthquakes are inevitable disastrous phenomena. Not only are they unavoidable, but the incredible difficulty in forecasting them renders these disasters even more hazardous and catastrophic. Finding an accurate seismic precursor is one of the greatest challenges for the scientific community worldwide. Seismic forecasting research dates back more than fifty years and, arguably, earthquake forecast still lacks in terms of accuracy and reliability. There is evidence that pre-seismic electromagnetic radiation and radon concentration observations can be utilised in earthquake prediction, taking into account specific measurable features of the associated earthquake process. More specifically, if such observations are made near the geological rupture, some measurable precursory activity prior to the seismic event can be expected. However, the reader should note that it is difficult to establish any direct correlation between prior observations and earthquake occurrences, especially in short-term forecasting [1,2].

Reducing the uncertainty in the estimation of the occurrence time and location or even the size of a forthcoming massive seismic event is the main goal of earthquake forecasting [3]. Seismic forecasting usually falls into four categories [4]: long term (10 years); intermediate term (1 year); short term (${10}^{-1}$ to ${10}^{-2}$ years); and immediate term (${10}^{-3}$ years or less). Hayakawa and Hobara [5] classify earthquake forecasting into three categories: long term (time scale of 10 to 100 years); intermediate term (time scale of 1 to 10 years) and short term (from several weeks to a few hours). Each forecast category relies on different data sources (e.g., instrument measurements for short-term forecasting, historical records and geological studies for long-term forecasting) and is dictated by the public need for earthquake preparedness, which includes an array of safety procedures for each level of forecast [1].

In seismic-prone countries, short-term early warning in a time window of weeks, days, or hours is deemed most important but is inherently more difficult than long-term forecasting. Short-term earthquake forecasting involves the study of short-term precursory activity based on observations of physical quantity irregularities acquired near and before earthquake events and can be further supported by serendipitous findings in observations not purposed for earthquake monitoring but that are nonetheless acquired near the earthquake location [6]. Abnormalities in electromagnetic fields, anomalous variations in radon concentration in the soil, groundwater, surface water and atmosphere, erratic gas emissions, uneven surface distortions caused by pressure differentials, ionospheric perturbations and anomalies detected by satellite and remote sensory devices are among these physical quantity irregularities [6].

Observations of pre-seismic electromagnetic disturbances (in the radiofrequency–RF range) are one of the most promising tools for short-term earthquake forecasting. The related subject is termed seismo-electromagnetism [7]. As it has been shown by many studies (see, e.g., the reviews [5,6,8,9,10,11] and the references therein), pre-seismic electromagnetic emissions occur in a wide frequency range, i.e., for frequencies below 10 Hz (Ultra-Low Frequencies—ULFs), frequencies in the kHz range and up to several MHz (High Frequencies—HFs) and frequencies between 100 MHz and 300 MHz (Very High Frequencies—VHFs). The research originated back in the 1970s, when the first successful seismic forecast was reported for an earthquake of magnitude M = 2.6, occurring on 3 August 1973, near Blue Mountain Lake, New York [12]. Following this, the M = 7.4 Heicheng earthquake of 4 February 1975 was correctly anticipated by seismologists, boosting the prospect that credible earthquake forecasting may be feasible. This forecast led to the issuance of a warning within a period of 24 h before the primary shock, perhaps avoiding more casualties than the 1328 deaths that the event resulted in. A major setback to the earthquake forecast endeavour was the 1976 M = 7.8 Tangshan earthquake, which struck 18 months later and was not anticipated. The number of deaths caused by this earthquake reached the hundreds of thousands [6,8]. Research has recently focused on short-term forecasting rather than long-term forecasting [13]. Pre-seismic electromagnetic observations and abnormality recordings have been documented by several study teams throughout the globe as precursors of earthquakes. Electromagnetic variations are recorded by ground stations, remote sensory devices [14,15] and satellites [14,16].

Radon precursors of pre-seismic activity are also intriguing. Due to its importance, research on radon monitoring has become a rapidly growing topic in the search for warning signs before earthquakes [2,6,8,17,18,19,20,21,22,23,24,25]. This is due to the fact that radon may travel great distances from the host emitting rocks [26] and can be detected at very low levels [27]. Anomalous radon concentration variations in the soil, groundwater and atmosphere may be observed prior to earthquakes [6,17,19,24,28,29]. Before earthquakes, anomalous radon fluctuations are addressed in the soil, groundwater, atmosphere and thermal spas [6,17,19,24,28,29]. The time series features, such as the range, length, number of radon anomalies, precursory time and epicentral distance, vary greatly [6,17,30,31]. However, radon emissions are influenced by seasonal variations, rainfall and barometric pressure alterations [6,8,21,23,26,30,31] and for this reason, radon time series are preprocessed accordingly [6,17,21,22,23]. The majority of the associations between radon and earthquakes involves events of small and intermediate magnitudes. Large magnitudes earthquake associations with radon observations also exist [2,32,33,34,35,36].

Ionospheric studies, satellite measurements and remote sensing measurements have recently gained significant interest in earthquake precursor research following the publication of the Lithosphere–Atmosphere–Ionosphere Coupling (LAIC) [37]. Due to the widespread availability of GPS data, many studies report GPS-based total electron content (TEC) data of the ionosphere, providing valuable information and convincing evidence of precursory activity [38,39,40,41,42]. Other researchers have studied the lower ionosphere extensively in relation to earthquakes by exploring possible ionospheric perturbations in the upper and lower regions [43,44,45,46,47].

Despite continuous research efforts, the understanding of the fracturing mechanisms of the crust remains limited [1,6,8,18,19,20,21,22,48,49,50,51,52,53,54,55,56]. Given that the fracture of heterogeneous materials has still not adequately been explained, despite significant effort at the theoretical and experimental level [1], the interpretation of earthquake generation is still problematic. In addition, each earthquake is unique with its own spatial and temporal distribution and any candidate precursor cannot be generalised to other earthquakes in a straightforward manner. Eftaxias et al. [1] stated that before the final catastrophe, several seismic, geochemical, hydrological and electromagnetic changes occur as the earthquake’s preparatory process unfolds. This suggests that a credible earthquake precursor should be derived on the basis of an acceptable physical model which adequately explains the earthquake process [6] and especially its preparatory stage, and this further complicates the search for earthquake precursors.

2. Electromagnetic Precursors

2.1. ULF Emissions

In 1964, seismogenic electromagnetic emissions with frequencies lower than 10 Hz were first observed [57]. It has been found that ground electric potential variations, atmospheric ULF electromagnetic waves and other related phenomena occur prior to earthquakes [8,9,10,58,59,60,61,62]. Monitoring and recording ULF emissions on the ground, arising directly from the lithosphere, is one of the several widely used seismo-electromagnetic methods. This is because ULFs (f = 0.01 Hz–10 Hz) have great skin depth, low attenuation, less contamination and less penetration through the magnetosphere and ionosphere [63] and as a result, ULF waves can travel up to an observation point close to the Earth’s surface with little attenuation [5]. Although most ULF precursors are electric, nowadays, researchers also study magnetic ULF precursors [6,8,9]. Nevertheless, it should be noted that concerns on the existence of a direct relation between ULF signals and impending earthquakes have also been raised [9].

The VAN method (from Varotsos, Alexopoulos, Nomicos) for ULF emissions has a long track record of more than forty years [10,61,62]. The method introduced the concept of Seismic Electric Signals (SES). SES are ULF disturbances of frequencies f < 1 Hz. The most significant physical property of SES is selectivity [8,9,10], which means that SES prefer specific paths and consequently, a ULF station is sensitive to SES from certain seismic areas only, namely from specific focal areas. The map showing these potential areas is called the selectivity map of a station. Due to their selectivity, SES can be detected even from hundred kilometres away from the epicentre. By installing two dipoles of length 50 m, 100 m, 200 m and preferably 1000 m in a cross configuration, magnetotelluric variations are discriminated from anthropogenic disturbances. The VAN method has successfully forecasted seismic events within a precursory window of some days or weeks in both Greece [61,62] and Japan [5,10]. Nowadays, SES (among other signals) are incorporated in the natural time method (see Section 5), which has demonstrated improved forecasting potential on several earthquakes (e.g., [64] and references therein). Nevertheless, the applicability of the VAN method has divided the scientific community into those supporting [10] and those rejecting it [6].

The 1988 Spitak M = 6.9 earthquake [65], the 1993 Guam M = 8.0 earthquake [59], the 1996 Hetian M = 7.1 event [66] and the 1997 Kagoshima M = 6.5 earthquake [67] were successfully forecasted using ULF electromagnetic emissions. Using the cumulative daily sum of the local energy of the earthquakes weighted by the squared distance from the measurement station, a method was suggested by Hattori et al. [68,69] for ULF-based forecasting of earthquakes. Employing this method, Han et al. [70] reported an increased probability of ULF magnetic anomalies 1–2 weeks before medium and strong shallow earthquakes, hence confirming previous findings published by Hattori et al. [69] and emphasising additionally that the perturbations were better associated with stronger and closer earthquakes. According to data gathered from 17 stations in Japan, statistically significant diurnal geomagnetic anomalies were found two months before the $M_w$ = 9.0 2011 Tohoku earthquake. Comparable estimations were published by Han et al. [71] and Xu et al. [72]. Prior to the catastrophic earthquakes that occurred in September 2015 at Coquimbo, Chile, September 2017 at Chiapas, Mexico, and September 2020 at Vrancea, Romania, ground-based stations recorded pre-seismic ULF anomalous geomagnetic disturbances [64,73,74,75,76,77,78,79]. Pre-seismic perturbations in the spectral density ratio between the horizontal and vertical ULF components were reported by Hirano and Hattori [80] and Ouyang et al. [81]. ULF magnetic field emissions were continuously measured in Agra station in India with the help of three-component coil magnetometers with promising forecasting results [63]. ULF geomagnetic data from the Panagjurishte and Surlari stations in Romania were successfully utilised for the forecast of a $M_w$ = 6.4 earthquake which occurred in the coastal zone of Albania on 26 November 2019 [78].

2.2. HF Emissions

In the range between a few kHz to several MHz, a number of HF emission disruptions prior to earthquakes have been reported [1,23,34,35,82,83,84,85,86,87]. As stated by Hayakawa and Hobara [5], there are two methods to detect the seismic precursors, i.e., a direct measurement of the electromagnetic emissions radiated from the hypocentre of earthquakes in the lithosphere, or the indirect detection of anomalous propagation disturbances in the atmosphere and ionosphere caused by transmitter signals already in place. The identification of HF electromagnetic disturbances can aid in determining the source of seismic activity. According to Eftaxias et al. [1,82], the various frequencies of the HF electromagnetic precursors, in conjunction with the detected time lag between events and impending earthquakes, indicate distinct stages and mechanisms of the earthquake preparation processes. It is also believed that cracking and faulting events at different length scales occur prior to an earthquake, as an abrupt mechanical breakdown of the Earth’s heterogeneous crust, thus the complex operations that occur may be the initial source of numerous electromagnetic precursors [1,82,86].

Several publications (e.g., [1,34,82,83,86] and references therein) suggest that the high persistency and organisation in an electromagnetic time sequence points to the development of a positive feedback mechanism regulating the sudden fracto-electromagnetic process that occurs during earthquake preparation. There is increasing evidence that such a feedback mechanism might be a sign of the earthquake fracture process. However, HF electromagnetic precursors have not been fully understood yet. There is still much to be done to comprehensively delineate the HF electromagnetic precursors. Separating two events that happened at different times, like an earthquake and its potential HF electromagnetic precursor, is a challenging task. New methods may provide more data that would enable one to acknowledge the seismogenic source of detected HF electromagnetic abnormalities and connect them to the pivotal phase of the earthquakes production.

Apart from persistency, the strong antipersistent properties of an electromagnetic time series, as well as the change between persistency and antipersistency, are also evidence of an underlying nonlinear feedback of the system initiating the crack-opening process that leads the system to be out of equilibrium ([34,35,88] and references therein). The reader should note that according to Eftaxias et al. [1], the antipersistent behaviour is comparable to that of systems that experience a continuous phase transition at equilibrium. Stationary-like features possibly observed in antipersistent sections of preseismic electromagnetic time series might also be attributed to the heterogeneous part of the fracturing media. According to Contoyiannis et al. [89], Kapiris et al. [84,85] and Eftaxias et al. [90,91], the precursory electromagnetic antipersistent anomalies are associated with a continuous thermal phase transition with strong critical characteristics. Although observing an HF anomaly is a necessary condition for the anticipation of a forthcoming event, it is far from considering it as a sufficient condition, i.e., as a prerequisite for the occurrence of this event [1,88]. Notably, it is important to rule out any potential relationship of these anomalies with magnetic storms, artificial electromagnetic sources or solar flares [1], with the note that the latter may trigger seismicity and impact the earthquake preparation zone [92]. In relation with this, Anagnostopoulos et al. [93] considered that the sun was an agent provoking seismic activity through coronal holes driven by high-speed solar wind streams.

2.3. VHF Emissions

VHF emissions have also been employed in the search for electromagnetic earthquake precursors. According to Pullinets [94], one of the two authors of the LAIC model [37], LAIC describes a system made up of subsystem interactions and a synergy of several processes, one of which is the electromagnetic emissions which may alter the characteristics of the atmosphere and ionosphere. Although some scientists dispute the precursory usability of VHF emissions (e.g., [8] and references therein), recent scientific interest has been stimulated by this subject. For example, Sorokin et al. [95] report a theoretical physical model that explains the over-horizon propagation of pulsed VHF radiation and the origin of the related seismic phenomena in association with the troposphere’s origins, the thermal effects, the associated IR emissions and the modification of the distribution of plasma in D, E and F ionospheric layers. Ouzunov et al. [96] report atmospheric variations in the intensity of broadband wireless signal propagation correlated with pre-earthquake processes. Since 2012, these authors have continued to conduct ground observations in Bulgaria in the VHF band between 1.8 GHz and 3.5 GHz, discovering phenomena related to a signal’s amplification days or hours before the seismic occurrences, with notable examples the M = 5.6 earthquake of 22 May 2012 in Bulgaria, M = 5.1 earthquake of 12 August 2018 in Albania, the M = 4.1 earthquake of 2 August 2018 in Southern Bulgaria and the M = 5.5 earthquake of 28 October 2018 in Romania. A VHF early warning system is utilised among other systems in Mexico [97]. Moriya et al. [98], on the basis of designing a data-collection system, report several anomalous VHF-band radio-wave propagation events prior to earthquakes, with the most significant being the Tokachi-oki earthquake ($M_j$ = 8.0, $M_j$ being a magnitude defined by the Japan Meteorological Agency) on September 26 2003 and the southern Rumoi sub-prefecture earthquake ($M_j$ = 6.1) on 14 December 2004. Devi et al. [99] states that the VHF emissions indicate unusual atmospheric parameters brought on by earthquake precursor processes, which may allow for the reception of VHF communications at distances of more than 1000 km. According to the authors, the lower VHF TV transmissions of less than 70 MHz are linked to modifications in the tropospheric environment and the ionospheric mode of propagation. According to Erickson [100], anthropogenic electromagnetic emissions are primarily caused by mobile communications, car ignition systems, industrial equipment and radio and television broadcasting stations. Eftaxias et al. [101] report VHF disturbances prior to earthquakes in Greece showing that the related features are possibly correlated with the fault model characteristics of the associated earthquake and the degree of geotectonic heterogeneity within the focal zone.

2.4. Remote Sensing and Satellite Techniques

The application of spaceborne remote sensing has grown in popularity and effectiveness within the field of natural disasters [102]. Improved quality data with repeated spatio-temporal coverage over large areas in rough geomorphological and geological conditions can be obtained through the development of geospatial technologies and advanced data processing [103,104,105]. The post-disaster visualisation of remote sensing images helps in knowledge production, emergency intervention thinking and decision-making during earthquakes [106]. Consequently, seismo-electromagnetic research has entered a new phase with the development of remote sensing tools. It is now possible to simultaneously monitor a number of locations throughout the globe, including seismic occurrences taking place in tectonic systems with differing geomagnetic conditions. That is essential to the related research. For the remote sensing data to yield reliable findings, a worldwide coverage with sufficient spatial and temporal resolution is needed [8].

The remote sensing of the co-seismic effects of earthquakes is of importance. Co-seismic effects occur as around 100 m long Earth’s cracks, which are followed by landslides, lateral spreading and changes to urban and suburban areas, which, in turn, may cause human casualties and affect agriculture and the food supply chain. For instance, landslides and substantial lateral spreading were noted following the earthquakes in Kashmir in 2005 and Mirpur in 2019 [107,108]. The seismically generated structural, stratigraphic and hydrogeological side effects are a subject of remote sensing and satellite studies. For example, the Landsat satellite imagery has been used to study the dynamic relationship between observed seismicity and lineament density [109]. Ground-based remote sensing techniques are efficient non-destructive geophysical methods that provide high-resolution subsurface images to detect several co-seismic features. Remote sensing data from a number of satellites and sensors are also useful tools for co-seismic mapping.

Synthetic Aperture Radar (SAR) is among the best remote sensing techniques for mapping the co-seismic changes. Interferometric SAR (InSAR) is one of the most powerful remote sensing technique of the SAR family used to detect several surface deformations over large areas with high accuracy [110]. InSAR allows low-speed surface deformations to be detected over vast areas with centimetre to millimetre precision [111]. The permanent Scatterers InSAR is a robust remote sensing technique for mapping co-seismic deformation and movements of urban infrastructures [112,113] for which the spaceborne remote sensing techniques are less effective because they cannot provide complete near-surface information needed for such estimations.

The Ground Penetrating Radar (GPR) remote sensing technique has been applied to shallow subsurface seismic investigations due to its high-resolution, time and cost-effective nature [114,115,116]. GPR has gained popularity in studies related to the detection of faults and fracture networks [117], slope instabilities [116] and landslides [118]. GPR is a reliably accurate mapping tool to study a single site and a localised subsurface deformation but is difficult to perform over extensive earthquake’s Karst depression–landslide-affected areas. However, only a few studies using field GPR measurements report data regarding the co-seismic liquefaction and the related ground failure [119,120].

2.5. TEC

Total Electron Content, or TEC, is the electron density of a 1 ${\mathrm{m}}^2$ cylinder that is vertically stacked from a ground point to the ionosphere [121,122]. One TECU is the TEC measurement unit and equals ${10}^{16}$ electrons per square meter vertically arranged up to the ionosphere. By definition, TEC is associated with the LAIC model. GPS receivers and ionosondes are used to continuously monitor TEC at various locations across the world [122]. The corresponding data are accessible through a number of repositories and URLs [122,123,124,125,126] via the Ionosphere Exchange (IONEX) data file structure [127].

To investigate seismically generated TEC fluctuations in the ionosphere, researchers have used a variety of schemes and approaches [39,122,128,129,130,131,132,133]. There has been much discussion about the anomalous variations in the ionospheric $F_2$ peak electron density $N_m{F}_2$ (plasma frequency $f_o{F}_2$), which are recorded by ionosondes and TEC, which, in turn, are determined by ground-based GPS receivers and appear prior to earthquakes [39]. Based on 184 $M\ge$ 5.0 earthquakes which occurred in Taiwan over a 6-year period between 1994 and 1999, Liu et al. [39] conducted a statistical investigation that showed anomalous decreases in the ionospheric $N_m{F}_2$ in the afternoon within 1–5 days prior to the earthquakes and pronounced reductions in the ionospheric GPS TEC in the afternoon and late afternoon periods within 5 days prior to 20 $M\ge$ 6.0 earthquakes in Taiwan.

According to Liu et al. [129], while pre-earthquake ionospheric anomalies may occur almost at any local time, TEC over a possible epicentre region typically decreases or increases significantly in the afternoon and/or evening periods, one to six days prior to the occurrence of an earthquake. According to these authors, during the period of earthquake preparation, the generated seismoelectric fields may permeate the ionosphere and induce TEC variations within it, hence affecting the seismo-electromagnetic environments surrounding the epicentre. A few days prior to the 12 May 2008 $M_w$ = 7.9 Wenchuan earthquake, Zhao et al. [134], Liu et al. [39] and Pulinets and Ouzounov [37] reported that ionospheric GPS TEC enhancement and/or reduction in anomalies simultaneously appeared above the epicentre and its magnetic conjugate point.

Increased ionosphere observations from space and on Earth clearly show that there is a coupling mechanism between lithosphere-based seismic activity and ionosphere-based deviations or disturbances in electron concentrations, particularly prior to major earthquakes [121,122,124]. The measurements include variability in the critical frequency of the $F_2$ layer, $f_o{F}_2$ and TEC [124,135]. Compared to costly and sparse $f_o{F}_2$ observations using Earth- or space-based ionosondes, TEC measurements are more readily acquired with the use of global GPS TEC [124]. The impact of the earthquakes on the ionosphere is bigger for earthquakes with higher magnitudes and greater depths, hence affecting the relative density of TEC anomalies, within an radius area of 1000 km around the earthquake’s hypocentre [125]. Gulaeva and Arikan [125] suggest that the number of positive TEC storm anomalies is twice that of non-storm values, and that this observation supports dominant post-earthquake TEC enhancement with the ionosphere peak decreasing for 12 h during daytime but growing by night-time 6 h after the earthquake, followed by a gradual recovery afterwards.

According to Sorokin et al. [95], there are two possible causes for the TEC ionospheric anomalies: variations brought on by acoustic gravity waves and variations created by electric fields. Variations in the density of TEC are caused by a variety of natural events, including dust storms, thunderstorms, solar radiation, volcanic activity, radioactive gases and thunderstorms [136,137,138]. For instance, TEC increased during the 2014–2015 high-solar-radiation cycle, which was caused on by high-energy solar particles interacting with the Earth’s ionosphere resulting in TEC shifting [121]. Therefore, these parameters should be taken into account when studying TEC ionospheric variations.

3. Radon Precursors

3.1. Radon Properties

Radon (²²²$Rn$) is a natural radioactive noble gas. It is produced when radium (²²⁶$Ra$) decays. According to Nazaroff and Nero [26], there are thirty-nine known isotopes of radon, ranging from ¹⁹³$Rn$ to ²³¹$Rn$. Radon has a half-life of 3.823 days and is the most stable isotope. ²²⁰$Rn$, or thoron, has a half-life of 54.5 s. Due to its short half-life, thoron decays rapidly and because of this, it is often detected at low concentrations. However, that depends on the concentration of its parent nucleus (²²⁴$Ra$), especially in comparison to that of ²²⁶$Ra$. Radon is primarily responsible for the radioactivity present in the atmosphere at sea level [139].

Radon emissions mostly originate from soil [26]. About 10% of the radon that is diluted in soil gets released into the atmosphere [139]. In addition to soil, radon may be found in surface and underground waters, as well as fragmented rock [139,140]. While all radon atoms produced are diluted in fluids, only a portion of radon emerges from porous media and fractured rock, enters the volume of the pores and dissolves within the pore’s fluid [140]. Once there, either convection, advection or molecular diffusion can cause a macroscopic transport [26]. Interconnected pores and water aquifers allow this movement to appear [140]. Radon dissolves into the water present in the pores of soil and rock and is carried away by it [26]. The most crucial elements for these processes are the pressure differentials, the temperature gradients, and the permeability of soil [6].

Radon is a significant radiological risk factor since it contributes to over the half of the effective dose equivalent of the population’s exposure to natural sources and is the leading natural cause of lung cancer [141,142,143]. Due to this, radon is a subject of extensive research worldwide [141,142,143,144,145,146,147,148,149,150]. In addition to the above health risks, radon offers several beneficial uses in a variety of applications. In meteorology, the amount of uranium is calculated from the changes of radon’s emission in soil, and the obtained information is then utilised to monitor air masses. When assessing how accurate chemical transport models are in estimating greenhouse gas emissions, radon can be a helpful tracer for understanding how the atmosphere functions [151]. Radon has been employed in hydrological studies and for faults identification [152,153,154]. The shift in radon concentrations near faults and the anomalous variations in groundwater levels [155,156] have been employed as well, because the corresponding concentration variations reflect the water–rock interactions [157] and the pathways generated by active faults [158]. Radon has also been studied for co-seismic effects and tidal strain [159]. Radon’s half-life in association with its inert nature provides it with the ability to travel long distances without significant loss [27]. Because of this, radon has been extensively used to study tectonic activity [6,8,17,160,161]. Under this perspective, radon is the best among the various hydrogeological tracers for earthquake forecast.

Radon combines hydrological, geological and environmental properties. Hydrologically, it dilutes in water [26] molecules and water aquifers. It is present in surface and, most importantly, underground waters [139,140]. Geologically, it is easily transferred within soil and rock reaching areas away [2]. Environmentally, it is naturally emitted and present in atmosphere, soil and rock, surface and ground water. It is easily detected. All these properties combined have made radon one of the best precursors of seismic activity and the one with the longest history in earthquake-related studies [8,17,86,160,161].

3.2. Pre-Seismic Radon Anomalies

Abnormal radon changes before earthquakes have been found in groundwater, soil gas, atmosphere and thermal spas (e.g., [6,8,17,29,52,147,160,161,162,163,164,165,166,167,168,169,170,171,172,173]) and recently, have been associated with TEC variability (please see Section 2.5) [131,132,133,174]. There are considerable variations in the relationships between magnitude, precursory time and epicentral distance in connection with the range and number of radon anomalies and other features of the associated time series (e.g., [6,8,17,30,31]). For instance, the epicentral distances of earthquakes identified with the aid of radon vary from 10 km to 100 km, whereas the recorded precursory durations span from three months to a few days before the earthquake’s occurrence. Comparable ranges have also been published by Cicerone et al. [6], Ghosh et al. [17], Petraki et al. [161], Conti et al. [8] and Huang et al. [160]. Several precursory signals have been obtained with passive techniques (no external power needed), which offer rough time series estimations, since these methods integrate the radon concentrations over extended periods of time (of at least > 1 week), necessitated for the measurement. This coarseness poses significant bias to the precursory estimations. Nowadays, radon precursory signals are monitored with active techniques (external power required). The active techniques are implemented with portable monitors which allow for high sampling rates (typically between 1 ${\min }^{-1}$ and 1 ${\mathrm{h}}^{-1}$). As a consequence, these techniques offer detailed signals of radon and fine estimations (e.g., [2,6,8,17,160,161]). It is crucial to mention that additional factors influence the estimates of radon and earthquakes. For instance, seasonal fluctuations, geological and geophysical conditions, rainfall and changes in barometric pressure all have an impact on radon concentration levels [6,21,22,23,26,30,31,139,160]. Because of this, the associated time series data are typically shown alongside the precursory signals of radon. Most of the correlations between radon and earthquakes are based on small-to-intermediate-sized magnitudes. This further limits the calculations since as of right now, neither for mild earthquakes nor for powerful earthquakes does there exist a universal radon model that can be used to explain a particular impending seismic event ([1,35,82,83,175] and references therein).

3.2.1. Soil

The release of radon from soil is important for research on earthquake forecasting. Because of this, one of the key elements in forecasting strong earthquakes is the monitoring of radon emissions, and this is accomplished by various research groups [2,19,24,25,27,31,36,48,51,52,54,147,162,163,164,165,166,167,168,169,170,171,172,173,174,176,177,178,179,180,181]. The variability of radon’s emission linked to the seismic occurrences at a monitoring station determines how successful these investigations are. Radon concentration in soil depends on a number of parameters and thus varies between different natural environments. The objectives of the detection of radon pre-earthquake precursors are guided by the certain prospects of each region [49,50,161,182]. Traditionally, because of well-investigated relationships between radon and environmental parameters, any radon concentration deviations are believed to be indicative of changes brought by tectonic force during the earthquake preparation. In general, twice the standard deviation or more from the average soil radon concentration at a site of observation is thought to reflect appreciable anomalies. The radon anomalies are attributed to earthquake-related stress–strain changes underneath the Earth’s crust, but this has been a subject of significant argumentation [2,8]. Two approaches have been taken towards the study of soil radon as an earthquake precursor: one involves simulation experiments in the lab and the other involves monitoring the concentration of radon in soil gas over an extended period of time, in comparison to unusual emission changes with respect to seismic occurrences. In order to understand the gas dynamics underlying the ascent of radon from deep below the Earth’s crust to the surface, a number of in situ and laboratory experiments as well as models have been suggested [6,8,17,161].

Based on observations and conclusions drawn from all of the aforementioned studies, it has been determined that radon gas, which originates from the decay of radium in rocks inside the crust, partly stays in the crustal matrix, while the remainder migrates away through interconnected pores and aquifers using diffusion, fluid flow and alpha recoil. Changes in the strain field are caused by the accumulation of tectonic stress before an earthquake. According to Fleischer and Mogro-Campero [183], the deformation of rock mass under stress creates new channels that allow deep Earth gasses to ascend to the surface.

3.2.2. Groundwater

Although the idea that radon anomalies in groundwater may be connected to earthquakes was initially put up in 1927, the Great Tashkent earthquake of 1966 produced the first indication of an abnormally high radon concentration in groundwater [184]. Subsequently, a number of groups employed the concentration of radon in groundwater to study earthquakes [20,28,30,162,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200]. Groundwater radon concentrations frequently increase before earthquakes (e.g., [36,178,179,193,201]). However, in some cases, the amount of radon in groundwater decreased [20,186,187,193]. Significant earthquakes may be related to groundwater radon’s peculiar behaviour as it offers information about subsurface dynamics [179], particularly in areas where high-stress build-up occurs in the crust [202]. The processes driving seismic activity can be better understood by examining the links between seismicity and geochemical signal variability [179].

The route that groundwater follows underground, or the kinds of rocks and soil it encounters, determines the amount of radon that is present in the groundwater [203] or escapes from it [204]. The measurement of the underground water baseline is crucial because radon gas permeates the water from these rocks and soils and alters the amount of radon in these waters. Baseline radon concentrations in groundwater vary greatly. According to Kandari et al. [205], radon concentrations in 15 water samples from the Dehradun region, which is close to an active fault, ranged from 1.70 Bq ${\mathrm{L}}^{-1}$ to 7.57 Bq ${\mathrm{L}}^{-1}$. In southern Catalonia, 15 hot springs had groundwater radon levels ranging from 1.4 Bq ${\mathrm{L}}^{-1}$ to 105 Bq ${\mathrm{L}}^{-1}$ [206]. Using an AlphaGUARD system, Spanish researchers measured 28 groundwater samples collected from northeastern Gran Canaria (Canary Islands, Spain). They found that the highest and lowest levels of dissolved radon concentration were 76.9 Bq ${\mathrm{L}}^{-1}$ and 0.3 Bq ${\mathrm{L}}^{-1}$, respectively [207]. More range values are provided globally [141,142,143,208]. Significantly more radon is found in groundwater in thermal spas [139,209,210].

The seasonal fluctuation in groundwater radon concentrations may be attributed to temperature, precipitation, and other climatic conditions, but its anomalies may also be linked to shifts in tectonic stress [199,211]. While it is now well accepted that radon anomalies may be associated with earthquakes, anomalies are typically exceedingly hard to locate since variations in radon concentration frequently exhibit the features of nonlinear dynamic fluctuations. Thus, the development of efficient identification techniques is necessary. To some extent, the conventional statistical techniques are erroneous and subjective. A few data mining techniques, such as artificial neural networks and machine learning, have had some success recently [166,199,200,211].

3.2.3. Atmosphere

The primary source of atmospheric radon concentration is the exhalation from the Earth and to a lesser extent the escape from surface and subsurface water [26]. Numerous processes are involved, and meteorological elements have a significant impact on them [143]. Therefore, detecting anomalies in air radon in relation with earthquakes is significantly more challenging than detecting them in groundwater or soil radon. Prior studies computed anomalies in atmospheric radon concentrations by establishing a threshold level for the anomalies based on a normal variation period and removing the seasonal component anticipated from a sinusoidal model [164,212]. The results of these conventional methods depend on how the seasonal component is determined because the assessment is based on departures from the assumed sinusoidal model and the selected normal period of average fluctuations [213].

Japan is the primary source of studies on earthquake forecasting using atmospheric radon. Iwatata et al. [164] reported that anomalies in the atmospheric radon concentration were linked to the moment releases of large earthquakes based on ten years of continuous observation of the concentration over north-eastern Japan and Hokkaido. Yasuoka and Shinogi [214] reported that two months before the main shock of the 1995 Kobe earthquake ($M_w$ = 6.9; 17 January 1995, ${34.6}^{\circ }$ N, ${135.0}^{\circ }$ E), an increase in atmospheric radon concentration was noticed at Kobe Pharmaceutical University. Goto et al. [215] reported anomalous atmospheric radon concentrations associated with a shallow inland earthquake ($M_j$ = 5.5, depth = 7 km; 5 July 2011, ${34.0}^{\circ }$ N, ${135.2}^{\circ }$ E) in northern Wakayama. Yasuoka et al. [29] reported that the residual values for each day could be fitted very well to a log-periodic oscillation model by applying the exponential smoothing method to the fluctuations in the residual values. The authors stated that the residual values stopped increasing on 31 December 1994, and they concluded that this corresponded to the critical point of the best-fit model. These authors stated that rather than the main direct stresses causing the Kobe earthquake, the local stresses were responsible for the unusual ²²²$Rn$ fluctuation as well. Using the irreversible thermodynamic model, Kawada et al. [216] proposed that the preseismic radon shift was caused by a small change in crustal strain. Furthermore, a quantitative study by Omori et al. [213] revealed that the unusually high radon concentration (about 10 Bq ${\mathrm{m}}^{-3}$) before the Kobe earthquake increased air conductivity and was sufficient to produce ionospheric disturbances. Yasuoka et al. [217] claimed that further mechanically induced precursors were seen prior to the Kobe earthquake. Igarashi et al. [189], for instance, described such precursory variations in groundwater radon concentration. Tsunogai and Wakita [218] documented further preseismic variations in crustal strain, groundwater discharge rate and chloride ion content in groundwater. Because of the mechanical behaviour of the crust, these preseismic fluctuations should be related to one another [29,189,216]. The fact that the temporal change in atmospheric radon concentration has not been compared with that in other preseismic events was noted by Igarashi et al. [189]. The linkage between preseismic fluctuations in the subsurface, atmosphere and ionosphere could have been substantially verified if radon activity had been clearly linked to the earthquake preparation process [189]. Additionally, current research supports the link between atmospheric radon and the Kobe earthquake [215,219].

4. Models

4.1. Electromagnetic Precursors Models

4.1.1. Models for the ULF Precursors

There are three main models that have been proposed for the interpretation of the magnetic component of the preseismic ULF disturbances:

• Magneto-hydrodynamic model [220]: According to this model, an electrically conducting fluid flowing through a magnetic field causes an additional induced field to be created. If B is the magnetic field, Maxwell’s equations indicate that the induced magnetic field $B_i$ can be given by the equation $B_i=R_m·B$, where $R_m$ is the magnetic Reynolds number, comparable to the hydrodynamic Reynolds number, which determines the relative significance of the convective and diffusive components.
• Piezomagnetic model [221]: this model suggests that an applied stress causes ferromagnetic rocks to shift in magnetisation, which in turn, induce a secondary magnetic field.
• Electrokinetic model [222]: this model suggests that electric currents flowing in the Earth due to electrified interfaces present at solid–liquid boundaries induce magnetic fields.

Varotsos et al. [223] proposed a theory about the current produced by charged distortions and currents induced by piezo-electric effects. The electrokinetic theory served as the foundation for this theory. In water-saturated media with fluid-filled channels, electrokinetic currents can be found [224,225]. In order to model the parameters of these electrokinetic currents, Surkov et al. [226] assumed that an earthquake hypocentre was surrounded by water-saturated porous rocks with fluid-filled pore channels, where cations from the fluid were adsorbed by the walls of pores and cracks in the solid material. According to this author, the fluid moving along the channel carries anions and as a consequence produces an extrinsic electric current between the fluid and the surrounding walls.

When an earthquake is imminent, the seismic hypocentre within the Earth’s crust is surrounded by cracks and fractured material, where new fractures are continuously produced forming the so-called fracture zone. The fracture zone can range in size from a few hundred metres to several kilometres. Feder [227] postulated that there was a fractal structure present in the pore’s space within the fracture zone. Newly developed cracks are sealed off as soon as they arise under reduced pressure, as a result of the pressure release that is caused by cracking. This, in turn, allows water from the uncracked outside zone to enter as soon as a network of linked channels or fractal clusters is formed. This can be seen, alternatively, as a grid of new cracks that are closed as the water sinks from the nearby locations of greater pressure. According to Surkov et al. [226], during cluster formation, the porosity and permeability of rocks decrease from the centre of the fracture zone towards the perimeter. An interior area manages to surpass the percolation threshold and due to this, the permeability outside the fracture zone tends to zero. In actuality, there is a limited permeability since crustal rocks have a large variety of interconnectible small cracks. Furthermore, the rock’s conductivity together with the surface and bulk conductivities of the tiny fluid-filled cracks contribute to the non-zero conductivity of the surrounding space. However, according to Surkov et al. [226], the conductivity beyond the fracture zone is minimal. This indicates that because of the recently formed fluid-filled cracks, the conductivity’s value is more closely tied to the conductivity of the percolation threshold. It is important to note that only the percolation hypothesis can adequately explain the range of fracture diameters. Surkov et al. [226] limited the study by using a basic percolation hypothesis that ignored the crack-channel size distribution. The correlation length $\xi$ is $\xi =1{\left|p-p_c\right|}^{\nu }$ with p being the probability that a channel can conduct the fluid, $p_c$ is the critical probability in the percolation threshold and $\nu$ = 0.88.

The aforementioned three ULF models have successfully described major earthquakes identified with ULF data: the $M_w$ = 9.0 earthquake at Tohoku, Japan; the $M_w$ = 8.3 earthquake at Coquimbo, Chile; the $M_w$ = 8.1 earthquake at Chiapas, Mexico; and the Vrancea seismicity in Romania [29,49,165,186,189,190,214,217,219,228].

4.1.2. Models for the HF Precursors

The behaviour of a stressed rock is comparable to that of an electromagnetically strained rock [175]. The crack propagation is the basic process responsible for the material’s failure [83]. The release of photons, electrons, ions, and neutral particles is observed when fracture, deformation, wearing and peeling cause new surface characteristics to appear in various materials [1,82,83,90,229]. The total of these emissions are collectively referred to as fracto-emissions [83]. The significant charge separation brought on by the rupture of the inter-atomic ionic bonds is the source of the electric charge between the micro-crack faces. An electric dipole or a more intricate system is created by the electric charges on the surfaces of freshly developed micro-cracks. It has been shown that a dynamical instability controlling the oscillations in the velocity and shape of a crack on the fracture surface controls the crack’s mobility [83].

According to experimental data, micro-fracturing events repeat and intensify until a multi-crack state occurs, indicating that local branching is the instability mechanism at work. It is important to note that laboratory research has identified strong fracto-emissions during unstable crack propagation [22,23,34,83,230]. Because of the intense wall vibrations of the cracks during the micro-branching instability stage, the cracked material functions as an efficient emitter. As a result, the opening of cracks in a material can be seen as a potential precursor of general fracture because electromagnetic emissions occur in a wide frequency range from kHz to MHz when the material is stretched. These electromagnetic precursor are detected during in-field measurement and in laboratories under controlled conditions [1,34,82,90,91,175]. Consequently, the main technique for forecasting earthquakes is to record the electromagnetic emissions from potential microfractures in the focal region prior to the final break-up [1].

As stated in several papers (e.g., [82,89,175,229] and references therein), a “symmetry breaking” is linked to a thermal second-order phase transition. For non-equilibrium irreversible processes, the evolution of the “symmetry breaking” with time was reported in order to obtain an understanding of the catastrophic nature of the fracture events. The investigation revealed that the system’s balance was progressively lost. This allowed for the estimation of the duration beyond which the process responsible for the preseismic electromagnetic emissions could continue as a non-equilibrium instability.

The analysis indicated three key periods: (i) the crucial epoch, also known as the critical window, in which the short-range correlations transit to long-range ones; (ii) the “symmetry breaking” epoch; and (iii) the integration of the “symmetry breaking.” It is widely acknowledged that a notable rise in localisation and directionality occurs at the terminal phase of the earthquake preparedness procedure. Therefore, it is critical to identify distinctive epochs in the precursory electromagnetic activity progression and to connect these to the corresponding final phases of the earthquake preparation process.

Tracing “symmetry-breaking” could indicate that the focal area’s heterogeneous component, which encircles the fault plane’s strong asperities’ backbone, has reached the point of microfracture propagation completion. At this point, the rupture becomes blocked at the boundary of the strong asperities’ backbone. Asperities are already under “siege” [83].

4.2. Radon Precursors Models

Scholz et al. [231] presented the Dilatancy–Diffusion model, which connects anomalous radon changes to the mechanical crack development rate in the volume of a dilatancy, so as to simulate the underlying dynamics of radon prior to earthquakes. This model states that the first medium is a porous, fractured, submerged rock. Favourably placed fractures open when tectonic forces grow because the cracks expand and disengage close to the pores. As a result, the preparation zones’ overall pore pressure decreases, allowing water from the surrounding medium to enter the zone. Radon emission may fluctuate suddenly as a result of the pore pressure returning and the number of cracks growing. The crack-avalanche model [176,232] states that the increase in tectonic stress forms a fractured focal rock zone. This zone gradually changes in volume and form over time. According to the hypothesis of stress corrosion [233], the slow crack propagation, which is controlled by stress corrosion in the rock matrix saturated by groundwater, may be linked to the unusual behaviour of radon concentration. The LAIC model [37,94] also describes radon’s stress accumulation in the ground. This is attributed to the relative movement of tectonic blocks, which, in turn, lead to the formation of micro-cracks, cracks, and fractures. Radon gas released from microfractures combines with water and travels via various media to the earth. Water and carrier gases are often responsible for the transportation of radon from the Earth’s deep strata to the surface [234].

Nikolopoulos et al. [2,21,32,35,235], Petraki et al. [22,23], Alam et al. [36,178,179,201] and Petraki [86] proposed the asperity model [83] (please see Section 4.1.2) to explain radon emanation during the preparation of earthquakes. Preseismic radon anomalies are attributed to variations in fractional Brownian (fBm) profile movements. In the views described in Section 4.1.2, the focal area consists of a backbone of strong and large asperities that sustain the system and a strongly heterogeneous medium which surrounds it. The fracture of the heterogeneous system in the focal area obstructs the backbone of asperities. As the fracture becomes persistent, strong antipersistent radon concentration anomalies occur interchangeably. This has been associated with several earthquakes in Greece and China.

Other aspects have been expressed by other investigators. For example, Talwani et al. [236] reported that the anomalous behaviour of radon gas could be because of the opening of pores’ spaces during rock fracturing as a result of seismic events. Explosion tests were performed to identify the relationship between the dynamic loading effect and the observed concentrations of radon [18]. The experimental results revealed that the increase in radon values was a consequence of seismic waves applied to the rock. According to other investigators [237,238,239], crustal activities have been identified as one of the reasons for radon emission.

5. Analysis Methods

Several investigations on earthquake forecast have been based on visual observations [6,11,161]. Despite providing some indications, the visual observations are not enough to support the preseismic nature of the derived signals (e.g., [2,86] and references therein). Due to this, the analysis nowadays rely on the physical background of the related earthquake processes. The mainstream analysis comprises fractal methods [1,34,35,36,55,60,84,85,90,91,101,147,158,178,230,235,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255], methods based on the theory of information and entropy [82,89,175,229], symbolic dynamics [21,23,86,256,257,258,259] and natural time methods [64,260,261,262,263,264]. Within the above framework, several metrics have been utilised as adequate for the related analysis. These metrics comprise exponents from the spectral power law (e.g., [34,84,254]), Detrended Fluctuation Analysis (DFA) (e.g., [1,235]), Rescaled-Range Analysis (R/S) (e.g., [265]), Multifractal Detrended Fluctuation Analysis (MFDFA) (e.g., [36,55]), fractal dimensions from Katz’s, Sevcik’s and Higuchi’s methods [34,147], Hurst exponents and entropy values from (i) entropy per letter; (ii) conditional entropy; (iii) entropy of the source; (iv) t-entropy; (v) Tsallis entropy; (vi) perturbation entropy; (vii) normalised Tsallis entropy and parameters for critical phenomena (e.g., [1,82]).

Spectral power-law analysis and Hurst exponent analysis have been utilised in all ULF, HF and radon precursors. DFA, fractal dimensions from Katz’s, Sevcik’s and Higuchi’s methods and the $R/S$ analysis have been used with success both for HF and radon precursors. Symbolic dynamics with entropy per letter, Tsallis entropy and normalised Tsallis entropy have also been employed for both HF and radon precursors but to a lesser degree. Natural time has been employed mainly in ULF signals. The remaining techniques and metrics have been used mainly for HF precursors. Multifractal Detrended Fluctuation Analysis (MFDFA) has been employed in all types of precursors [55,179,241,255,266,267,268] but is not presented here due to its complicated interpretation [269].

Due to their importance in both electromagnetic and radon precursors, the important properties of fractal behaviour, long memory and Hurst exponent analysis are first given in the following sub-sections. Then, DFA is presented because it is a robust method that has been used in both LF, HF and radon precursors. The fractal dimension calculations through Katz’s, Sevcik’s and Higuchi’s methods are given next because they have been utilised both in HF electromagnetic and radon precursors and finally, the $R/S$ analysis, because it is the main direct method to calculate Hurst exponents and has been employed both in HF electromagnetic and radon precursors.

5.1. Important Properties: Fractal Behaviour, Long Memory and Hurst Exponents

5.1.1. Fractal Behaviour

Many physical systems in nature display a fractal behaviour, which is reflected when these systems are stretched, translated or rotated in space. Based on their mathematical characteristics, these systems are classified as either self-similar or self-affine. These systems are fractals because each component of the system is a large-scale imitation or representation of the system as a whole due to the self-affinity and self-similarity that define all system components. This characteristic allows for the investigation of fractal systems through a part-by-part analysis. System fractals can exhibit self-similarity or self-affinity. While self-affine systems behave almost in this way, self-similar systems have exact inter-parts representations.

The system’s complexity [270], which indicates whether the system is driven by linear mechanisms and order [271,272], is also connected with the scaling and fractal behaviour. The correlations are strong because a system’s complex behaviour may be predicted by its fractal behaviour and vice versa.

5.1.2. Long Memory

The long memory [270,273,274] of a system can show if the system has long-range interactions or is random. Specifically, it may reveal if a geo-system has strong persistent and antipersistent behaviour or if the long-range interactions are rather loose. If a system exhibits long memory, then the past, present and future states of the system are linked together in a manner that the present state of the system is not only derived from its past (Markovian behaviour) but also defines its future (non-Markovian behaviour) [82,175]. This behaviour is characteristically seen when the fracture of the Earth’s crust yield to the inevitable general breakdown during the unstoppable approaching of an ensuing earthquake [1,83,90,91]. Precisely there, the past determines the presence and also the inevitable future breakdown of the system.

5.1.3. Hurst Exponent

Because it may depict enduring connections in space or time, the Hurst exponent (H) provides a straightforward technique for assessing a system’s long memory [275,276]. Time-evolving fractal events may be identified with the Hurst exponent, and the corresponding time series’ roughness can be evaluated [277]. Important details about the time series are revealed by the Hurst exponent’s value [241,275,276,278,279]:

(i)

The series has positive long-range autocorrelation if $0.5<H\le 1$. A series’ high value is followed by another high value and vice versa. High Hurst exponents suggest persistent interactions that are anticipated to remain until the series’ remote future;

(ii)

Low values of the time series follow high values if $0\le H<0.5$, and vice versa. In the future of the time series, there is a persistent transition between low and high values for low H values (antipersistency);

(iii)

If $H=0.5$, the time series is completely uncorrelated, i.e., the related processes are random.

5.2. Significant Analysis Methods for Electromagnetic and Radon Precursors

5.2.1. Power-Law Analysis

In the event that a temporal fractal is present in the time series, the power spectral density, $S\left(f\right)$, exhibits a power-law behaviour:

$$S\left(f\right)=a·f^{-\beta }$$

(1)

In Equation (1), a represents the spectral density amplification, f denotes a transform’s frequency, and $\beta$ is the power-law exponent, which measures the strength of the power-law associations. This transform can be the wavelet transform [84] or the FFT of the signal [254,255]. Given its perceived benefits, the wavelet transform based on the Morlet base function is the most frequently employed [2,34,82,84,175,230,235,245]. In particular, f represents the central frequency of the Morlet wavelet.

Equation (1)’s logarithmic transformation yields:

$$logS\left(f\right)=loga+\beta ·logf$$

(2)

Given that Equation (1) is a straight line, $\beta$ and a may be found by using the least-squares approach to fit the associated data.

The technique has been utilised mostly in sliding windows of various lengths moved one sample forward. Independent windows are also utilised as well, under the restriction that the square of Spearman’s ($r^2$) coefficient in each window should have $r^2\ge$ 0.95 for the power-law fit to be acceptable.

5.2.2. DFA

The original time signal is first integrated in order to apply DFA. Then, within a window of size n, the integrated signal’s fluctuations, $F\left(n\right)$, are found. The linear $log\left(F\left(n\right)\right)-log\left(n\right)$ transformation is then fitted using least squares to obtain the integrated time series’ scaling exponent (self-similarity parameter), $\alpha$. Depending on the dynamics of the system, the $log\left(F\left(n\right)\right)-log\left(n\right)$ line may show one crossover at a scale n where the slope displays an abrupt shift, two crossovers at two distinct scales $n_1$ and $n_2$ [86] or nothing at all.

The following process may be used to construct the DFA of a one-dimensional temporal signal $y_i$, ($i=1,\dots ,N$) [34,86,280]:

(i)

First, the original time series is integrated:

$$y\left(k\right)=\sum _{i=1}^k\left(y\left(i\right)-\langle y\rangle \right)$$

(3)

In Equation (3), the symbols <…> represent the total average value of the time series, whereas k represents the different time scales.

(ii)

Next, the integrated time series $y\left(k\right)$ is divided into equal-length bins, n, which do not overlap.

(iii)

The trend in the bin is subsequently expressed by the function $y\left(k\right)$, which is then fitted. Simple linear trends or polynomials of order two or a higher order may be used. The notation $y_n\left(k\right)$ indicates the y coordinate of this linear function in each box n.

(iv)

Next, each box of length n is detrended in the integrated time series $y\left(k\right)$ by subtracting the local linear trend, $y_n\left(k\right)$. In this way, and for every bin, the detrended time series $y_d^n\left(k\right)$ is calculated as follows:

$$y_d^n\left(k\right)=y\left(k\right)-y_n\left(k\right)$$

(4)

(v)

Next, for each bin of size n, the root-mean-square (rms) of the integrated and detrended time series fluctuations is calculated as

$$F\left(n\right)=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{k=1}^N}{\left\{y\left(k\right)-y_d^n\left(k\right)\right\}}^2}$$

(5)

where $F\left(n\right)$ are the rms fluctuations in the detrended time series $y_d^n\left(k\right)$.

(vi)

The technique steps (i)–(v) are repeated for different sizes $\left(n\right)$ of the scale boxes. This indicates the precise of a kind of relationship that exists between $F\left(n\right)$ and n. An exponential relationship exists between $F\left(n\right)$ and n if the time series contains long-term associations.

$$F\left(n\right)\sim n^{\alpha }$$

(6)

The DFA scaling exponent $\alpha$ of Equation (6) assesses the strength of the time series long-term relationships.

(vii)

Equation (4)’s logarithmic translation yields a linear relationship between $logF\left(n\right)$ and $log\left(n\right)$. A strong linear relationship implies that the accompanying fluctuations have a long memory since they are long-lasting. This study used the square of the Spearman’s ($r^2$) to assess the linear fit’s accuracy. According to Nikolopoulos et al. [34,230,235,280], good linear fits are considered as having $r^2\ge$ 0.95 or higher.

As with Section 5.2.1, DFA has also been utilised in sliding windows of various lengths moved one sample forward.

5.2.3. Fractal Dimension Analysis with Katz’s Method

The transpose array ${[s_1,s_2,\dots ,s_N]}^{⊺}$ of the series $s_i$, $i=1,2,\dots ,N$, is first determined in accordance with Katz’s method, where $s_i=(t_i,y_i)$ and $y_i$ are the measured series values at the time instances $t_i$ [281,282]. This process yields the fractal dimension D.

The two subsequent points of the time series ($s_i$ and $s_{i+1}$) are represented by the value pairs ($t_i,y_i$) and ($t_{i+1},y_{i+1}$), for which the Euclidean distance is:

$$dist(s_i,s_{i+1})=\sqrt{\left(t_i^2-t_{i+1}^2\right)+\left(y_i^2-y_{i+1}^2\right)}$$

(7)

The distances in Equation (7) add up in a curve whose total length is:

$$L=\sum _{i=1}^{i=N}dist(s_i,s_{i+1})$$

(8)

This curve stretches in the planar to d, if it does not cross itself, where d is as follows:

$$d=max\left(dist(s_i,s_{i+1})\right),i=2,3,\dots ,N$$

(9)

By combining Equations (7)–(9), the Katz fractal dimension, D, becomes

$$D={\displaystyle \frac{log\left(n\right)}{log\left(n\right)+log(d/L)}}$$

(10)

where $n=L/\overline{\mathrm{a}}$, and $\overline{\mathrm{a}}$ is the average value of the distances of the points.

5.2.4. Fractal Dimension Analysis with Higuchi’s Method

To determine a time series’ fractal dimension, D,

$$y\left(1\right),y\left(2\right),y\left(3\right),\dots ,y\left(N\right)$$

(11)

recorded at $i=1,2\dots N$ intervals, the following is the construction of a new sequence, $y_m^k$ [246,247,283]:

$$y_m^k:y\left(m\right),y(m+k),y(m+2k),\dots ,y(m+\left[{\displaystyle \frac{N-m}{k}}\right]k)$$

(12)

The length of the curve associated with the time series is given by [283]:

$$L_m\left(k\right)={\displaystyle \frac{1}{k}}\left(\sum _{i=1}^{\left[{\displaystyle \frac{N-m}{k}}\right]}y(m+ik)-y(m+(i-1)k)\right)\left({\displaystyle \frac{N-1}{{\left[\frac{N-m}{k}\right]}^k}}\right)$$

(13)

In both equations, m and k are integers that specify the time interval between the series samples and are connected by the formula $m=1,2\dots k$, where $\left[\dots \right]$ is the Gauss notation, namely, the bigger integer part of the included value.

The following normalisation factor is inserted.

$${\displaystyle \frac{N-1}{{\left[\frac{N-m}{k}\right]}^k}}$$

(14)

The lengths of Equation (14) show an average value, $〈L\left(k\right)〉$, that displays a power law of the following form:

$$〈L\left(k\right)〉\propto k^{-D}$$

(15)

Higuchi’ s fractal dimension, D, is finally calculated by the slope of the linear regression of the logarithmic transformation of $〈L\left(k\right)〉$ versus k, where $k=1,2,\dots ,k_{max}$. It must be noted that the time intervals are $k=1,\dots ,k_{max}$ for $k_{max}\le 4$, i.e., $k=1,2,3,4$, for $k_{max}=4$ and $k=\left[2^{(j-1)/4}\right]$, $j=11,12,13\dots ,$ for $k>4$ ($k_{max}>4$). Again, $\left[\dots \right]$ is the Gauss notation [282].

5.2.5. Fractal Dimension Analysis with Sevcik’s Method

Using Sevcik’s approach [284], the fractal dimension of a time series is estimated from the Hausdorff dimension, $D_h$, as [282].

$${\displaystyle D_h=\underset{ϵ\to 0}{lim}\left[-\frac{log\left(N\right(ϵ\left)\right)}{log\left(ϵ\right)}\right]}$$

(16)

where $N\left(ϵ\right)$ is the total number of $ϵ$-length segments that together form a curve related to the time series. $N\left(ϵ\right)=L/2ϵ$ [282] and $D_h$ are as follows if the length of the curve is L:

$${\displaystyle D_h=\underset{ϵ\to 0}{lim}\left[-\frac{log\left(L\right)-log\left(2ϵ\right)}{log\left(ϵ\right)}\right]}$$

(17)

The N points of the curve L can be mapped to a unit square of $N\times N$ cells of the normalized metric space by twice performing a linear transformation. Equation (18) yields Sevcik’s fractal dimension with this transformation [282,284]:

$${\displaystyle D_h=\underset{N\to \infty }{lim}\left[1+\frac{log\left(L\right)-log\left(2ϵ\right)}{log\left(2\right(N-1\left)\right)}\right]}$$

(18)

The calculation improves as $N\to \infty$.

5.2.6. Rescaled Range Analysis

In order to identify trends that could recur in the future, the $R/S$ analysis uses two variables: the range, R, and the standard deviation, S, of the data [275,276]. In accordance with the $R/S$ technique, the average, ${\langle x\rangle }_N={\displaystyle \frac{1}{N}}{\sum }_{n=1}^{N}x\left(n\right)$, over a period of N time units transforms a time series $X\left(N\right)=x\left(1\right),x\left(2\right),\dots ,x\left(N\right)$ into a new variable $y(n,N)$ in a specific time period $n,(n=1,2,\dots ,N)$. The so-called cumulative deviation of the time series, $y(n,N)$, has the following formula:

$$y(n,N)=\sum _{i=1}^n(x\left(i\right)-{\langle x\rangle }_N)$$

(19)

The rescaled range is calculated as [86,275,276]:

$$R/S={\displaystyle \frac{R\left(n\right)}{S\left(n\right)}}$$

(20)

The distance between the lowest and largest value of $y(n,N)$ defines the range $R\left(n\right)$ in:

$$R\left(n\right)=\underset{1⩽n⩽N}{\mathit{max}}y(n,N)-\underset{1⩽n⩽N}{\mathit{min}}y(n,N)$$

(21)

The standard deviation $S\left(n\right)$ is calculated as follows:

$$S\left(n\right)=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{n=1}^N}{(x\left(n\right)-{\langle x\rangle }_N)}^2}$$

(22)

$R/S$ exhibits a power-law dependence on the bin size n

$${\displaystyle \frac{R\left(n\right)}{S\left(n\right)}}=C·n^H$$

(23)

where H is the Hurst exponent, and C is a proportionality constant.

The final equation’s log transformation is a linear relationship:

$$log({\displaystyle \frac{R\left(n\right)}{S\left(n\right)}})=log\left(C\right)+H·log\left(n\right)$$

(24)

This is used to directly calculate the Hurst exponent H, which is the slope of the best line fit. It is important to note that the only direct method to calculate Hurst exponents is via the $R/S$ analysis.

6. Precursors and Earthquake-Related Parameters

Several attempts have been made to link earthquake-related parameters and data derived from precursors. There is a variety of empirical relationships between earthquake magnitudes, preparation zone areas, precursory time and other earthquake-related characteristics. Some of these empirical relations are given in the following.

Rikitake [285] proposed a model showing the relations between anomaly, the precursory time T in days, the magnitude of an earthquake m and distance from epicentre R in km. According to this model:

$$logT=0.76·M-1.83$$

(25)

Talwani [286] suggested an empirical earthquake forecast model as follows:

$$M_L=logD-0.07$$

(26)

where $M_L$ is the local magnitude of an earthquake, and D is the forecasting period in days.

Guha [287] provided another model associating the precursory time T in days and the magnitude m of an earthquake:

$$logT=A+B·M$$

(27)

where A and B are statistically determined coefficients.

Dobrovolsky et al. [288] proposed an empirical relationship for the calculation of the earthquake preparation zone $R_D$ (km) and the magnitude (m) of the ensuing earthquake:

$$R_D={10}^{0.43·M}$$

(28)

Fleischer [289] suggested that the epicentral distance D in km and the magnitude m of an earthquake were associated as follows:

$$D=(1/1.66)·{10}^{A·M}$$

(29)

where A = 0.813 for m < 3 and A = 0.480 for m > 3.

Fleischer and Morgo Campero [290] suggested that

$$X_M={10}^{0.48·M}$$

(30)

where $x_M$ is the dislocation range in km, and m is the magnitude of an earthquake, where M ≥ 3.

Virk [291] proposed a different relation that combined the epicentral distance D in km and the magnitude m of an earthquake:

$$D={10}^{A·M}$$

(31)

where A = 0.32 for 10 km < D < 50 km, A = 0.43 for 50 km < D < 100 km, A = 0.56 for 100 km < D < 500 km and A = 0.63 for 500 km < D < 1250 km.

The epicentral distance, $R_E$, in km between a monitoring site and the earthquake’s epicentre can be calculated by

$$R_E=D·R$$

(32)

where R is the Earth’s radius (6370 km) and

$$D=cos{\alpha }_i·cos{\alpha }_j+sin{\alpha }_i·sin{\alpha }_j·(cos({\beta }_i-{\beta }_j)$$

(33)

where $({\alpha }_i,{\beta }_i)$ are the coordinates of the earthquake, and $({\alpha }_j,{\beta }_j)$ are the coordinates of the monitoring station [179].

Chetia et al. [163] used multiple linear regressions to examine the greatest variability caused by pressure, temperature and rainfall in soil gas radon. They suggested that the precursory time T (days), epicentral distance D (km) and magnitude m ($M_w$) were connected with the relationship

$$log\left(DT\right)=0.79·M+b$$

(34)

where b equals 0.18, a is approximately 3.51, and D equals roughly to $D100.58M$.

The reader may recall in relation to the estimations given in this section that there is no one-to-one correspondence between recorded anomalies and occurrence of an earthquake [1]. Moreover, the earthquake generation processes are multi-facet [1] and therefore, a combination of techniques is needed [1,34,36,82,269] to increase the scientific evidence. In view of these references, the estimations presented in this section have significant limitations. On the other hand, several papers of the previous decades, but also modern ones, make use of these estimations. For several scientist, these are considered as adequate and sufficient.

7. Table of Papers

Table 1. Papers on electromagnetic precursors. The papers are presented in chronological order from the oldest to the newest. The precursory time also includes the aftershock data presented in some papers. ED is the effective-sensitive distance between the monitoring site and the epicentre of the earthquake. Blank cells indicate there is no information available in the reference(s).

Location	Magnitude	Date(s)	Emission Type	Measurement Frequency	Instrumentation	Method(s)	Precursory Time	ED	Reference
Chile	9.5	22/05/1960	Radio	18 MHz	Radioastronomy receiver	Visual observation	6 days	Worldwide	[292]
Hollister, California	5.2	28/11/1974	ULF magnetic		Array of 7 proton magnetometers	Visual observation	7 weeks–several months	11 km	[293]
Tangshan, China	7.8	28/07/1976	Resistivity			Visual observation	2–3 years	<150 km	[294]
Tangshan, China	7.8	28/07/1976	Self-potential and magnetotelluric			Visual observation	3 months	<120 km	[294]
Sungpan-Pingwu, China	7.2	16/08/1976	Telluric currents			Visual observation	1 month	<200 km	[295]
Sungpan-Pingwu, China	6.8	22/08/1976	Telluric currents			Visual observation	1 month	<200 km	[295]
Sungpan-Pingwu, China	7.2	23/08/1976	Telluric currents			Visual observation	1 month	<200 km	[295]
Kyoto, Japan	7.0	31/03/1980	LF electric	81 kHz	Electric antenna	Visual observation	0.5 h	250 km	[43]
Tokyo, Japan	5.3	25/09/1980	LF electric	81 kHz	Electric antenna	Visual observation	1 h	55 km	[43]
Tokyo, Japan	5.5	28/01/1981	LF electric	81 kHz	Electric antenna	Visual observation	3/4 h	50 km	[43]
Kalamata, Greece	6.2	13/09/1986	Electric			Visual observation	3–5 days	200 km	[296]
Spitak, Armenia	6.9 ($M_s$)	07/12/1988	ULF magnetic	0.01–1 Hz	3-Axis magnetometers	Visual observation, statistical analysis	4 h	128 km	[7]
Spitak, Armenia	6.9 ($M_s$)	07/12/1988	ULF magnetic	0.005–1 Hz	3-Axis magnetometers	Visual observation, statistical analysis	4 h	120 km, 200 km	[65]
Loma Prieta, California	7.1 ($M_s$)	18/11/1989	ULF magnetic	0.01 Hz		Visual observation, statistical analysis	3 h	7 km	[65]
Loma Prieta, California	7.1 ($M_s$)	19/11/1989	ULF, HF electromagnetic	0.01 Hz, 32 kHz	Ground-based magnetometers	Visual observation	3 h	52 km	[58]
Spitak, Armenia	6.9 ($M_s$)	23/01/89	LF to HF electromagnetic	140, 450, 800, 4500, 15,000 Hz	COSMOS-1809 satellite with 12 satellite orbits of f < 450 Hz	Visual observation, FFT	<3 h		[297]
Upland, California	4.3	17/04/1990	ULF magnetic	3–4 Hz	Vertical magnetic sensor	Power law, FFT	1 day	160 km	[298]
West Iran	7.5	20/06/1990	Ionospheric radiowave	0–8 kHz, 10–14 kHz	INTERCOSMOS-19 satellite	Visual observation, modelling	16 days	250–2000 km	[299]
Watsonville, California	4.3	23/03/1991	ULF magnetic	3.0–4.0 Hz	North–south magnetic sensor	Statistical analysis, power law with FFT	Data averaged over 2 days	600 km	[298]
Watsonville, California	4.3	23/03/1991	ULF magnetic	3.0–4.0 Hz	Vertical magnetic sensor	Power law-FFT	Data averaged over 2 days	600 km	[298]
NW Crete, Greece	6.0	21/11/1992	HF electric	41, 53 MHz	Electric dipole antennas	Visual observation	1–3 days	20–150 km	[300]
Coalinga, California	4.0	15/01/1992	ULF magnetic	3.0–4.0 Hz	Vertical magnetic sensor	Power law–FFT	Data averaged over 2 days	400 km	[298]
Hokkaido, Japan	7.8	12/07/1993	foF2 ionospheric			Visual observation, statistical analysis	3 days	290 km, 780 km, 1280 km	[135]
Guam	7.1 ($M_s$)	08/08/1993	ULF magnetic	0.02–0.05 Hz	3-Axis ring core type fluxgate magnetometers	Fractal analysis, FFT	1 month	65 km	[60,301]
Guam	8.3 ($M_J$)	08/08/1993	ULF magnetic	0.02–0.05 Hz	3-Axis ring core type, fluxgate magnetometers	Multifractal Detrended Fluctuation Analysis	1 month	65 km	[241]
Hokkaido, Japan	8.2 (MJMA)	07/12/1993	SES	≤1 Hz	Electric antennas	Natural time analysis	1 month	$\delta$lat and $\delta$long < 30	[79]
Hokkaido-Toho Oki, Japan	8.1 ($M_W$)	04/10/1994	HF electric		Borehole antenna	Visual observation	20 min	>1000 km	[278]
Hokkaido, Japan	7.6 (MJMA)	04/10/1994	SES	≤1 Hz	Electric antennas	Natural time analysis	1 month	$\delta$lat and $\delta$long < 30	[79]
Hokkaido, Japan	7.4 (MJMA)	28/12/1994	SES	≤1 Hz	Electric antennas	Natural time analysis	1 month	$\delta$lat, $\delta$long < 30	[79]
Hyogo-ken Nanbu (Kobe), Japan	7.2 (MJMA)	17/01/1995	HF electric	22.2 MHz	Phase-switched interferometer polarized antennas		1 h	77 km	[302]
NE Samos, Greece	5.0	07/05/1995	HF electric	41, 53 MHz	Electric dipole antennas	Visual observation	1–3 days	20–150 km	[300]
Kozani-Grevena, Greece	6.6 ($M_W$)	13/05/1995	HF electric, LF magnetic				2 weeks	70 km, 200 km	[303]
Kozani-Grevena, Greece	6.6 ($M_W$)	13/05/1995	HF electric	41, 54 MHz, magnetic 3, 10 kHz	Electric dipole and magnetic loop antennas	Fractal analysis	20 h	284 km	[304,305]
Kozani-Grevena, Greece	6.6 ($M_W$)	13/05/1995	HF electric	41, 54 MHz, magnetic 3, 10 kHz	Electric dipole and magnetic loop antennas	Fractal analysis and statistical methods.	20 h	284 km	[305]
Kozani-Grevena, Greece	6.6 ($M_W$)	13/05/1995	HF electric	41, 54 MHz, magnetic 3 kHz	Electric dipole and magnetic loop antennas	Fractal analysis and statistical methods.	20 h	284 km	[84]
Kozani-Grevena, Greece	6.6 ($M_W$)	13/05/1995	HF electric	41 MHz	Electric dipole and magnetic loop antennas	Fractal analysis and statistical methods.	20 h	284 km	[306]
Kozani-Grevena, Greece	6.6 ($M_W$)	13/05/1995	HF electric and LF magnetic	41, 54 MHz and 3, 10 kHz	Electric dipole and magnetic loop antennas	Intermittent dynamics of critical fluctuations	20 h	284 km	[307]
Kozani-Grevena, Greece	6.6 ($M_W$)	13/05/1995	SES	≤1 Hz	Electric antennas	Visual and mathematical analysis	4 weeks	70–80 km	[308,309]
Kozani-Grevena, Greece	6.8 ($M_S$)	13/05/1995	SES	≤1 Hz	Electric antennas	Visual and mathematical analysis	24, 25 days	70–80 km	[309]
Kozani-Grevena, Greece	6.8 ($M_S$)	13/05/1995	SES	≤1 Hz	Electric antennas	Visual and mathematical analysis	22 min	70–80 km	[310]
SE Crete, Greece	5.0	29/07/1995	HF electric	41, 53 MHz	Electric dipole antennas	Visual observation	1–3 days	20–150 km	[300]
Hyogo-ken Nanbu (Kobe), Japan	7.2 (MJMA)	11/06/1996	DC potential, LF radio waves and MF and HF	223 Hz and 77.1 MHz and 1–20 kHz, 163 kHz	LF Omega transmitter and receiver	Visual, statistical analysis	<7 days	>100 km	[311]
Hyogo-ken Nanbu (Kobe), Japan	7.2 (MJMA)	11/06/1996	HF radio waves	10.2 kHz	LF Omega transmitter and receiver	Statistical analysis, modelling	2 days	70 km	[7]
Akita-ken Nairiku-nanbu, Japan	5.9	11/08/1996	LF and HF electric	10 kHz and 1 MHz	Vertical-dipole ground electrodes	Visual analysis and analysis of related parameters	6 days	<100 km	[311]
Chiba-ken Toho-oki, Japan	6.6	11/09/1996	Electric	10 kHz, 1 MHz	Vertical-dipole ground electrodes	Visual analysis and analysis of related parameters	3 days	320 and 430 km	[311]
Umbria–Marche, Italy	5.5	26/03/1998	LF radiowaves,	0.006 Hz	Radio wave vertical antenna		1.5 months	818 km	[312]
San Juan Bautista, California	5.1 ($M_W$)	12/08/1998	UHF magnetic	0.01–10 Hz	3-Component magnetic field inductor coils	Power spectrum analysis	2 h	3 km	[313]
Egio, Eratini, Greece	6.6 ($M_W$)	07/09/1999	LF electric and HF magnetic	41, 54 MHz and 3, 10 kHz	Electric dipole, magnetic loop antennas	Fractal analysis, block entropy	12–17 h	<300 km	[314]
Athens, Greece	5.9 ($M_W$)	07/09/1999	SES and LF electric and HF magnetic	1 Hz and 41, 54, 135 MHz and 3, 10 kHz	ULF, electric dipole and magnetic loop antennas	Fractal analysis, block entropy	<3 h	247 km	[90]
Athens, Greece	5.9 ($M_W$)	07/09/1999	HF magnetic	3, 10 kHz	Magnetic loop antennas	Delay times method, block entropy, spectral fractal analysis	12–17 h	247 km	[315]
Athens, Greece	5.9 ($M_W$)	07/09/1999	HF magnetic	3, 10 kHz	Magnetic loop antennas	Fractal analysis	12–17 h	247 km	[316]
Athens, Greece	5.9 ($M_W$)	07/09/1999	HF magnetic	3, 10 kHz	Magnetic loop antennas	Symbolic dynamics	12–17 h	247 km	[316]
Athens, Greece	5.9 ($M_W$)	07/09/1999	HF magnetic	3, 10 kHz, HF electric 41, 54 MHz	Electric dipole antennas, magnetic loop antennas	Wavelet power spectrum analysis	12–17 h	247 km	[304,305]
Athens, Greece	5.9 ($M_W$)	07/09/1999	HF magnetic	10 kHz	Electric dipole antennas, magnetic loop antennas	Block entropy	12–17 h	247 km	[259]
Athens, Greece	5.9 ($M_W$)	07/09/1999	HF magnetic	3, 10 kHz	Magnetic loop antennas	Block entropy	12–17 h	247 km	[316]
Athens, Greece	5.9 ($M_W$)	07/09/1999	HF magnetic	3, 10 kHz, electric 154 MHz	Electric dipole and magnetic loop antennas	Intermittent dynamics of critical fluctuations	20 h	247 km	[307]
Athens, Greece	5.9 ($M_W$)	07/09/1999	LF electric and HF magnetic	135 MHz and 3, 10 kHz,	Electric dipole and magnetic loop antennas	Intermittent dynamics of critical fluctuations	>3 h	247 km	[307]
Athens, Greece	5.9 ($M_W$)	07/09/1999	HF magnetic	10 kHz	Magnetic loop antennas	Tsallis entropy	12–17 h	247 km	[317]
Chi-Chi, Taiwan	7.6 ($M_W$)	20/09/1999	foF2 ionospheric		IPS-42 ionosonde	Visual observation	3–4 days	120 km	[318]
Chia-Yii, Taiwan	6.4 ($M_W$)	22/10/1999	foF2 ionospheric		IPS-42 ionosonde	Visual observation	b1–3 days	179 km	[318]
Izu-Penisula, Japan	6.4 (MJMA)	01/07/2000	ULF magnetic	0.001–1 Hz	3-Axis ring core-type fluxgate magnetometers	Fractal analysis with FFT, Higuchi, Bulgara–Klein methods	<1 month	80 km–1160 km	[249]
Izu-Penisula, Japan	6.4 (MJMA)	01/07/2000	ULF magnetic	0.001–1 Hz	3-Axis ring core-type fluxgate magnetometers	Fractal analysis with FFT, fractal dimension	<1 month	80 km–1160 km	[254]
Izu-Penisula, Japan	6.1 (MJMA)	09/07/2000	ULF magnetic	0.001–1 Hz	3-Axis ring core-type fluxgate magnetometers	Fractal analysis with FFT, Higuchi, Bulgara–Klein methods	<1 month	80 km–1160 km	[249]
Izu-Penisula, Japan	6.1 (MJMA)	09/07/2000	ULF magnetic	0.001–1 Hz	3-Axis ring core-type fluxgate magnetometers	Fractal analysis with FFT, fractal dimension	<1 month	80 km-1160 km	[254]
Izu-Penisula, Japan	6.3 (MJMA)	15/07/2000	ULF magnetic	0.001–1 Hz	3-Axis ring core-type fluxgate magnetometers	Fractal analysis with FFT, Higuchi, Bulgara–Klein methods	<1 month	80 km–1160 km	[249]
Izu-Penisula, Japan	6.3 (MJMA)	15/07/2000	ULF magnetic	0.001–1 Hz	3-Axis ring core-type fluxgate magnetometers	Fractal analysis with FFT, fractal dimension	<1 month	80 km–1160 km	[249]
Izu-Penisula, Japan	6.4 (MJMA)	18/08/2000	ULF magnetic	0.001–1 Hz	3-Axis ring core-type fluxgate magnetometers	Fractal analysis with FFT, Higuchi, Bulgara–Klein methods	<1 month	80 km–1160 km	[254]
Lefkas, Greece	5.9 ($M_W$)	14/06/2003	LF electric and HF magnetic	41, 54 MHz and, 3, 10 kHz	Electric dipole and magnetic loop antennas	Fractal analysis, block entropy	12–17 h	<300 km	[314]
Andaman, Sumatra, Indonesia	9.0 ($M_W$)	26/12/2004	ULF magnetic	1 Hz	3-Axis ring core-type, fluxgate magnetometers	Spectral density ratio analysis, transfer functions analysis, fractal dimension	<1.5 month	<750 km	[319]
Andaman, Sumatra, Indonesia	8.7	26/12/2004	ULF magnetic	1 Hz	CHAMP satellite vector magnetic antennas	Wavelet power spectrum analysis	2 h	700 km	[320]
Nias, Sumatra, Indonesia	8.7 ($M_W$)	28/03/2005	ULF magnetic	1 Hz	3-Axis ring core-type, fluxgate magnetometers	Spectral density ratio analysis, transfer functions analysis, fractal dimension	<1.5 month	<750 km	[319]
Nias, Sumatra, Indonesia	8.7 ($M_W$)	28/03/2005	ULF magnetic	1 Hz	CHAMP satellite vector magnetic antennas	Wavelet power spectrum analysis	2 h	700 km	[320]
Miyagi-ken oki, Japan	7.2 ($M_W$)	16/08/2005	Electric	49.5 MHz	Discon-type antenna from 25–1300 MHz	Multifractal Detrended Fluctuation Analysis	2–3 weeks, few days for Kunimi station	90–140 km	[240]
Mid Niigata prefecture	6.8 (MJMA)	16/08/2005	DC and ULF magnetic and HF electromagnetic	0.02–0.05 Hz and 40 kHz	3-Axis ring core-type fluxgate magnetometers, Discon-type antennas from 25–1300 MHz	Signal analysis with FFT	17–21, 5–7 days	<220 km	[321]
Greece	5.2 ($M_L$)	18/01/2007	SES	≤1 Hz	Electric and magnetic antennas	Natural time analysis	3 min	<150 km	[322]
Greece	5.8 ($M_L$)	03/02/2007	SES	≤1 Hz	Electric and magnetic antennas	Natural time analysis	22 min	<150 km	[322]
Vanuatu, Japan	7.1 (MJMA)	25/03/2007	TEC		DEMETER satellite	Statistical analysis	15 days		[323]
Honshu, Japan	6.7 (MJMA)	25/03/2007	TEC		DEMETER satellite	Statistical analysis	15 days		[323]
Lesvos, Greece	6.1 ($M_L$)	12/06/2007	LF electric and HF magnetic	41, 54 MHz and 3, 10 kHz	Electric dipole and magnetic loop antennas	DFA, power law	10–12 days	30 km	[230]
Wenchuan, China	8.0 ($M_s$)	12/05/2008	DC, ULF	≤1 Hz	Cr18Ni9C electrodes	Visual observations	3 days	1000 km	[324]
Greece	6.4 ($M_W$)	08/06/2008	SES	≤1 Hz	Electric and antennas	Natural time analysis		<30 km	[261]
L’Aquila, Italy	6.3	06/04/2009	LF electric and HF magnetic	41, 54 MHz and 3, 10 kHz	Electric dipole and magnetic loop antennas	Fractal analysis, block entropy, DFA, R/S analysis, Hurst analysis,	<3 h	816 km	[1,82]
Oran, Algeris	5.5 ($M_w$)	06/06/2008	Rinex, $F_2$ disturbances, TEC		Geodetic stations	Seismological, spectral analysis	Several days		[325]
Tokachi, Japan	8.0 (MsMA)	26/09/2003	SES	≤1 Hz	Electric antennas	Natural time analysis	1 month	$\delta$lat, $\delta$long < 30	[79]
Yutian, China	7.3 ($M_s$)	20/03/2008	TEC and ULF electric field data		Onboard DEMETER, Swarm and China’s seismo-electromagnetic satellites	Statistical, visual analysis	3 min–2 days		[326]
Lake Baikal, Siberia	6.3	27/08/2008	Electromagnetic signals from thunderstorms	VLF range	Single-point lightning direction finder-rangefinder	Visual observations	Hours		[327]
Indonesia	5.0	07/01/2009	Electromagnetic signals from thunderstorms	VLF range	Single-point lightning direction finder-rangefinder	Visual observations	7 days		[327]
Chichi-jima, Japan	7.8 (MJMA)	22/10/2010	SES	≤1 Hz	Electric antennas	Natural time analysis	1 month	$\delta$lat, $\delta$long < 30	[79]
Conception, Chile	8.8 ($M_W$)	27/02/2010	$N_m{f}_2$ ionospheric anomalies		FORMOSAT-3/COSMIC satellite	Kriging interpolation, global $N_m{f}_2$ map	5 h	epicentre area	[328]
Tohoku, Japan	9.0 (MJMA)	11/3/2011	SES	≤1 Hz	Electric antennas	Natural time analysis	1 month	$\delta$lat, $\delta$long < 30	[79]
Tohoku, Japan	9.0 (MJMA)	11/3/2011	GPS TEC		Modified single layer mapping function at the ionospheric pierce points at 350 km	GPS satellites (PRN 18, PRN26)	40–50 min	500–600 km	[329,330]
Tohoku, Japan	9.0 (MJMA)	11/03/2011	Ionospheric measurements	HF 3–25 MHz	Ionosonde detection network combined with Digisondes and COSMIC satellite	HF Doppler, planar ionospheric disturbances	6 h after	2000 km	[331]
Japan	6.0	14/03/2012	Electromagnetic signals from thunderstorms	VLF range	Single-point lightning direction finder-rangefinder	Visual observations	10 days	3000 km	[327]
India	5.6	25/04/2012	HF electric field	3.012 kHz	GPS terrestrial vertical antenna	Visual observations	1–13 days	2671 km	[332]
India	5.6	27/04/2012	HF electric field	3.012 kHz	GPS terrestrial vertical antenna	Visual observations	1–13 days	3284 km	[332]
Dholavira, India	5.1 ($M_w$)	20/06/2012	ULF magnetic and $K_p$,$D_{st}$ data	0.001–0.5 Hz	Digital fluxgate magnetometer	Visual and fractal dimensions	7 days	around, above epicentre	[333]
Yutian, China	6.3 ($M_s$)	12/08/2012	ULF electric field data, TEC	≤1 Hz	Onboard DEMETER, Swarm and China’s seismo-electromagnetic satellites	Statistical, visual analysis	10–20 days		[326]
India	5.9	22/07/2013	HF electric field	3.012 kHz	GPS terrestrial vertical antenna	Visual observations	1–13 days	2642 km	[332]
India	5.7	20/09/2013	HF electric field	3.012 kHz	GPS terrestrial vertical antenna	Visual observations	1–13 days	1905 km	[332]
India	5.7	02/10/2013	HF electric field	3.012 kHz	GPS terrestrial vertical antenna	Visual observations	1–13 days	2766 km	[332]
Yutian, China	7.3 ($M_s$)	12/02/2014	TEC and ULF electric field data		Onboard DEMETER, Swarm and China’s seismo-electromagnetic satellites	Statistical, visual analysis	Same days		[326]
Greece	6.9	24/05/2014	SES and geomagnetic signals	0.5–40 Hz and 0.0001–100 kHz	Mikhnevo GPO (seismometric, radiophysical, magnetometric, electrical) equipment				[334]
Ileia, Greece	4.4 ($M_L$)	30/08/2015	HF magnetic	3, 10 kHz	Magnetic loop antennas	Fractal analysis	3 days	24 km	[35]
Illapel, Chile	8.3 ($M_w$)	16/09/2015	Co-seismic ionospheric TEC	0.1–1 Hz	Global Navigation Satellite System	Wave perturbation ionosphere model with seismic source		1500 km	[335]
Ileia, Greece	4.5 ($M_L$)	12/12/2015	HF magnetic	3, 10 kHz	Magnetic loop antennas	Fractal analysis	3 days	24 km	[35]
Sumatra	7.8 ($M_w$)	02/03/2016	TEC	3.012 kHz	GPS terrestrial vertical antenna	3D tomography method	11–16 min after	$1^{\circ }$, 75 km	[336]
Afghanistan	6.6	10/04/2016	Seismic, geomagnetic and acoustic signals	0.5–40 Hz and 0.0001–100 kHz and ${10}^{-4}$–20 Hz	Mikhnevo observatory, LEMI-018 triaxial fluxgate magnetometer	Visual observations		2000–3000 km	[334]
Italy	6.6	30/06/2016	Seismic, geomagnetic and acoustic signals	0.5–40 Hz and 0.0001–100 kHz and ${10}^{-4}$–20 Hz	Mikhnevo observatory, LEMI-018 triaxial fluxgate magnetometer	Visual observations		2000–3000 km	[334]
Chiapas, Mexico	M8.2	06/07/2017	SES	≤1 Hz		Natural time analysis	Few hours		[64]
Greece	6.6	20/07/2017	Seismic, geomagnetic and acoustic signals	0.5–40 Hz and 0.0001–100 kHz and ${10}^{-4}$–20 Hz	Mikhnevo observatory, LEMI-018 triaxial fluxgate magnetometer	Visual observations		2000–3000 km	[334]
Mexican flat slab	M7.1	19/09/2017	SES	≤1 Hz		Natural time analysis	Several hours		[64]
Iraq	7.3	12/11/2017	Seismic, geomagnetic and acoustic signals	0.5–40 Hz and 0.0001 Hz–100 kHz and ${10}^{-4}$–20 Hz	Mikhnevo observatory, LEMI-018 triaxial fluxgate magnetometer	Visual observations		2000–3000 km	[334]
Ileia, Greece	4.5 ($M_L$)	07/05/2018	HF magnetic	3, 10 kHz	Magnetic loop antennas	Fractal analysis	3 days	24 km	[35]
Lombok, Indonesia	6.4	28/07/2018	Ne, Te and TEC	Onboard sensors	China’s seismo-electromagnetic satellites	$dTEC$, Statistical analysis	1–5 days	2000 km	[337]
Lombok, Indonesia	6.8	05/08/2018	Ne, Te and TEC	Onboard sensors	China’s seismo-electromagnetic satellites	$dTEC$, Statistical analysis	1–5 days	2000 km	[337]
Lombok, Indonesia	5.9	09/08/2018	Ne, Te data and TEC	Onboard sensors	China’s seismo-electromagnetic satellites	$dTEC$, Statistical analysis	1–5 days	2000 km	[337]
Lombok, Indonesia	6.9	19/08/2018	Ne, Te data and TEC	Onboard sensors	China’s seismo-electromagnetic satellites	$dTEC$, Statistical analysis	1–5 days	2000 km	[337]
Indonesia	7.5 ($M_w$)	28/09/2018	Physical properties of atmosphere and NeTe, ionospheric disturbances		China’s seismo-electromagnetic satellites	Seismological, climatological analysis	3.7, 6 months and 2.7 months	$3^{\circ }$	[338]
Zakynthos, Greece	6.6 ($M_L$)	25/10/2018	LF electric and HF magnetic	41, 54 MHz and 3, 10 kHz	Electric dipole and magnetic loop antennas	Fractal analysis, block entropy, DFA, R/S analysis, Hurst analysis	Post-activity	40 km	[34]
Ileia, Greece	4.3 ($M_L$)	04/02/2019	HF magnetic	3, 10 kHz	Magnetic loop antennas	Fractal analysis	3 days	24 km	[35]
Ridgecrest, Mexico	M7.1	06/072019	SES	≤1 Hz		Natural time analysis	Several hours		[64]
Indonesia	6.9 ($M_w$)	07/07/2019	VLF	48.83–366.21 Hz	Electric field detector of China’s seismo-electromagnetic satellites	Electric field PSD	Before and after	near the epicentre	[14]
Indonesia	7.2 ($M_w$)	14/07/2019	VLF	48.83–366.21 Hz	Electric field detector of China’s seismo-electromagnetic satellites	Electric field PSD	Before and after	near the epicentre	[14]
Laiwui, Indonesia	7.2 ($M_w$)	14/07/2019	TEC, plasma, global ionospheric map		China’s seismo-electromagnetic satellites	Cross-validation analysis and moving-mean method	1, 3, 8 days		[335]
Jiashi, China	6.4 ($M_s$)	19/01/2020	Electron density and rock temperature		Zhangheng-1 electromagnetic satellite	Visual observations	15 days	150 km	[339]
Yutian, China	6.5 ($M_s$)	25/06/2020	ULF, TEC, Global ionospheric Map	≤1 Hz	Onboard DEMETER, Swarm and China’s seismo-electromagnetic satellites	Statistical, visual analysis	Same days		[326]
Turkey	7.8 ($M_w$)	06/02/2023	TEC		Global Navigation Satellite System, ionosondes	Statistical, visual analysis	22–25 min after	750 km	[16]
Turkey	7.5 ($M_w$)	06/02/2023	TEC		Global Navigation Satellite System, ionosondes	Statistical, visual analysis	22–25 min after	750 km	[16]

presents a collection of papers for electromagnetic precursors.

Table 2. Papers of radon precursors. The papers are presented in chronological order from the oldest to the newest. The precursory time also includes the aftershock data presented in some papers. RA stands for the relative amplitude of the radon anomalies and AD for the anomaly duration. SSNTDs stands for solid-state nuclear track detectors. ED is the effective-sensitive distance between the monitoring site and the epicentre of the earthquake. Blank cells indicate there is no information available in the reference(s). Russian Federation is used as the successor state of the former USSR.

Location	Magnitude	Date(s)	RA	AD (days)	Instrumentation	Methodology	Precursory Time	ED	Reference
Pohai Bay, China	7.4	18/07/1969	60%	170 days	Instruments of Kutzan station for radon in water	Visual observations		200 km	[188]
Szechwan Luhuo, China	7.9	06/02/1973	120%	9 days	Instruments of Tangku station for radon in water	Visual observations		170 km	[186,188]
Markansu, Russian Federation	7.3	04/02/1975	38% and 17%	270 days and 50 days	Instruments of Alma-Ata station for radon in water	Visual observations		530 km	[188]
Liaoning, Haicheng, China	7.3	04/02/1975	38% and 17%	270 days and 50 days	Instruments of Tangangzi station for radon in soil	Visual observations		50 km	[188,340]
Liaoning, Haicheng, China	7.3	04/02/1975	10%	1 day	Instruments of Liaoyang station for radon in soil	Visual observations		85 km	[188,341]
Gazli, Russian Federation	7.3	17/05/1976	220%	4 days	Instruments of Tashkent station for radon in water	Visual observations		470 km	[188]
Yunnan Lungling, China	7.5	29/05/1976	20%	510 days	Instruments of Lungling station for radon in soil	Visual observations		190 km	[186,188]
Yunnan Lungling, China	7.5	29/05/1976	8%	160 days	Instruments of Erhyuan station for radon in soil	Visual observations		470 km	[186,188]
Szechwan Songpan Pingwu, China	7.2	16/08/1976	29%	480 days	Instruments of Erhyuan for radon in soil	Visual observations		40 km	[186,188]
Szechwan Songpan Pingwu, China	7.2	16/08/1976	70%	7 days	Instruments of Kutzan station for radon in soil	Visual observations		320 km	[188,341]
Hopeh Tangshan, China	7.8	27/07/1976	30%	5 days	Instruments of Tangshan station for radon in water	Visual observations		5 km	[188,342]
Hopeh Tangshan, China	7.8	27/07/1976	50%	15 days	Instruments of Antze station for radon in water	Visual observations		100 km	[188,342]
Isferi Batnen, Russian Federation	6.6	31/01/1977	−30%	60 days	Instruments of Tashkent station for radon in water	Visual observations		190 km	[188]
Hopeh Chienan, China	6.0	04/03/1977	70%	3 days	Instruments of Peking station for radon in water	Visual observations		200 km	[188,341]
Hopeh Lutai, China	6.7	12/03/1977	30%	1 day	Instruments of Tungchao station for radon in water	Visual observations		115 km	[188,341]
Isferi Batnen, Russian Federation	6.6	24/03/1977	−20%	125 days	Instruments of H-O-Garm station for radon in water	Visual observations		200 km	[188]
Alma-Ata, Russian Federation	7.1	04/02/1978	32%	50 days	Instruments of Alma-Ata station for radon in water	Visual observations		65 km	[188]
Zaslai, Russian Federation	6.7	01/11/1978	−30%	470 days	Instruments of Obi-Garm station for radon in water	Visual observations		270 km	[188]
Zaslai, Russian Federation	6.7	01/11/1978	−40%	470 days	Instruments of Yavros station for radon in water	Visual observations		300 km	[188]
Izu-Oshima, Japan	6.8	14/01/1978	7%	230 days	Instruments of SKE-1 station for radon in water	Visual observations		25 km	[186,188]
Izu-Oshima, Japan	6.8	14/01/1978	−8%	7 days	Instruments of SKE-1 station for radon in water	Visual observations		25 km	[186,188]
Imperial valley, California, USA	6.6	15/10/1979	400%	116 days and 50 days	Instruments of KPAS station	Radon in water		335 km	[187,188]
Irpinia, Italy	6.5	23/11/1980	170%	5–6 months	Instruments of Rieti station for radon in groundwater	Visual observations	4 months	150 km	[343]
Japan	7.9	06/03/1984	few days		Instruments for radon in groundwater	Bayesian statistics, ±2$\sigma$	1 week	1000 km	[344]
Japan	6.7	06/02/1987	few days	4	Instruments for radon in groundwater	Bayesian statistics, ±2$\sigma$	3 days	130 km	[344]
Equador	6.9	06/03/1987	230%	30 days	Radon in soil, SSNTDs	Visual observations	50 days	200 km	[345]
Uttarkashi, India	7.0 ($M_s$)	20/10/1991	180%	7 days	Radon in soil, SSNTDs	Visual observations	1 week	450, 330 km	[346,347]
Mindoro, Philippines	7.1	11/04//1994	600%	7 days	BARASOLVDG	Visual observations	22 days	48 km	[348]
Kobe, Japan	7.2	1/17/1995	−2%	4 months	Radon in atmosphere, flow ionisation chamber at 18 m	Daily min data analysis	4 to 0 months	130 km	[29,214,349]
Chamoli, India	6.5 ($M_s$)	29/03/1999	200%	2 days	Radon in soil, water with emanometric technique	±2$\sigma$	1–7 days	393 km	[347]
Hiwacho-Mitsugaichi, Shobara, Japan	7.3 (MJMA)	06/10/2000	16–20%	>6 months	Gas flow ionisation chamber	Residual analysis		207 km	[215]
Scotia sea, Antarctica	7.5 ($M_s$)	04/08/2003	400–700%	16 days	CR-39, TASTRAK	Visual, power law	6	1176 km	[350]
Chengkung, Taiwan	6.8	10/12/2003	−13%	6 months	Radon in water, liquid scintillation counter, wells 167–187 m deep	30 km	65 days	20 km	[193]
Yura, Hidaka, Japan	7.4 (MJMA)	05/10/2004	16–20%	>6 months	Gas flow ionisation chamber	Residual analysis		22 km	[215]
Indonesia	9.1	26/12/2004	60%	4–6 days	Radon and progeny in gases from thermal springs at Bakreswar, India	±2$\sigma$, visual observations		2275 km	[351]
Middle Kurils, Simushir Island, Kamchatka Peninsula	8.1($M_w$)	20/04/2006	33–35%		Gas-discharge counter for radon progeny	Visual observations	8 months–3 years	800 km	[153]
Olutorsk, Kamchatka Peninsula	7.6 ($M_w$)	/20/04/2006	33–35%	33–35%	Gas-discharge counter for radon progeny	Visual observations	8 months–3 years	1035 km	[153]
Middle Kurils Kamchatka Peninsula Simushir Island, Pacific Ocean	8.3 ($M_w$)	13/01/2007	33–35%		Gas-discharge counter for radon progeny	Visual observations	8 months–3 years	800 km	[153]
Wenchuan, China	8 ($M_s$)	12/05/2008	10 times the baseline	12 days	SD-3 A, automatic radon instrument, Guzan station	Statistical analysis		155 km	[204]
Wenchuan, China	8 ($M_s$)	12/05/2008	5 times the baseline	Scattered days	FD-125, ZnS(Ag)	Sliding window power law, DFA, fractal dimension, 13-method combination analysis	1–2 months	150–500 km	[36]
Kato Achaia, Peloponnese, Greece	6.5 ($M_L$)	06/08/2008	20 times the baseline	12 h	Alpha GUARD, CR-39, radon in in soil	Sliding window power law, statistics, outliers	2 months	40 km	[2]
Kato Achaia, Peloponnese, Greece	6.5 ($M_L$)	06/08/2008	20 times the baseline	12 h	Alpha GUARD radon in in soil	Sliding window power law, DFA, spectrogram, scalogram	2 months	40 km	[23]
Kato Achaia, Peloponnese, Greece	6.5 ($M_L$)	06/08/2008	20 times the baseline	12 h	Alpha GUARD radon in in soil	Sliding window fractal dimension analysis, Hurst exponents	2 months	40 km	[23]
Kato Achaia, Peloponnese, Greece	6.5 ($M_L$)	06/08/2008	20 times the baseline	12 h	Alpha GUARD radon in in soil	Sliding window $R/S$, DFA and block entropy analysis, R-L, variogram methods, fractal dimensions	2 months	40 km	[21]
Aegean Sea, Lesvos area, Greece	5.0 ($M_L$)	19/03/2008	20 times the baseline	1 h	Alpha GUARD radon in soil	Sliding window $R/S$, DFA and block entropy analysis, R-L, variogram methods, fractal dimensions	3 months	40–70 km	[21]
Tohoku, Japan	9.0 (MJMA)	11/03/2011	80–160 times the baseline	>16 days	Radon, thoron instrumentation at Seongryu Cave	Statistical, visual analysis	1 month		[204]
PhekN agaland, India	5.8	29/07/2012	2–3 times the baseline	1 month	LR-115 in soil	±2$\sigma$, visual observations	16–31 days	224 km	[352]
Myanmar, India	6.0	29/07/2012	2–3 times the baseline	1 month	LR-115 in soil	±2$\sigma$, visual observations	16–31 days	132 km	[352]
Awaji Island, Japan	6.7 (MJMA)	13/04/2013	16–20%	>6 months	Gas flow ionisation chamber	Residual analysis		44 km	[215]
Luhsan, Cina	7 ($M_s$)	20/04/2013	10 times the baseline	20 days	SD-3 A, automatic radon instrument, Guzan station	Statistical analysis		82 km	[204]
Gansu, China	6.6 ($M_s$)	22/07/2013	10–20%	2 months	FD-125 instrument, radon in groundwater	Monofractal, multifractal DFA		688 km	[179]
Evia Island, Greece	5.0 ($M_L$)	15/11/2014	−5 times the baseline	10 min	VDG BARACOL, radon in soil	Sliding window $R/S$, DFA, scalograms	10–12 days	100 km	[32]
Nepal	7.8	25/04/2015	4 times the baseline	15 days	LR-115 in soil	±2$\sigma$, visual observations	5 days	722 km	[353]
West Bengal, India	7.8	26/04/2015	3.5 times the baseline	15 days	LR-115 in soil	±2$\sigma$, visual observations	6 days	612 km	[353]
Kalamei, Nepal	7.8	12/05/2015	3 times baseline	15 days	LR-115 in soil	±2$\sigma$, visual observations	5 days	618 km	[353]
Lesvos Island, Greece	4.1 ($M_L$)	10/09/2015	8–20 times the baseline		Alpha GUARD radon in soil	Sliding window $R/S$, DFA, scalograms		50 km	[235]
Lesvos Island, Greece	4.6 ($M_L$)	26/10/2015	8–20 times the baseline		Alpha GUARD radon in soil	Sliding window $R/S$, DFA, scalograms		50 km	[235]
Zhupanovo, Kamchatka Peninsula	7.2 ($M_w$)	30/01/2016	33–35%		Gas-discharge counter for radon progeny	Visual observations	8 months–3 years	110 km	[153]
Jiuzhaigou	7 ($M_s$)	08/08/2017	±3 times	>2 months	SD-3 A, automatic radon instrument, Songpan station	Statistical analysis		67 km	[204]
Uglovoye Podnyatiye, Kamchatka Peninsula	7.3 ($M_w$)	20/12/2018	33–35%		Gas-discharge counter for radon progeny	Visual observations	8 months–3 years	490 km	[153]
North Kurils, Kamchatka Peninsula	7.5 ($M_w$)	25/03/2020	33–35%		Gas-discharge counter for radon progeny	Visual observations	8 months–3 years	449 km	[153]

shows the paper collection for radon precursors. The papers in both tables are presented chronologically and relatively old events are also included. To avoid unessential records, historical earthquakes are limited to very strong and extremely strong ones. Although the knowledge and methodologies have evolved, the techniques addressing older earthquakes are definitely not new, but their publication time is contemporary with the investigated earthquakes. The historical electromagnetic precursors also include the great papers that pioneered and were breakthroughs in seismic analysis. Especially for the radon precursors, the table also includes the traditional treatment of radon anomalies which refers to the statistical dislocation of radon’s concentration with respect to the corresponding baseline values and the duration of the anomaly, since both are still used nowadays. Modern methods include in both cases fractal behaviour and self-organisation analysis.

The collection of papers is the most significant part of this review, because it gathers the knowledge and may assist the related research. Due to the great number of papers, a special presentation approach was selected according to which every row presents all the data of each earthquake and the maximum available information. Since every paper is special and the published information is not uniform, a variety exists in the data of the papers presented in both Tables as column entries. Because of this, effort was put to present specific information as well. All papers were accessed from the site of each journal and the available information was downloaded as a BibTEX file, or converted to BIB format from the corresponding RIS record of each journal. Digital Object Identification (doi) data were also searched and inserted wherever available.

In the next pages both tables are given. Table 1 and Table 2 support the view that the majority of publications are based on visual observations of collected data and a subsequent statistical analysis. This can be explained by the fact that it is very difficult to collect data by several nearby stations. As mentioned by Cicerone et al. [6], it is a serendipitous finding to have a strong earthquake and a station which collects data during the seismic rupture and is installed in the broader epicentral area. As mentioned in several publications and expressed collectively in Eftaxias [175], there is no one-to-one correspondence between an earthquake occurrence and an anomaly detection. Moreover, even the most advanced methods fail to deliver a very effective forecasting of earthquakes. These facts complicate the analysis even today (2024).

The most advantageous methods seem nowadays to reflect the fractal and self-organisation nature of the rupturing crust of the Earth during the preparation of earthquakes. A very robust method is the natural time analysis, which produces promising results. Satellites are now numerous and can be accessed conveniently. This provides new insights into the related research. Remote sensing and SAR techniques are very powerful as well, not neglecting the great number of installed stations worldwide. These new tools boost modern approaches which are also multi-facet and necessitate the collaboration of different groups.

Historically, radon gas claims the majority of publications in relation to earthquakes with many radon papers suggesting associations with very strong earthquakes. Nowadays, there is a balance between radon and electromagnetic precursors, with the latter providing more options due to the different frequency bands and the remote sensing and satellite methods. Radon has also led to new approaches, and therefore, both precursors are very significant. In fact electromagnetic and radon precursors seem to be the subject of many papers to date.

The collaboration between scientists and the use of multilevel approaches with different methodologies will be the key point for seismic research in the following years. This research is ongoing and in a continuous search for credible and powerful precursors.

8. Conclusions

This paper presented a review on electromagnetic and radon precursors for earthquake forecasting. The electromagnetic precursors emerge in diverse frequency bands ranging from ultra-low to very high frequencies. Nowadays, electromagnetic data are collected from satellites, whereas remote sensing techniques are increasingly used as well. Within the electromagnetic spectrum, TEC measurements and the modern approach of SAR studies are also found. Various investigators are still working independently; nevertheless, there is great space for collaborations. The traditional approach for earthquake prediction is still recordings from ground stations, with the precursors of the ULF range having the greater history and potential. MHz and kHz frequencies provide very good estimations as well. On the other hand, radon precursors are those with the oldest usage. Many great earthquakes have been studied with the help of radon stations worldwide. Radon is easily detected and may travel far due to its inert nature. For this reason, it is suitable for forecasting earthquakes occurring at relatively long distances.

The majority of the reported precursory anomalies have been and still are visually observed. Several statistical approaches have been utilised in the papers. Especially for radon, the ±2$\sigma$ criterion is the one most frequently used. In recent years, advanced methods have been published and used in several new publications. Special mention is given to the modern approach of natural time which has great potential and many future earthquakes to be applied to. Power-law as well as monofractal and multifractal Detrended Fluctuation Analysis have been used in both electromagnetic and radon precursors. Considerable attention has been given to $R/S$ analysis, fractal dimension analysis and Hurst exponents. Block entropy and several entropy measures have been used as well. A combinational analysis between different monofractal methods has been used with success. All these modern methods attempt to outline the fractal and self-organised critical features of the fracturing parts of the Earth’s crust during the preparation of earthquakes. Much research needs to be carried, and new approaches are still in demand.

Several models have been proposed for the interpretation of the collected precursory data. The LAIC model has been in great use by many papers. The theory of asperities has been employed both in electromagnetic and radon precursors. In radon research, other models have also been utilised. Since each earthquake is a special event, it is difficult to find a universal model which covers all aspects of the research outcomes. The main problem is that many precursors have been characterised as such, after the occurrence of the earthquakes. This is a disadvantage that will be overcome as the pertinent research progresses. There are papers that forecast earthquakes prior to their occurrence, and this their most distinguishing feature. The work of researchers from different sub-disciplines of electromagnetic and radon precursors will hopefully provide better forecasting results in the near future.