Electromagnetic and Radon Earthquake Precursors

1. Introduction

Earthquakes, volcanic eruptions and tsunamis are all inevitable disastrous phenomena. Not only that they are 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 to more than fifty years. There is evidence that pre-seismic electromagnetic radiation or radon concentration observations can be utilised for forecasting, taking into account specific measurable features of the associated earthquake process. More specifically, if such an observation takes place near the geological rupture, some measurable precursory activity prior to the seismic event can be expected.

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 [1]. Seismic forecasting usually falls into four categories [2]: 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 [3] 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. The separation into several stages is determined by the features of the processes that generate a massive earthquake and the needs for earthquake preparedness-which include a range of safety procedures for each level of forecast [1]. The reader should note, that there is no rarely any direct correlation between abnormalities in the measurements and earthquake occurrences,especially in short-term forecasting [4,5]. In seismic-prone countries, short-term earthquake pre-warning in a time window of weeks, days, or hours is deemed as most important, although being significantly more difficult than the long-term forecasting. The science of short-term earthquake forecasting is the study of short-term precursory activity occurring through systematic observations of physical quantities taking place near and before earthquake occurrences and can be further supported by serendipitous findings in observations not purposed for earthquake monitoring but are nonetheless acquired near the earthquake location [6]. Abnormalities in electromagnetic fields, anomalous variations of radon concentration in soil, groundwater, surface water and atmosphere, erratic gas emissions, uneven surface distortions caused by pressure differentials, irregular adjustments to ionospheric parameters, ionospheric perturbations, anomalies detected in satellite devices and other remote sensory devices and excess Total Electron Content (TEC) are among these physical quantities [6].

Observations of pre-seismic electromagnetic disturbances (of the Radio Frequency-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 [3,6,8,9,10,11] and references therein) pre-seismic electromagnetic emissions occur in a wide frequency range from frequencies well below 10 $Hz$ (Ultra Low Frequencies-ULF), between some $kHz$ range up to several $MHz$ (altogether characterised hereafter as High Frequencies-HF) and between 100 $MHz$ up to 300 $MHz$ (Very High Frequencies-VHF). The research originated back in the 1970s where the first successful seismic forecast was reported for an earthquake of magnitude 2.6 occurring on August 3, 1973, near Blue Mountain Lake, New York [12]. Following this, the M=7.4 Heicheng China earthquake of February 4, 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 hours 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]. Seismologists’s research has recently been focused on short-term forecast rather than long-term forecast [13]. The pre-seismic electromagnetic observations and abnormality recordings have been documented by several study teams throughout the globe as precursors of earthquakes. The EM 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 premonitory signs before to earthquakes [5,6,8,17,18,19,20,21,22,23,24,25,25]. This is due to the fact that radon may travel great distances from the host rocks where it is created [26] and can be detected at very low levels [27]. Anomalous radon concentration variations in soil gas, groundwater and atmospheremay be observed prior to earthquakes [6,17,19,24,28,29]. Before earthquakes, anomalous radon fluctuations are addressed in soil gas, groundwater, atmosphere, and thermal spas [6,17,19,24,28,29]. The time-series features, such as the range, length, and number of radon anomalies, as well as the precursory time and epicentral distance, vary greatly [6,17,30,31]. However, the amounts of 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 usually screened for atmospheric parameter influences [6,17,21,22,23]. The majority of the associations between radon and earthquakes are based on events of small and intermediate magnitudes. Large magnitudes earthquakes associations with radon observations also exist [5,32,33,34,35,36].

Ionospheric studies, satellite measurements and remote sensing devices have gained significant recent international interest in earthquake precursory investigations, after the realisation 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 [38,39,40,41,42]. Other researchers have studied the lower ionosphere extensively in relation to earthquakes by using different techniques to identify the precursory characteristics of earthquakes, as well as, the perturbation on the upper and lower regions [43,44,45,46,47]

Despite the massive scientific efforts, the causes of earthquake generation remain unknown. One important element is our inadequate understanding of the fracturing mechanisms of the crust [4,6,8,18,19,20,21,22,48,49,50,51,52,53,54,55,56]. Given that the fracture of heterogeneous materials is still not adequately described, despite a significant recent effort at the laboratory, theoretical, and numerical level [4], it is clear why the explanation of the genesis of earthquakes is still limited [4,6,18,19,20,21,22,48,49,50,51,52,53,54,55,56]. That explains why the related research is very significant, while the questions are far from being solved. Moreover, each earthquake is unique and has a wide distribution. Thus, not only is the understanding of the complex related processes mandatory, but is also sedulous to identify credible earthquake precursors and meticulous work is needed for that [4]. This is complicated further when considering that the successful selection of an earthquake precursor might probably need an acceptable physical model to explain its existence [6]. Eftaxias et al. [4] stated that before the final catastrophe, several geological, geochemical, hydrological, and environmental parameters generate different features and scales that are associated and describe the earthquake’ss precursory processes.

2. Electromagnetic Precursors

2.1. ULF Emissions

In 1964, seismogenic electromagnetic emissions with frequencies lower than 10 $Hz$ were first seen [57]. It has been found that variations in ground electric potential, ULF electromagnetic waves in the atmosphere, and other known phenomena occur prior to earthquakes [8,9,10,58,59,60,61,62]. Monitoring of ULF emissions directly recorded from the lithosphere is one of the several widely used seismo-electromagnetic methods. This is because ULF (f=0.01 $Hz$ – 10 $Hz$) has 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 [3]. Although most ULF precursors are electric, nowadays researchers study also magnetic ULF precursors [6,8,9]. It is noteworthy, nevertheless, that some contentious claims on ULF-range signals associated with earthquakes have also been reported[9].

The VAN method (from Varotsos, Alexopoulos, Nomicos), for ULF emissions has a long history 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 properties of SES i s selectivity [8,9,10]. That reflects the fact that the SES choose preferable paths and, consequently, a ULF station is sensitive to SES from certain seismic areas only, namely from some specific focal areas. The map showing these potential areas is called selectivity map of a station and has to be determined in order that the station produces useful data. Due to selectivity, SES can be detected even from hundred kilometres away of the epicentre. By installing proper dipoles in a cross of 50 m, 100 m, 200 m and, preferably also, 1000 m, magnetotelluric variations are discriminated from anthropogenic disturbances. The VAN method has successfully predicted events within a precursory window of some days or weeks both in Greece [61,62] and Japan [3,10]. Nowadays the SES ULF signals (as other signals as well) are incorporated in the modern method of Natural Time (see Section 5) which increases the successful forecasting of several earthquakes [e.g., [64], and references therein]. The discussion on the VAN method has divided the scientific community in those supporting it [10] it 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 all successfully predicted using ULF. Using the cumulative (daily sum) of the local energy of the earthquakes weighted by the squared distance from the measurement station method, which was suggested by Hattori et al. [68] and Hattori et al. [69], Han et al. [70] reported an increased probability of ULF magnetic anomalies 1-2 weeks before medium and strong shallow earthquakes, confirming previous findings published by Hattori et al. [69] and emphasizing that the perturbations are better associated with stronger and closer earthquakes. From 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 are reported by Hirano and Hattori [80] and Ouyang et al. [81]. ULF magnetic field emissions at are continuously measured in Agra station in India with the help of 3-component search coil magnetometers with promising forecasting results[63]. Two large earthquakes of magnitudes M=7.4 andM=6.8 occurred in Pakistan have been successfully predicted within 16 days with fractal methods (see section Section 5). ULF geomagnetic data from the Panagjurishte and Surlari stations in Romania have been successfully utilised for the forecast of a $M_w$=6.4 earthquake occurred at the coastal zone of Albania on 26 November 2019 [79].

2.2. HF Emissions

In the range between a few $kHz$ to several $MHz$, a number of HF emission disruptions have been reported prior to earthquakes. [4,23,34,35,82,83,84,85,86,87]. As stated by Hayakawa and Hobara [3], the two methods used to detect the seismic precursors are the direct measurement of the electromagnetic emissions radiated from the hypocentre of earthquakes in the lithosphere, or the indirect detection of propagation anomalous 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. [4,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 events at different scales and features will occur prior to an earthquake,as it is an abrupt mechanical breakdown in the heterogeneous earth’s crust and thus, the multiplexed operations that occur may be the initial source of numerous electromagnetic precursors of the widespread earth ’s crust collapse [4,82,86]. These are transient phenomena and a way to examine them is by analysing the observed pre-seismic time series, however, including sequences of discrete, brief time periods.The goal is to identify a clear shift in dynamical properties as the catastrophic event approaches near. In order to develop a quantitative identification of electromagnetic precursors, several mathematical concepts are utilised (see section , so as to set the criteria, to isolate detected anomalous electromagnetic emissions versus the noisy abnormal statistical patterns [84,85,86,88,89,90,91].

Several publications [92,93,94,95,96] suggest that the high persistency and organisation in a launched electromagnetic anomaly 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 a feedback mechanism similar to this might be a sign of the earthquake fracture process.Naturally,there is no study that can establish the high precursory value of a particular abnormality on its own There is still much to be done to comprehensively address the HF electromagnetic precursors. That is also valid for the ULF precursors also. Separating two events that happened at different times, like an earthquake and its potential HF electromagnetic precursor, is a challenging task. It is still to be determined if alternative methods may provide more data that would enable one to acknowledge the seismogenic source of the detected HF electromagnetic abnormalities and connect them to a pivotal phase of earthquake production.

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

2.2.1. VHF Emissions

VHF have been also employed in the search of electromagnetic earthquake precursors. According to Pullinets [95], one of the two authors of the LAIC model [37], in actuality, the LAIC is a complex system made up of subsystem interactions and a synergy of several processes, one of which is the VHF electromagnetic emission frequency band, which further functions by altering the characteristics of the atmosphere and ionosphere. Although some scientists dispute with the precursory usability of VHF emissions (e.g., [8], and references therein), the scientific interest is stimulated in the recent years on this subject. For example, Sorokin et al. [96] report a a full-fledged theoretical physical model, according to which, the over-horizon propagation of pulsed VHF radiation, can be explained. as also, the origin of such seismic related phenomena in link to the generation in the troposphere, the thermal effects and associated IR emissions, as well as, the modification of plasma distribution in the D, E and F layers of the ionosphere. Ouzunov et al. [97] 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 natural amplification of the signal’s strength days or hours before the seismic occurrences, even distant from the observation zones, such as the M=5.6 earthquake of May 22, 2012, in Bulgaria, M=5.1 earthquake of August 12, 2018, in Albania, the M=4.1 earthquake of August 2, 2018, in Southern Bulgaria and the M= 5.5 earthquake of October 28, 2018, in Romania. A VHF early warning system is utilised among other systems in Mexico [98]. Moriya et al. [99], on the basis of a designed a data-collection system, report several anomalous VHF-band radio-wave propagation prior to earthquakes, with most significant, the Tokachi-oki earthquake ($M_j$= 8.0, Mj, $M_j$ a magnitude defined by the Japan Meteorological Agency) on 2003 September 26 and the southern Rumoi sub-prefecture earthquake ($M_j$= 6.1) on 2004 December 14. Devi et al. [100] 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 more than 1,000 $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. Regarding VHF or higher frequencies that are pertinent to observations in radio astronomy. According to Erickson [101], anthropogenic electromagnetic emissions are primarily caused by mobile communications, car ignition systems, industrial equipment, and radio and television broadcasting stations. Eftaxias et al. [102] 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 and the degree of geotectonic heterogeneity within the focal zone.

2.3. Remote Sensing and Satellite Techniques

The application of space-borne remote sensing has grown in popularity and effectiveness within the field of natural disasters [103]. Improved quality data with repeated spatio-temporal coverage covering large areas in rough geomorphological and geological conditions can be obtained through the development of geospatial technologies and advanced data processing [104,105,106]. The post-disaster visualisation of remote sensing images helps in knowledge production, emergency intervention thinking and decision-making during the earthquakes [107]. In fact, the seismo-electromagnetic research has entered a new phase with the development of the remote sensing tools. This is because they make it possible to monitor a number of locations throughout the globe, including various 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 as well. Co-seismic effects occur as around 100 m-long earth cracks, which are followed by landslides, lateral spreading and changes to urban and suburban areas. For instance, landslides and substantial lateral spreading were noted following the earthquakes in Kashmir in 2005 and Mirpur in 2019 [108,109]. Human casualties from earthquakes are brought on by landslides, which also significantly alter agriculture and the food supply chain. These overall structural, stratigraphic and hydrogeological seismically-generated side-effects are significant features in remote sensing and satellite studies. As another example, the Landsat satellite imagery has been used to study the dynamic relationship between observed seismicity and lineament density[110]. 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 such hazard co-seismic mapping.

The Synthetic Aperture Radar (SAR) remote sensing techniques are among the best candidates for mapping co-seismic changes Interferometric SAR (InSAR) is one of the most powerful remote sensing techniques of the SAR family, that has been used to detect several surface deformations over large areas with high accuracy [111]. InSAR-based remote sensing methods allow low-speed surface deformations to be detected over vast areas with centimeter to millimeter precision [112]. The permanent Scatterers InSAR method is also accepted as a robust technique for mapping co-seismic deformation and in-field conditions, as well as, movements of urban infrastructures [113,114]. Space-borne remote sensing techniques are less effective towards this directions because they cannot provide complete information on the near-surface features produced as a consequence of an earthquake with the potential to damage the built environment.

Like the non-destructive near-surface geophysical remote sensing methods, the ground penetrating radar (GPR) method has been applied to shallow subsurface seismic investigations due to its high-resolution, time and cost-effective nature [115,116,117]. GPR has gained popularity in studies related to the detection of faults and fracture networks [118], slope instabilities [117] and landslides [119]. The GPR is one of the reliably accurate mapping tools to study a single site and imaging of a localised subsurface deformation but difficult to perform such surveys over an extensive earthquake Karst depressions-landslide affected area to detect the near-surface target features. Among the aforementioned geohazards, few studies focusing on co-seismic liquefaction and related ground failure have been conducted using field GPR measurements [120,121].

2.3.1. TEC

Total Electron Content, or TEC, is the electron density of a 1 $m^2$ cylinder that is vertically stacked from a ground point to the ionosphere [122,123]. 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 [123]. The corresponding data is accessible through a number of repositories and URLs [123,124,125,126,127] via the Ionosphere Exchange (IONEX) data file structure [128].

To investigate seismically generated TEC fluctuations in the ionosphere, researchers have used a variety of schemes and approaches [39,123,129,130,131,132,133,134]. There has been much discussion about 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. [130], 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 perturb TEC alterations within it, hence modifying the seismo-electromagnetic environments surrounding the epicentre. A few days prior to the May 12, 2008 $M_w$=7.9 Wenchuan earthquake, Zhao et al. [135], Liu et al. [39], and Pulinets and Ouzounov [37] report that ionospheric GPS TEC enhancement and, or, reduction anomalies simultaneously appear 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 [122,123,125]. The measurements include variability in the critical frequency of the $F_2$ layer, $f_o{F}_2$ and TEC [125,136]. 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 [125]. The efficiency of the impact of earthquakes on the ionosphere is growing with earthquake magnitude and depth representing relative density TEC anomalies within area of 1000 $km$ radius around the earthquake’s hypocentre [126]. Gulaeva and Arikan [126] suggest that the positive TEC storm anomalies are twice as much as those of non-storm values and that this observation supports dominant post-earthquake TEC enhancement with ionosphere peak decreasing during 12 h for daytime but growing by night-time during 6 h after the earthquake and followed by gradual recovery afterwards.

According to Sorokin et al. [96], there are two possible causes of 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 [137,138,139]. For instance, TEC increased during the 2014-2015 high solar radiation cycle, which was brought on by high-energy solar particles interacting with the earth’s ionosphere resulted in TEC shifting [122]. Therefore all these parameters should be taken into account when studying TEC ionospheric variations.

3. Radon Precursors

3.1. Radon Properties

Radon (${222}^{Rn}$) is a natural radioactive noble gas. It is produced when radium (${226}^{Ra}$) decays. According to Nazaroff and Nero [26] there are thirty-nine known isotopes of radon, ranging from ${193}^{Rn}$ to ${231}^{Rn}$. Radon has a half-life of 3.823 days and is the most stable isotope. ${220}^{Rn}$, or else, thoron, has a half-life of 54.5 seconds, ${220}^{Rn}$. Due to its short half-life, thoron decays rapidly and because of this it is often detected at low concentrations. That depends however on the concentration of its parent nucleus (${224}^{Ra}$), especially in comparison to that of ${226}^{Ra}$. Radon is primarily responsible for the radioactivity present in the atmosphere at sea level [140].

Radon emissions mostly originate from soil [26]. About 10% of the radon that is diluted in soil gets released into the atmosphere [140]. In addition to soil, radon may be found in surface and underground waters, as well as fragmented rock [140,141]. 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 [141]. Once there, either convection or molecular diffusion advection can cause a macroscopic transport [26]. Interconnected pores and water aquifers allow this movement to appear [141]. Radon dissolves into the water present in the pores of soil and rock and is carried away by it [26]. The migration of radon in soil and fragmented rock is implemented by all fluids present there, enclosed air included [17,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 significantly to the effective dose equivalent and makes up over half of the population’s exposure to natural sources and the leading natural cause of lung cancer [142,143,144]. Due to this, radon is a subject of extensive research worldwide [142,143,144,145,146,147,148,149,150,151]. 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 utilized 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 [152].

Radon is also among the various hydrological, geochemical, geological and environmental species that have been employed in hydrological studies and for faults identification [153,154,155]. The shift in concentrations of $C{O}_2$ [156,157] near faults and the anomalous variations in groundwater levels [158], have been employed as well because the corresponding concentration variations, reflect the water-rock interactions [159] and the pathways generated by active faults [156]. Due to these properties, several hydro-geological species that have been utilised as tracers of pre-seismic activity [6,8]. Especially radon has been studied as well for co-seismic effects and tidal strain [160]. 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,161,162]. Under this perspective, radon is the best among the various hydro-geological species for earthquake forecast.

Radon combines hydrological, geological and environmental properties. Hydrologically, it dilutes to water [26] molecules and water aquifers. It is present in surface and, most importantly, underground waters [140,141]. Geologically, it is easily transferred within soil and rock reaching areas away [5]. Environmentally, it is naturally emitted and present in aspects of the environment, i.e., atmosphere, earth, water. It is naturally radioactive and easily detected. All these combination properties 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,161,162].

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,148,161,162,163,164,165,166,167,168,169,170,171,172,173,174,170,52,171,172,173,174]) and, recently, between radon TEC (please see Section 2.3.1) [132,133,134,175]. 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. [162], Conti et al. [8] and Huang et al. [161]. Several precursory signals have been obtained with passive techniques (no electricity needed), which offer rough time-series estimations, since these methods integrate the radon concentrations over extended periods of time (of at least >1-4 weeks), necessitated for the measurement. This roughness poses significant bias to the precursory estimations. Nowadays, radon precursory signals are monitored with active techniques (require electricity). The active techniques are implemented with portable monitors which allow for high rates of radon monitoring (typically between 1 $mi{n}^{-1}$ and 1 $hou{r}^{-1}$). As a consequence, these techniques offer detailed signals of radon and fine estimations (e.g., [5,6,8,17,161,162]). 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 concentrations levels [6,21,22,23,26,30,31,140,161]. 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, it appears that neither for mild earthquakes nor even for powerful earthquakes is there a universal model that can be used as a hallmark of a particular impending seismic event ([4,82,83,176], 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 emanations and this is done by various research groups [5,19,24,25,27,31,36,48,51,52,54,148,163,164,165,166,167,168,169,170,171,172,173,174,175,177,178,179,180,181,182,181,182,175]. The stability of the emission response of radon to seismic occurrences at the monitoring station, determines how successful these investigations are. Radon concentration in soil depends on a number of parameters and is, hence, 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,162,183]. Traditionally, because of well investigated relationships between radon and environmental parameters, the 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 though 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 [5,8].

Two approaches have been taken to the study of soil radon as an earthquake precursor: one involves doing simulation experiments in the lab and the other, involves monitoring the concentration of radon in soil gas at a specific place, over an extended period of time, in comparison to unusual emission changes in 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,162].

Based on observations and conclusions drawn from all of the aforementioned worldwide studies, it has been determined that some radon gas, which originates from the decay of radium in rocks inside the crust, 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 to an earthquake. According to Fleischer and Mogro-Campero [184], the deformation of rock mass under stress creates new channels that allow gasses from deep earth 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 [185]. Subsequently, a number of groups employed the concentration of radon in groundwater to study earthquakes [20,28,30,163,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201]. Groundwater radon concentrations frequently increase before earthquakes (e.g., [36,179,180,194,202]). A few times, before an earthquake, the amount of radon in groundwater’s has decreased [20,187,188,194]. Significant earthquakes may be related to groundwater radon’s peculiar behaviour as it offers information about subsurface dynamics [180], particularly in areas where high-stress buildup occurs in the crust [203]. The processes driving seismic activity can be better understood by examining the links between seismicity and geochemical signal variability [180].

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 [204] or escapes from it [205]. 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. [206], radon concentrations in 15 water samples from the Dehradun region, which is close to an active fault, range from 1.70 $Bq{L}^{-1}$ to 7.57 $Bq{L}^{-1}$. In southern Catalonia, 15 hot springs had groundwater radon levels ranging from 1.4 $Bq{L}^{-1}$ to 105 $Bq{L}^{-1}$ [207]. 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{L}^{-1}$ and 0.3 $Bq{L}^{-1}$, respectively [208]. Range values are provided globally [142,143,144,209]. Significantly more radon is found in groundwater in thermal spas [140,210,211,212].

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 [200,213]. 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 artificial neural networks and machine learning, have had some success recently [164,200,201,213].

3.2.3. Atmosphere

The primary source of atmospheric radon concentration is the exhalation from the earth and with a lesser extent the escape from surface and subsurface water [26]. Numerous processes are involved and meteorological elements have a significant impact on them [144]. Therefore, detecting anomalies in air radon in relation to 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 [167,214]. 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 [215].

Japan is the primary source of studies on earthquake forecasting using atmospheric radon. Iwatata et al. [167] reported that anomalies in the atmospheric radon concentration have been 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 [216], reported that two months before to the main shock of the 1995 Kobe earthquake ($M_w$=6.9; January 17, 1995, N${34.6}^{\circ }$, E${135.0}^{\circ }$), an increase in atmospheric radon concentration was noticed at Kobe Pharmaceutical University. Goto et al. [217] reported anomalous atmospheric radon concentrations associated with a shallow inland earthquake ($M_j$=5.5, depth=7 $km$; 5 July 2011, N ${34.0}^{\circ }$, E ${135.2}^{\circ }$) 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 of the residual values. The authors stated that the residual values stopped increasing on December 31, 1994, and, they concluded that this corresponded to the critical point of best fit model. These authors contended that rather than main stresses causing the Kobe earthquake directly, local stresses are responsible for the unusual ${222}^{Rn}$ fluctuation as well. Using the irreversible thermodynamic model, Kawada et al. [218] proposed that the preseismic radon shift was caused by a little change in crustal strain. Furthermore, a quantitative study by Omori et al. [215], revealed that the unusually high radon concentration (about 10 $Bq{m}^{-3}$) before the Kobe earthquake increased air conductivity and was sufficient to produce an ionospheric disturbances. Yasuoka et al. [219] claimed that further mechanically-induced precursors were seen prior to the Kobe earthquake. Igarashi et al. [190], for instance, described such precursory variations in groundwater radon concentration. Tsunogai and Wakita [220] documented further pre-seismic variations in crustal strain, groundwater discharge rate, and chloride ion content in groundwater. Because of the mechanical behaviour of the crust, these pre-seismic fluctuations should be related to one another [29,190,218]. 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. [190]. The linkage between preseismic fluctuations in the subsurface, atmosphere, and ionosphere could have been substantially verified if radon activity was clearly linked to the earthquake preparation process [190]. Additionally, current research supports the link between atmospheric radon and the Kobe earthquake [217,221].

4. Models

4.1. Electromagnetic Precursors Models

4.1.1. Models for the ULF Precursors

Three are the main models that have been proposed for the interpretation of the magnetic component of the pre-seismic ULF disturbances:

• Magneto-hydrodynamic model [222]. 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, the 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, the latter determining the relative significance of the convective and diffusive components.
• Piezomagnetic model [223]. This model states that an applied stress causes ferromagnetic rocks to shift in magnetisation, which in turn, induces a secondary magnetic field.
• Electrokinetic model [224]. 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. [225] established a theory about the current produced by charged distortions and currents induced by piezo-electric effects. The electrokinetic theory serves as the foundation for this theory. In water-saturated media with fluid-filled channels, electrokinetic currents can be found [226,227]. In order to model the parameters of these electrokinetic currents, Surkov et al. [228] assumed that an earthquake hypocentre is surrounded by water-saturated porous rocks with fluid-filled pore channels where cations from the fluid are 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 prepared, 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 [229] postulated that there is 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. [228] during the 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 inter-connectible small cracks. Furthermore, the rock’s conductivity together with the surface and bulk conductivities of the tiny fluid-filled cracks contribute to a non-zero conductivity of the surrounding space. However, according to Surkov et al. [228], 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. [228] limited the study by using a basic percolation hypothesis that ignored the crack-channel size distribution. The correlation length $\xi$ where $\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 three main ULF models have described successfully 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, at Romania [29,49,168,187,190,191,216,219,221,230].

4.1.2. Models for the HF Precursors

The behaviour of a stressed rock is comparable to that of an electromagnetically strained rock [176]. 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 [4,82,83,89,91]. The total of these emissions is 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 get more irritating until a multi-crack state occurs, indicating that local branching is the instability mechanism at work. It’s important to note that laboratory research has identified strong fracto-emissions during unstable crack propagation [22,23,34,83,231]. 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, 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 in-field during measurement and at laboratories under controlled conditions [4,34,82,89,90,176]. Owing to the previously indicated viewpoint, the main technique for forecasting earthquakes is to record the electromagnetic emissions from potential micro-fractures in the focal region prior to the final break-up [4].

As stated in by several papers (e.g., [82,88,91,176], 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 has been reported in order to obtain an understanding of the catastrophic nature of the fracture events. The investigation revealed that the system’s balance is progressively lost. This allowed for the estimation of the duration beyond which the process responsible for the pre-seismic electromagnetic emissions could continue as a non-equilibrium instability.

The analysis has indicated three key periods (i) the crucial epoch, that 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’s 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 has become blocked at the boundary of the strong asperities’backbone. Asperities are already under "siege" [83].

4.2. Radon Precursors Models

Scholz et al. [232] 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 zone’s 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 [177,233] states that the growing of tectonic stress forms a fractured focal rock zone. This zone gradually changes in volume and form over time. 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, according to the hypothesis of stress corrosion [234].

The LAIC model [37,95] describes also 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 micro-fractures 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 [235].

Nikolopoulos et al. [5,21,32,35,236], Petraki et al. [22,23], Alam et al. [36,179,180,202] and Petraki [86] proposed the asperity model [83] (please see Section 4.1.2) to explain radon emanation during preparation of earthquakes. Pre-seismic radon anomalies are attributed to variations of 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 begins critical persistent, strong anti-persistent and interchanged persistent, antipersistent radon anomalies occur. This has been associated with several earthquakes in Greece and China.

Other aspects have been expressed by other investigators. For example, Talwani et al. [237] reported that the anomalous behaviour of radon gas could be because of the opening of pore’s spaces during rock fracturing as a result of seismic events. Explosion tests have been 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 is a consequence of seismic waves applied to the rock. According to other investigators [238,239,240] 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,162]. Despite providing some indications, the visual observations are not enough to support the pre-seismic nature of the derived signals (e.g., [5,86], and references therein). Due to this, the analysis nowadays rely on the physical background of the related earthquake processes. The analysis’mainstream comprises the fractal methods [4,34,35,36,55,60,84,85,89,90,102,148,156,179,231,236,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257], methods based on the theory of information and entropy [82,88,91,176], symbolic dynamics [21,23,86,258,259,260,261,262] and natural time methods [64,263,264,265,266,267]. 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,256], Detrended Fluctuation Analysis (DFA) (e.g., [4,236]), Rescaled-Range Analysis (R/S) (e.g., [268]), Multifractal Detrended Fluctuation Analysis (MFDFA) (e.g., [36,55]), fractal dimensions from Katz’s, Sevcik’s and Higuchi’s methods [34,148], 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., [4,82]).

Spectral power-law analysis and Hurst exponent analysis have been utilised in all ULF, HF and radon precursors. DFA, fractal dimensions from the Katz’s, Sevcik’s and Higuchi’s methods and $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 been also 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,180,242,257,269,270,271] but will not presented here due to its complicated interpretation [272].

Due to their importance in both electromagnetic and radon precursors, the important properties of fractal behaviour, long-memory and Hurst exponent analysis are given in the following sub-sections, firstly. DFA is presented thereafter because it is a robust method that has been used in both LF, HF and radon precursors. The fractal dimension calculations through the 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 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 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 [273], which indicates whether the system is driven by linear mechanisms and order [274,275], is also connected with the scaling and fractal behaviour. The correlations are strong because a system’s complex behavior may be predicted by its fractal behaviour and vice versa.

5.1.2. Long-Memory

The long-memory [273,276,277] of a system can show if the system has long-range interactions or is random. In specific it may reveal if a geo-system has strong persistent and antipersistent behaviour ore if the long-range interactions are rather loose. If system exhibits long-memory, then the past, present and future states of the system are linked together in a manner that the presence of the system is not only derived from its past (Markovian behaviour) but defines also its future (non-Markovian behaviour) [82,176]. 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 [4,83,89,90]. 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 [278,279]. Time-evolving fractal events may be identified with the Hurst exponent and the corresponding time series’roughness can be evaluated [280]. Important details about the time series are revealed by the Hurst exponent’s value [242,278,279,281,282]:

(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 (anti-persistency);

(iii)

If $H=0.5$ the time series 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)$, will exhibit 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 [256,257]. Given its perceived benefits, the wavelet transform based on the Morlet base function are most frequently employed [5,34,82,84,176,231,236,246]. 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 forwards. Independent windows are also utilised as well, under the restriction that the square of the 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 get 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$, $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,....,N$) [34,86,283]:

(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, that 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 2 or 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{\frac{1}{N}{\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 of 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 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 uses the square of the Spearman’s ($r^2$) to assess the linear fit’s accuracy. According to Nikolopoulos et al. [34,231,236,283], good linear fits are considered as having $r^2\ge$ 0.95 or higher.

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

5.2.3. Fractal Dimension Analysis with Katz’s Method

The transpose array ${[s_1,s_2,...,s_N]}^{⊺}$ of the series $s_i$, $i=1,2,...,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$ [284,285]. 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 will stretch in the planar to d, if it does not cross itself, where d is :

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

(9)

By combining equations (7), (8) and (9), the Katz fractal dimension, D, becomes

$$D=\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),...,y\left(N\right)$$

(11)

recorded at $i=1,2...N$ intervals, the following is the construction of a new sequence, $y_m^k$ [247,248,286]:

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

(12)

The length of the curve associated to the time series is given by [286]:

$$L_m\left(k\right)=\frac{1}{k}\left(\sum _{i=1}^{\left[\frac{N-m}{k}\right]}y(m+ik)-y(m+(i-1)k)\right)\left(\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...k$, where $\left[...\right]$ is the Gauss notation, namely, the bigger integer part of the included value.

By inserting the normalisation factor

$$\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)

The Higuchi’ s fractal dimension, D, is finally calculated by the slope of the linear regression of logarithmic transformation of $〈L\left(k\right)〉$ versus k where $k=1,2,...,k_{max}$. It must be noted that the time intervals are $k=1,..,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...,$ for $k>4$ ($k_{max}>4$). Again $\left[...\right]$ is the Gauss notation [285].

5.2.5. Fractal Dimension Analysis with Sevcik’s Method

Using the approach of Sevcik [287], the fractal dimension of a time series is estimated from the Hausdorff dimension, $D_h$, as [285].

$${\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ϵ$ [285] 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 the fractal dimension of Sevcik’ with this transformation [285,287]:

$${\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 [278,279]. In accordance with the $R/S$ technique, the average, ${\langle x\rangle }_N=\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),...,x\left(N\right)$ into a new variable $y(n,N)$ in a specific time period $n,(n=1,2,...,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,278,279]:

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

(20)

The distance between the lowest and largest value of $y(n,N)$ a 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{\frac{1}{N}{\sum }_{n=1}^N{(x\left(n\right)-{\langle x\rangle }_N)}^2}$$

(22)

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

$$\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\left(\frac{R\left(n\right)}{S\left(n\right)}\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. Following some of these empirical relations are given.

Rikitake [288] 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 [289] suggested and empirical earthquake forecast model as

$$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 [290] provided another model associating the precursory time T in days and the magnitude M of an earthquake as

$$logT=A+B·M$$

(27)

where A and B are statistically determined coefficients.

Dobrovolsky et al. [291] 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 [292] suggested that the epicentral distance D in $km$ and the magnitude M of an earthquake are associated as

$$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 [293] 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 [294] proposed a different relation that combined the epicentral distance D in $km$ and the magnitude M of an earthquake as

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

(31)

whereA=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 as

$$R_E=D·R$$

(32)

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

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

(33)

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

Chetia et al. [165] 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\left(km\right)$ and magnitude M ($M_w$) are 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 [4]. Moreover, the earthquake generation processes are multi-facet [4] and therefore a combination of techniques is needed [4,34,36,82,272] 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, make use of these estimations. For several scientist these are considered as adequate and enough.

7. Table of Papers

Table 1 presents a collection of papers for electromagnetic precursors. Table 2 shows the papers collection for radon precursors. The papers in both tables are presented chronologically and therefore old events are also included. To avoid unessential records, historical earthquakes are limited to very strong and extremely strong. Although the knowledge and methodologies have evolved, the techniques of treatment of these old earthquakes are not definitely new, since they refer to the available ones of the publication time. The historical electromagnetic precursors also include the great papers that pioneered and opened breakthroughs in the seismic knowledge. Especially in radon precursors, the data include traditional anomaly treatment that hold up-to-date, that is the grade of anomalous behaviour in respect to the baseline and the anomaly duration. 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. For the reason that the papers are many, a special presentation approach was selected so as every row to present all the data of each earthquake and the maximum of available information. Since every paper is special and the published information is not uniform, a certain variety of data that are available in most of the papers have been included as column names, but effort has been put to present special information also. All papers have been accessed in the site of each journal and the available information were downloaded as 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.

Evaluating Table 1 and Table 2, it can be supported that the majority of publications are based on visual observations of collected data and some extend of statistical analysis. This can be explained by the fact that both for historical and new earthquakes it is very difficult to collect data from at least one station nearby. 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, generally mentioning, effective epicentre’s area. As mentioned in several publications and expressed collectively in Eftaxias [176], there is no one to one correspondence between earthquake occurrence and anomaly detection. Moreover, even the most advanced methods do not manage to have 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 preparation of earthquakes. Very robust method is the natural time analysis with has much more to give. The satellites are now several and can be accessed conveniently. This provides new insights in the related research. Remote sensing and SAR techniques are very powerful as well. Not omitting of course the number of installed stations worldwide. These new tools give boost to the modern approaches which are also multi-facet and with the collaboration of different groups.

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

The collaboration between scientists and the multi level approaches with different methodologies is the key point to the seismic reasearch in the following years. This research is ongoing and in a continuous search for credible and powerful precursors.

8. Conclusions

This paper is a review on electromagnetic and radon precursors for earthquake forecasting. The electromagnetic precursors emerge in a diverse frequency band ranging from ultra low to very high frequencies. Nowadays electromagnetic data are collected from satellites, whereas the remote sensing techniques are in increasing usage also. Within the electromagnetic spectrum is also the TEC measurements and the modern approach of SAR studies. The various investigators are still working independently, nevertheless there is great space for inter-collaborations. The traditional approach is still the recordings from ground stations, with the precursors of the ULF range to have the greater history and potential. $MHz$,$KHz$ frequencies provide also very good estimations. 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 assist in forecasting of earthquakes from far.

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 the recent years advanced methods have been published and used in several new publications. Special mention is in 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. Usage has been done to $R/S$ analysis, fractal dimension analysis and Hurst exponents. Block entropy and several entropy measures have been used as well. 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 preparation of earthquakes. Much research needs to be implemented and new approaches are still on 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 e 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 characterised so, after the occurrence of the earthquakes. This is a disadvantage that will be overcome as the research is progressing. There are papers that forecast earthquakes prior to their occurrence and this is still the most distinguishing issue. The work of researchers from different sub-disciplines of electromagnetic and radon precursors will provide better forecasting results in terms of science.