A Critical Review of Ground Based Observations of Earthquake Precursors

We aim at giving a short review of the seismo-associated phenomena detected on ground that in recent years have been investigated as possible earthquake precursors. The paper comes together with a companion article–published on this same volume by Picozza et al., 2021–devoted to summarize the space-based observation of earthquake–precursors by satellites missions. In the present work, we give an overview of the observations carried out on ground in order to identify earthquake precursors by distinguishing them from the large background constituted by both natural non-seismic and artificial sources. We start discussing the measurements of mechanical parameters and variations of geochemical fluids detected before earthquakes; then we review thermal and atmospheric oscillations; finally, observations of electromagnetic and ionospheric parameters possibly related to the occurrence of impeding earthquakes are discussed. In order to introduce a so large field of research, we focus only on some main case studies and statistical analyses together with the main hypotheses and models proposed in literature in order to explain the observed phenomenology.


INTRODUCTION
Earthquake prediction is the Saint Graal of seismology, but the study of possible earthquake precursors should be better regarded in the framework of fundamental geophysics more than in trying to guess the future Hough, 2020). Several authors have analyzed data and many papers have been published with claims ranging from (by chance) observations of spatial and temporal correlations -cautiously interpreted as possible earthquake precursors -up to the proposal of methods and procedures (never confirmed) aimed at forecasting earthquakes. The realm of studied earthquake precursors includes a variety of physical parameters ranging from mechanical deformation up to gas emissions; from variations of groundwater levels up to fluctuations of electromagnetic field (in a large spectrum of frequencies, possibly radiated, induced or generated as secondary effect of other primary perturbations); from variations of ground temperature up to fluctuations of ionospheric and magnetospheric parameters. The first reports about correlation between impending earthquakes and electromagnetic emissions date back to Varotsos (1981) (even though highly debated), Gokhberg et al. (1982) and Warwick et al. (1982). About 50 years ago, the successful short-term prevision of the strong Haicheng (China) earthquake (Lomnitz, 1994) as well as the failure in forecasting the event of Tangshan (China) (Lomnitz, 1994)-even though of comparable intensity -have started a swinging wave of hopes and delusions about the possibility that seismic events could be effectively forecasted (Turcotte, 1991;Uyeda et al., 2009b). Many measurements, claimed as earthquake precursors, have been carried out (occasionally or even by systematic observation campaigns) on single earthquakes (case studies). Unfortunately, only a few of them were repeated/reproduced in occasion of other seismic events (such as for example Hattori et al., 2013;Liu et al., 2010;Davidenko and Pulinets, 2019) so that the observations of these cases studies are in some aspect often "unique" and scattered in the panorama of earthquake precursor investigations. For example, some early observations have suggested that, before strong earthquakes, the focal area could generate and radiate detectable electromagnetic signals in a large range of frequency. Two hours before the M 9.2 Alaska earthquake of 1964 -one of the largest ever recorded seismic event in the era of regular seismological recording- Davies and Baker (1965) reported a strong ionospheric anomaly at about 4-5 MHz, recorded at Boulder, Colorado. Another "classical" case study was the measurement of fluctuations in magnetic field data (in ultra-low frequency range) carried out a couple of weeks (and then some hours) before the Loma Prieta earthquake (magnitude 7.1) on 1989 . A less numerous type of analysis is constituted instead by statistical studies, performed by applying a given data analysis procedure to an ensemble (more or less statistically significant) of seismic events. In this framework, a systematic effort have been carried out by USGS along the San Andreas Fault at Parkfield, CA, by installing creep and strain meters, groundwater level detectors, magnetometers, etc. (Bakun and Lindh, 1985), but gathered data and results have not allowed to forecast any seismic event (Langbein et al., 2005). The most decidedly adverse perspective, mainly from the seismology and geodesy point of view, was summarized by Geller according to whom:<< results in nonlinear dynamics are consistent with the idea that earthquakes are inherently (or actually) unpredictable because of the highly sensitive nonlinear dependence on initial conditions>> (Geller, 1997), while indirectly Uyeda replied that: <<There are reasons for this pessimism because mere conventional seismological approach is not efficient for this aim. Overturning this situation is possible only through multidisciplinary science>> (Uyeda et al., 2009b). This paper aims at presenting a critical overview of earthquake precursors observed on ground and the related analyses published in literature trying to point out their own characteristics (physical parameters, intensity, duration, background, etc.) and possible connection to the earthquake magnitude. We will also summarize some of the physical models proposed to reconcile the observed phenomenology with the physics of earthquake even though a scientific consensus of which could be preferable is still missing. We address the reader to the paper (Picozza et al., 2021) published in this same issue that provides a review of satellite-based observations and more in general of earthquake precursors in space.

SUMMARY OF CLAIMED GROUND-BASED EARTHQUAKE PRECURSORS
Even though our review does not pretend to be exhaustive of earthquake precursors on ground, we tried to select them between the most reliable ones published in literature. In the following, we will discuss several types of precursors that suggest an effective lithosphere-atmosphere-ionosphere-magnetosphere coupling, such as: 1) Seismicity that is the most extensively studied phenomenology before, during and after earthquakes (Mignan 2008;Hong et al., 2018;De Santis et al., 2019b). Also extreme low frequency acoustic emissions have been claimed (see for example Ihmlé and Jordan, 1994) to constitute earthquake precursors before the main rupture at a higher frequency. 2) Lithospheric mechanical deformations, such as those detectable with creep-and stain-meters (Niu et al., 2008;Langbein, 2015). 3) Variation of the groundwater level and composition, reported some weeks up to few hours before earthquakes (Hayakawa et al., 1997;Koizumi et al., 1999). 4) Gas exhalations, mainly (but not only) of radon or radioactive ions induced by gas-water release from earthquake preparation zone into the atmosphere (Khilyuk et al., 2000;Pulinets et al., 2003). 5) Fluctuations of temperature observed in temporal correlation with some earthquakes and possibly reconciled with variation of groundwater circulation and uplift or more recently with vapour condensation on surface (Tramutoli et al., 2005). 6) Propagation of acoustic gravity waves (AGW) (Molchanov and Hayakawa, 2008), A physical mechanism of seismoionospheric coupling including both AGW and radon exhalation has been recently suggested (Rapoport et al., 2020). 7) Fluctuation of electric and magnetic field components in a large range of frequencies [from ULF (Uyeda, et al., 2009b;Han et al. (2014)] to VHF (Sorokin et al., 2020). Many observations have been reported on ground and in space of (direct, induced and secondary) electromagnetic emissions localized on the earthquake area or measured along the related field line or spread around it. 8) Ground based observations of ionospheric parameters [such as Total Electron Content (TEC) (Liu et al., 2004;Liu, 2009), VLF reflection height (Hayakawa et al., 1996(Hayakawa et al., , 1997Rodger et al., 1999), whistler dispersion , critical frequency foF2 (Hobara and Parrot, 2005), etc.].
With the analysis of LEO satellite observations, the range of earthquake precursors investigations has been extended including measurements of disturbances of plasma parameters in the ionosphere-magnetosphere transition zone; thermal anomalies; AGW; precipitation of particles from the inner Van Allen belt (see Picozza et al., 2021 and references therein).
For each type of precursors listed above, we will try to discuss the main related information such as earthquake parameters (including time, location, magnitude and depth), time of the measured precursor, duration of the disturbance, amplitude, signal-to-noise ratio, distance from the epicentre, etc. Moreover, we have catalogued the physical models proposed in literature in order to verify which of them could better explain the phenomenology of the precursors and can reconcile the

Seismicity and Observation of Local Mechanical Deformations
In the seismic sequence-often associated to large earthquakes-the main shock is frequently preceded by foreshocks of lower magnitude (Jones and Molnar, 1976;Reasenberg, 1999). Foreshocks occur within a variable time interval from the main earthquake (even though often with an increased frequency before the main shock) and spatially near to the epicentre of the main event (Sornette and Sornette, 1989;Ben-Zion et al., 2003). Scholz (2019) suggested that foreshocks are part of the nucleation process resulting in the main seismic event and that, in some sequence, dilatancy could induce seismic energy release and explain short-term quiescence just prior to the main event (Press and Siever, 1982;Allegri et al., 1983;Lomnitz, 1994;Sgrigna and Malvezzi, 2003). A systematic worldwide study of foreshock catalogues for low magnitude earthquakes is a challenge due to the uneven threshold in seismic measurements around the world. It has been suggested that the foreshock occurrence and the claimed earthquake precursors (such as radon release, electromagnetic anomalies, groundwater level variations, etc.) could be correlated through some physical mechanism (De Santis 2014; Varotsos et al., 2019). After a strong earthquake, the aftershock sequence normally decays and the occurrence probability of a subsequent larger event is of a few per cent Gulia, and Wiemer (2019). Anyway, this probability is a function of the stress conditions due to previous earthquakes and long-term tectonic conditions. In this framework, foreshocks can give valuable information about the process in action (see Console et al., 1993;Avlonitis and Papadopoulos 2014; and references therein), but earthquake catalogues are still uncomplete and foreshock interpretation for earthquake forecasting is highly debated. Attempts to predict the next large earthquake, based on physical models and Coulomb stress transfer, have been unsuccessful due in part to incomplete knowledge of the location of faults Nanjo (2020). Gulia and Wiemer (2019) suggested that in some cases it would be possible to distinguish between decaying aftershock sequences and foreshocks preceding a large event. Gulia and Wiemer (2019) proposed that the probability of larger subsequent event is higher for seismic sequences diverting from the generally observed increase of b value after a mainshock. Anyway, the authors consider preliminary their results due to the reduced number of events analyzed (M > 6 with high dense seismic networks coverage). More recently, Trugman and Ross (2019) suggested that more than 70% of all mainshocks in Southern California was preceded by foreshocks, but van den Ende and Ampuero (2020) objected that only a percentage between 18 and 33% of mainshocks were preceded by significantly elevated seismicity rates. Several authors have reported surface deformations, such as tilt, strain, uplifts and downdrops, etc., measured before earthquakes of medium and high magnitude (Rikitake, 1987;Lomnitz, 1994). Some mechanical models Tse and Rice (1986) and Lorenzetti and Tullis (1989), etc. have been proposed in order to predict surface deformation possibly associated to the earthquake preparation phase. They need to describe the fault mechanical dynamics through constitutive relationships and to study the friction along the fault (Dieterich 1978;Dieterich, 1979;Ruina, 1983) in order to estimate if the induced strain rate, displacement, velocity, etc. are below the detectability threshold for the available instruments and experimental methodology. The results suggest that-due to signal-to-noise ratio, with respect to the intrinsic background-the strain rate is the most reliable observable. These investigations today can take advantage of temporal high-frequency sampling and spatial high-resolution measurements with GPS and SAR interferometry with satellites, in particular at low orbit. In the case of Amatrice earthquake, Panza et al. (2018) reported an increase of the deformation velocity for a transect (500 km wide) moving eastward, along the direction of tectonic extension (from Tyrrhenian to Adriatic coast). The velocity gradient has a peak localized in the Amatrice area (not observed for transects to the North and South). For further case studies see for example (Wright, 2016;Moro et al., 2017;Panza et al., 2018) and references therein. Even though the intensity of the measured deviations-carried out within some hundreds of kilometres from the epicentres-seems related with the earthquake magnitude, the claimed identification of these observations as earthquake precursors is not conclusive.

Variation of Geochemical Fluids: Release of Gas and Groundwater Level Fluctuations
Starting from the sixties, several reports have been published about a claimed increase of radon concentration before earthquakes (Lomnitz, 1994). The hypothesis is that an increase of compressional stresses could open and/or close micro-fractures and cracks facilitating radon exhalation of radon up to the surface together with the flow of groundwater carbon dioxide, methane, helium, etc. (Wakita et al., 1980;Teng et al., 1981;Biagi et al., 2001b). The 222 Rn is the most stable radon isotope, a noble radioactive gas generated by the alpha decay of Radium 226, with a half-life of about 3.8 days (Bé et al., 2011). 222 Rn is water soluble, with low concentration in surface water-due to the continuous release in the atmosphere-and higher concentration in deep groundwater. For this reason, the variation of concentration of 222 Rn has been studied as a marker of tectonic processes and proposed as a possible short-term precursor of seismic events (see for example Richon et al., 2003;Omori et al., 2007;Alam et al., 2020).

Physical Mechanisms for Correlating Gas Exhalation and Seismic Activity
The variety of models suggested in order reconciling the observed radon exhalations with physical mechanisms driven by seismicity can be summarized into 1) releases induced by ultrasonic vibrations; 2) pressure-induced change of solubility; 3) collapses of pores; 4) development of new pores; and 5) fluid mixing in depth (Thomas, 1988).
1) According to the vibrational model, gas releases would be induced, or facilitated, by ultrasonic vibrations of rocks. Even though this mechanism has been tested in laboratory and during artificial explosions in seismic explorations, the power density of natural seismic spectra at high frequencies could not be enough for explaining the phenomenology of the seismo-associated exhalations. Moreover, releases induced by seismicity are more intense than the explosion induced ones, and these, besides, can follow (instead of precede) the rupture events. 2) It has been suggested that gases emissions could be induced by variations of solubility due to an increase of the fluids pressure during the earthquake preparation phase, but the mechanism would not be effective because the required increase seems too high to be transferred to the fluid phase (Cicerone et al., 2009).
3) The collapse of the pores volume due to the stress increase of the incoming earthquake region has been also claimed to induce gases release in groundwater. Even though observed in some laboratory tests, this mechanism is questionable because, generally, high stress values on porous rocks are effective in increasing the pores volume and the observed periodicity in the gas releases intensity is not easily reconcilable with the irreversible compression of pores. 4) A more effective role in gas exhalation could be played by the rocks dilatancy that can increase of tens and hundreds the rock porosity percentage (Brace, 1978;Bernabé et al., 2003;LongJohn et al., 2018): microfracturing both facilitates gas release from the rocks and increases the reaction ratio with ground waters through the growth of microscopic surfaces (Holub and Brady, 1981). On the other side, significant pores volume increases have been observed only near the rock failure strength, which would mean that the claimed mechanism could be effective only in the small volume of the seismic fault experiencing the rupture process, whereas gases releases have been observed even far from the epicentre (Pulinets and Ouzounov 2011). 5) The phenomenology of spread exhalations could be reconciled by invoking cracking due to corrosion and subs stage occurring at low stress with higher fluid content. Finally, it has been proposed that a mix of ground fluids and chemical elements between different aquifer systems could be the most effective mechanism for generating fluctuations (both positive and negative) of gas exhalation even because it explains also the temperature fluctuations some time measured together with radon exhalations (Byerlee, 1993). In this framework, the mechanism would have a role in the electrokinetic generation of low frequency electromagnetic emissions (Fenoglio et al., 1995) that will be discussed later.
All the physical models claimed to explain precursors share the hypothesis that fast and non-linear processes in the rocks along the seismic fault (such as deformation, dilatancy, fluids flow changes, pores volume variation, etc.) could originate the anomalous variations of observed parameters (Press and Siever, 1982;Lomnitz, 1994).

Radon Exhalations
The majority of the reports about significant variations of radon concentration was for moderate and strong magnitude seismic events (about M ≥ 4.0), but fluctuations have been reported also for earthquakes of lower intensity. More than 80% of the measurements are constituted by increases of radon concentration with respect to the background reference value, with a distribution of the variations peaked between 50 and 100% and extending up to more than a thousand percent. The duration, as well as the beginning and ending time of the radon exaltation events, show a quite large variability-without a clear temporal distance with the claimed associated earthquake-that does not allow identifying time of radon exaltations as a sharp/reliable short-term earthquake precursor (Plastino and Bella 2001;Cicerone et al., 2009;Plastino et al., 2010;Sorokin et al., 2020). In fact, even though, in the majority of the cases, the anomalous fluctuations started within about one month before the earthquake and lasted less than about three months and a half, the statistics also includes radon exhalations ended before the earthquakes and/or continued after the seismic events (Pulinets and Ouzounov 2011). On the other side, Yasuoka et al. (2006) demonstrated that radon behaviour fits curve of the critical exponent (Sornette and Sammis, 1995).
The reported radon exhalations are more frequent near the epicentres, where also the highest variations have been measured. On the contrary, no significant correlation has been observed between the intensity and frequency of the gas exhalation from one side and the earthquake magnitude from the other one (Kissin and Grinevsky, 1990;Toutain and Baubron, 1999;Pulinets and Boyarchuk, 2005). That would suggest that, even though gas releases could be related to the seismic fault dynamics, their entity would not be correlated with the incoming earthquake intensity (Cicerone et al., 2009;Sorokin et al., 2020).
Both the advance in time of the radon exhalation (before the seismic event) and its duration are correlated with the earthquake magnitude, suggesting that bigger events are preceded by larger releases occurring more in advance (see for example the statistical study of Kissin and Grinevsky, 1990;Toutain and Baubron, 1999; as well as the most recent works Pulinets and Boyarchuk, 2005;Hartmann and Levy, 2005;Riggio and Santulin, 2015). This would furtherly support the Rikitake law of the empirical linearity between the logarithm of precursor time and the earthquake magnitude [see Figure 13 of Rikitake (1987)]. On the other side Pulinets and Ouzounov 2011 pointed out that the discrepancy between the claim that radon exhalations are a precursor to earthquakes (Toutain and Baubron, 1999;Omori et al., 2007;Pulinets, 2007) and the demonstration that radon releases are not a statistically reliable precursor (Geller, 1997) cannot be easily resolved with local measurement stations, due to the variety of measurement methods adopted worldwide and the uncertainty about the origin (crustal or mantle) of radon and the related transport models. Due to the difficulty of monitoring radon on the ground-with networks with high spatial resolution (İnan et al., 2008) over large areas (such as within and beyond the Dobrovolsky radius)-the authors suggested retrieving radon release (as well as methane, carbon dioxide and other geochemical fluids) indirectly, through remote sensing methods, from space-based observations of thermal infrared radiation, based on the LAIC model that correlates thermal anomalies with radon and other gas exhalations Surkov, 2015;De Santis, et al., 2019b). However, the claim that closer to and on the eve of an earthquake a higher release should be observed is not yet confirmed and a statistical assessment of the true/false vs. positive/negative cases, through a reliable confusion-matrix based classification, is desirable. The recent statistical analysis of correlation between groundwater radon variations and seismo-tectonic activity of time series about the Wenchuan earthquake (Alam A. et al., 2020) would show a persistent trend with a notable upsurge just before and during the Wenchuan earthquake in the near stations not observed in the response of more distant monitoring stations.
Based on the literature, at present, a quantitative estimation of the intensity of the correlation mechanism, as well as of the temporal and spatial distance between claimed anomalous gas releases and the occurrence of seismic events are still missing. Therefore, no firm conclusions can be drawn about the reliability of the hypothesis.

Aerosols and Bubbles Exhalations
Aerosol measurements showed anomalous variations before and after the Gujarat M 7.7 event on January 26th, 2001 (Okada et al., 2004)-that was also preceded by the anomalous disturbances of TEC, ions and electron density, electron temperature and VLF electric fieldand before the Chile M 8.8 earthquake on February 27th, 2010 (Akhoondzadeh, 2015). By analysing data of Aerosol Optical Depth (AOD) (measured through absorption of light at 550 nm by MODIS of Terra and Aqua satellites) Qin et al. (2014) have shown a significant fluctuation along the active Longmenshan faults, seven days before the Wenchuan earthquake (May 12th, 2008) before whom several have been the detections of anomalous variations of electromagnetic, atmospheric and ionospheric parameters including air temperature, outgoing longwave radiation, relative humidity, etc. Similarly, Akhoondzadeh and Jahani Chehrebargh (2016) reported anomalies in the AOD detected by MODIS before 16 high magnitude earthquakes; Boyarchuk (1997) shown exhalation of metallic aerosols (Cu, Fe, Ni, Zn, Pb, Co, Cr and radon); King (1986) reported emissions of Rn, He, H, Co 2 up to 500 km far from the epicentre (from few weeks to hours) before the seismic event; several analyses (Alekseev and Alekseeva, 1992;Heinicke et al., 1995;Pulinets et al., 1997;Biagi, 2009) claimed an increase of up two orders of magnitude of the charged aerosols density and an enhancement of the local radioactivity from weeks to days before earthquakes (correlated with possible exhalations of radon and other radioactive species).
Increased gas release in seismic regions during the earthquake preparation phase has been observed not only on land near faults (King, 1986), but also near submarine faults (McCartney and Bary, 1965;Lyon, 1974) with an intense rise of gas bubbles (including vapour, CO 2 , He, methane, etc.) due to volcanic activity (Marty et al., 1993;Nikolaeva et al., 2009). The bubbles of gas under the sea can carry small electric charge (10 -14 -10 -13 ) C (Gak 2013) that can originate electric field in the atmosphere (Harper 1957;Blanchard, 1963). Even though the detailed balance of air-sea fluxes is still not completely known (Zavarsky et al., 2018), marine spray aerosol strongly contributes to an increase in aerosol optical depth (Revell et al., 2019) and winds spread aerosols into the atmosphere both over the land and over the sea. Based on that Sorokin et al. (2020) advanced the hypothesis that the electromagnetic environment on the lands and over the sea could be similarly affected by aerosol releases, as would be confirmed by the analyses of electric field perturbations detected by satellite, that don't show a significant difference for seismic events on land and at sea.

Groundwater Level and Vapour Variations
Studies about level variation of ground waters in occurrence with seismic events have been carried since long time ago and are particularly numerous (Hamilton, 1975;Kovach et al., 1975;Raleigh et al., 1977;Cai and Shi, 1980,;Merifield and Lamar 1981;Golenetskii et al., 1982;Wakita et al., 1985;Asteriadis and Livieratos, 1989;Liu et al., 2006). Between the end of last century and the begin of the current, the investigations extended including test campaigns in other countries (Kissin and Grinevsky, 1990;Igarashi et al., 1992;Igarashi et al., 1995;Roeloffs and Quilty, 1997;Koizumi et al., 1999;King, 1986;Chadha et al., 2003;Koizumi et al., 2004). In many cases, the reports about level variation of ground waters in occurrence with seismic events are not systematic and "scattered", making difficult comparisons and a statistical analysis. Nevertheless, the trend of the main characteristics (such as distribution, spatial and temporal distance, intensity, etc.) seems qualitatively similar in the gas releases and ground water level measurements suggesting Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 676766 5 a common driver in the earthquake preparation phase that, anyway, asks for further confirmations. In order to point out the correlation between variation of ground water level and seismic events it is needed to reject the background effect (Roeloffs, 1988) due to not seismic processes such as tidal effects (Bredehoeft, 1967); rainfalls, pressure variations and seasonal contributions (Rice and Cleary, 1976;Roeloffs and Quilty, 1997), oil and gas extractions, etc. Even though measurements of water level fluctuations are spread up to values of some meters, the most frequently absolute variations are within 1 m with a large majority of decreases before the earthquake. The statistical distribution of advance time between ground water level variation and seismic events ranges between less than one day up to years with a peak at about one month and half, while the largest anomalies would occur nearer to the earthquake time. In addition, the largest variations have been observed nearer the epicentre (with the majority within 200 km), but there is not a clear evidence of a correlation between the entity of the variation and the earthquake magnitude. More recent observations (see for example Plastino, 2006;De Luca et al., 2018) have confirmed the general picture of the impact of seismic activity on the groundwater level fluctuation.İnan et al. (2010) reported hydro-geochemical anomalies lasting for more than a month before an earthquake magnitude 4.8, in Western Turkey, detected within few tens of kilometers from the epicenter. A longer preparation phase seems confirmed also by the analysis of De Santis et al. (2020) of a precursory anomaly lasting almost a year (from September 2018 to July 2019) of ground water level data from a borehole located about 200 km from the epicenter of the July 6, 2019, M7.1 Ridgecrest earthquake and six USGS groundwater-monitoring sites.

Thermal Anomalies Possibly Associated With Seismic Events
In the framework of earthquake precursor's research, several analyses have investigated the existence of anomalous temperature variations before earthquakes (Tronin 1996;Carreno et al., 2001;Tronin et al., 2002b;Jing et al., 2013;Ouzounov and Freund 2004;Saraf and Choudhury 2005;Choudhury et al., 2006;Ouzounov et al., 2007;Panda et al., 2007;Bi et al., 2009;Saraf et al., 2009;Qin et al., 2013;Lu et al., 2016). Relatively few are the reports about such anomalies detected on ground before earthquakes, the reason being the reduced number of test campaigns devoted to such a measurements and the complexity of the study aimed at discriminating positive effects from those not related to earthquakes. Now day, the satellite remote sensing has extended the available dataset as well as the capability of detection.

How Earthquakes Could Generate Thermal Anomalies
The first hypothesis advanced in order to reconcile the possible appearance of a thermal anomaly in conjunction with earthquakes has been that such an anomaly would originate in depth during the seismic preparation phase (for example by frictional heating on fault surfaces) and would propagate up to the surface. Due to the low thermal conductivity of the rocks, the direct diffusion of an anomalous thermal fluctuation is quite inefficient in reaching the surface asking for some years even for diffusing of few meters. An alternative idea was that thermal fluctuations could be generated by fluids flow variation in depth-such as uplift or interruption of the hot geothermal groundwater circulation in volcanic regions and/or due to the volumetric variation of porous rocks due to dilatancy that facilitates fluids and gases flow-and then transported to the surface by the groundwater motion. The speed of the flow, the volume of groundwater involved as well as the depth variation would drive the temperature fluctuation intensity and temporal trend. Being the geothermal gradient of about 1.5-3.5°C per 100 m (but even higher for volcanic regions), a temperature variation of several degrees can be generated by the circulation trough rocks at some hundreds of meters. The hypothesis of a fluid driven propagation would explain why the majority of the claimed observations observed on ground has been carried out near thermal sources (in Greece and Japan, where significant subductions and geochemical activity take place) (Zheng et al., 2020); within few tens of kilometres from the epicentre; with a peak of the intensity distribution at less than one degree; before and during the seismic events. Anyway, it must be noticed that no such geothermal anomalous fluctuation have been reported along the Saint Andreas fault (Lachenbruch and Sass, 1992). Moreover, because thermal anomalies have observed both on ground and over the sea surface, Pulinets and Ouzounov (2011) ruled out thermal waters as a source to explain observed anomalies.
A further hypothesis has been advanced for the generation mechanism according to which the variation of tidal force could trigger shallow thrust fault earthquakes (Cochran et al., 2004). In fact, Tanaka et al. (2002) have found a correlation between temporal trend of tidal force and seismic events. It appears reasonable that the tectonic stress due to tidal force can contribute to reach the rock critical failure point (Heaton, 1975;McNutt and Beavan, 1981;Kilston and Knopoff, 1983) and then could be correlated with the earthquake nucleation (Weiyu et al., 2018). In this framework, it has been claimed that tidal force could generate some seismo-associated thermal anomalies Ma et al., 2015;Yan Z. et al., 2017;Weiyu et al., 2018;Zhang et al., 2018). However, even this generation mechanism does not overcome the limits due to the inefficient propagation process up to the surface.
An interesting link could exist between thermal anomalies and radon exhalation. Even though Toutain and Baubron (1999) (see also references therein) have shown some systematic correlation between radon exaltation and seismic events, further investigations (for exampleİnan et al., 2008) have not confirmed previous results and other scientists (such as for example Geller, 1997) have demonstrated that the claim of using radon as earthquake precursor could be illusory. Pulinets and Ouzounov (2011) and Sorokin et al. (2020) proposed that the discrepancies between data, models and statistical accuracy could be due to the variety of instruments and methodology adopted for the measurements campaigns as well as to the mosaic-like distribution of radon exhalations that could appear randomly scattered in space and time. By assuming that radon exhalations would result into local increase of temperature, Pulinets and Ouzounov (2011) claim that the remote sensing of thermal anomalies (propagated from ground to upper atmosphere) would solve the difficulties of point like observations of radon on ground (Pulinets, 2007).

Detections of Thermal Anomalies
The prospective of investigating seismo-associated thermal anomalies has strongly changed with space-based observations. The Earth's surface naturally emits thermal radiation that is routinely measured by satellite sensors in the (8-14 μm) spectral range (Pergola et al., 2004;Genzano et al., 2007). Several different methods have been applied to identify thermal fluctuations. Tramutoli et al. (2001) suggested the RST (Robust Satellite Techniques) procedure; Ouzounov and Freund (2004) analyzed LST (land surface temperature) data; Bryant and Nathan (2003) suggested the NTG (Nighttime Thermal Gradient) method and finally Zhang et al. (2010) applied wavelet transformation to infrared data. On the other side, various thermal parameters and datasets, collected on ground and trough remote sensing have been studied (Heaton, 1975;Tanaka et al., 2002;Blackett et al., 2011). Several studies have analyzed TIR (Thermal InfraRed spectral radiations) measurements acquired by infrared sensors on board of several satellite missions such as MODIS (on Terra and Aqua), FY Chinese satellite, MSG-41 SEVIRI on MeteoSat missions, AVHRR on NOAA satellites, etc. (Weiyu, et al., 2018;Yan R. et al., 2017;Zhang and Meng, 2019;Zhang et al., 2021). A common feature of the cited processing methods is to statistically process long temporal series of remote sensing data in order to extract possible anomalies (Lisi, et al., 2010;Sannazzaro, et al., 2014). Even though this approach provides a better statistical robustness, on the other side it could prevent from detecting small but significant thermal anomalies (Weiyu et al., 2018). One of the last research paper on the geochemical and thermal seismo-associated anomalies by Martinelli et al. (2020a) provides also an interesting summary and further references on this subject.

Acoustic Gravity Waves
At the begin of 2000 several authors (Gokhberg and Shalimov, 2000;Molchanov et al., 2004;Korepanov, et al., 2009) suggested that the earthquake precursor's phenomenology could be explained by the emission of AGW oscillations in a range of 6-60 min generated by lithospheric oscillations or gas exhalations and amplified in their propagation from the Earth surface through the atmosphere up to the ionosphere. The model of Molchanov et al. (2004) is based on the emission of charged aerosol (including an increase of groundwater, stable/radioactive gases and heat flows) in the earthquake preparation phase that would generate local disturbances in the conductivity in the atmosphere-ionosphere circuit. The advantage of this approach would be its capability to explain the large variety of observed disturbance taking in to account anomalous electromagnetic ULF/ELF measurements on ground (Schekotov et al., 2006) and in space (Rozhnoi et al., 2004), small scale plasma irregularities (Molchanov, 2009) in the geomagnetic flux tubes possibly reconciled with earthquakes occurrence, thermal fluctuations  as well as the claimed preearthquake character of VHF measurements (Devi et al., 2012). The analyses of ultra-low frequency electric and magnetic field before earthquakes shown electromagnetic perturbations associated with growth of plasma density structures (with scales between 10 and 40 km) several days before earthquake in the magnetic tubes with the footprint over the epicentre (Zhang et al., 2009;Akhoondzadeh 2013) similarly to what occurs over tropical cyclones (Isaev et al., 2002;Sorokin et al., 2005). It has been suggested (Akhoondzadeh, 2013) that these seismo-associated plasma fluctuations would be generated by quasi-static electric field amplification due to AGWs instability for value of critical field greater than about 10 mV/m (Chmyrev et al., 1997;Sorokin et al., 1998Sorokin et al., , 2005. The plasma irregularities along the magnetic tube show variations of tens of percent and dimensions from few hundreds of meters to tens kilometres (Chmyrev et al., 2008) while amplitude of magnetic field fluctuation are of few nT and with frequencies from fraction of Hz up to few Hz (Chakraborty et al., 2018 and references therein). Beyond several advantages (Lizunov, et al., 2020), the models based on AGW and internal gravity waves suffer of the difficulty in explaining local disturbances (such as electromagnetic and plasma fluctuation over the epicentral zone) because the directional propagation of such a waves would shift the interaction zone with the ionosphere of a thousand of kilometres far from the epicentre Yang et al., 2019;Sorokin et al., 2020).

SEISMO-ELECTROMAGNETIC PERTURBATIONS DETECTED ON GROUND
As mentioned above, beyond the mechanical oscillations (of the ground surface or of the atmospheric layers), the geochemical fluids variations and the thermal anomalies, the list of proposed earthquake-precursors on ground includes many observations of electromagnetic disturbances. We will concentrate only on measurements on ground, addressing to the companion article by Picozza et al. (2021) on the same volume for a review of the space born detection of possible seismo-electromagnetic precursors.
One of the most cited, and now "classical", case study is the observation of the increase of the ULF (0.5-2.0 Hz) magnetic field amplitude stated one month before of the Loma Prieta (California) earthquake of October 17th, 1989 (Ms 7.1) registered 7 km far from the epicentre. The amplitude of the horizontal component intensified a couple of weeks before and then further strongly increased in the 0.01-0.5 Hz range some hours before the shock, when power was lost, . Nevertheless, a detector operating 50 km far did not measure any anomalous signal in the same period. The hypothesis that the anomalous observations could have been originated by atmospheric natural background has been rejected, but on 2009 Campbell (2009) (Fraser-Smith et al., 2011), confirmed that the large-amplitude fields observed before the Loma Prieta earthquake are fundamentally different from "natural magnetic disturbance fields" and that they may well have been precursor fields to the earthquake.

Looking for ULF and Higher Frequency Electromagnetic Precursors
On ground, Fraser-Smith et al. (1990), reported a SNR in the ULF frequency range up to 60. Instead, the SNR reported for satellite observations by Larkina et al. (1989), Serebryakova et al. (1992) and Parrot (1994) was limited to 10. Nevertheless, this value is significantly higher than the background noise suggesting the possibility to detect seismo-electromagnetic precursors also from space at least in some case of strong earthquakes. The alleged preearthquake magnetic ULF (lower than several Hz) signals have a skin depth (tens of kilometres for 1 Hz) larger than that at higher frequencies (hundreds of meters for 10 kHz); moreover, low frequency magnetic signals suffer of a reduced attenuation, which allows long-range detection. These are some of the reasons that have drawn much attention to low-frequency magnetic field anomalies as important earthquake precursors. In addition to the Loma Prieta observations, early studies of ULF precursors were done also on the 1988 Spitak M 6.9 earthquake (Kopytenko et al., 1993), 1993 Guam M8.0 (Hayakawa et al., 1996), 1996 Hetian M7.1 event (Du et al., 2002), 1997 Kagoshima M 6.5 earthquake (Hattori et al., 2002), etc. In order to take into account the dependence of the ULF magnetic anomalies from the magnitude and the hypocentre distance earthquake from the observation site (Hattori 2004;Schekotov et al., 2006), Hattori et al. (2006) and Hattori et al. (2013) suggested to consider the cumulative (daily sum) of the local energy of the earthquakes weighted by the squared distance from the measurement station. This method have been adopted also by Han et al. (2014), that have carried an interesting statistical study -based on a superposed epoch approach for a long time period by analysing data from 20 geomagnetic observatories in the period 2009-2012 pointing out a relevant perturbation about one month and half before the earthquake.
In the same years, several studies have been carried out looking for possible precursors on ground. in frequency bands higher than ULF, in particular investigating the variability of well know and highly stable radio signals used for communication and geolocalization networks. By studying the attenuation of radio LF and VLF point-to-point transmissions over long distances, several authors (see for example Hayakawa, et al., 1996;Bella et al., 1998;Biagi et al., 2001a;Biagi et al., 2004) have reported cases of attenuation of broadcasted signal some days before seismic events occurred near the radio propagation path and suggested that preseismic ionospheric disturbances could have affected the conditions of transmission/reflection in the Earthionosphere waveguide. More recent investigations have claimed similar effects (Rozhnoi et al., 2009;Fidani, 2010;Biagi et al., 2011;. One of the most controversial aspect of seismo-associated electromagnetic precursors is the large variety of observations even for the measurements possibly associated to the strongest events. For example, electromagnetic precursors have been detected, or claimed, only in some specific frequency bands and the spread in latitude and in longitude is different (or even opposite with a claimed clustering) in some analyses. We believe that, at present, many aspects of the seismoelectromagnetic phenomena must still be understood. Nevertheless, the studies seem to support the hypothesis that seismo-electromagnetic precursors can be observed several hours before medium and strong earthquakes and that their intensity seems higher near the epicentre.

How Seismo-Associated Electromagnetic Disturbances Could be Generated
Several models have been proposed in literature in order to explain the observed phenomenology of seismoelectromagnetic field precursors. We can divide them in two main categories depending on whether they could explain the lower frequency precursors (ULF) (Uyeda, et al., 2009b) or generate the disturbances at higher frequency (mainly ELF/ VLF but also higher bands up to HF).

Generation of Low Frequency Magnetic Fields
In order to explain the ULF precursory fluctuations, three main effects have been suggested such as 1) the magnetohydrodynamic (MHD) effect (Draganov et al., 1991); 2) the piezomagnetic effect (Sasai, 2001) (Parrot, 1995) and 3) the electrokinetic effect (Nourbehecht, 1936(Nourbehecht, -1963Fitterman, 1979). 1) In the MHD effect, a conducting fluid moving in a magnetic field generates a secondary induced field. We can define a magnetic Reynold number Rm such as the ratio between the convective term (due to the resistance to flux change of magnetic field) and the diffusive one (due to Ohmic dissipation) in the MHD equation. From the formula R m μ 0 σ v L (where σ is the conductivity, v the velocity of the fluid and L is the characteristic length scale of the source) the induced field B i can be estimated as a function of the external one: B i R m B. 2) In the piezomagnetic mechanism, the applied stress would change the magnetization M of the ferromagnetic rocks inducing a secondary magnetic field. By solving the differential equation, that, in an isotropic material, we can write as ΔM i β T ij M j -where β is the sensitivity to the stress and T ij is a function of the stress tensor, of the displacement vector and of the material's parameters-it is possible to estimate the magnetic field at the surface due to piezomagnetic mechanism. 3) Finally, in the electrokinetic effect an electric current is generated not by an electric field gradient, but by a pressure gradient at the electrified interface of a solidliquid boundary. The electrokinetic current would be responsible of an induced magnetic field according to the Biot-Savart law. Draganov et al. (1991) reconciled the  observations with an origin in the MHD effect. Nevertheless, this conclusion was obtained by using values of permeability and lithostatic pressure of about two orders of magnitude higher than that at the Loma Prieta hypocentre depth (Fenoglio et al., 1995). On the contrary, it can be argued that, due to the strong attenuation of the magnetic field intensity, which decreases as the third power of the distance, MHD effect would give a negligible contribution to the magnetic signal detected on ground. Even the contribution of the piezoelectric effect could be negligible (Fenoglio et al., 1995), being of about 10 −2 nT, i.e. two order of magnitude less than the two or about 7 nT reported by . Whereas the electrokinetic effect could be able to generate a signal of 5-10 nT of the same order of magnitude of the signals measured before the Loma Prieta earthquake (Fenoglio et al., 1995). Interesting is also the research activity that tries to apply/ extend the lesson learned in laboratory to research on field (see for example Vallianatos and Tzanis, 1999;Tzanis and Vallianatos, 2002;Vallianatos et al., 2004).

Electromagnetic Emissions at Higher Frequency
In order to explain the generation of electromagnetic emissions detected before some earthquakes (mainly in the ELF/VLF range but also up to HF) several physical models have been proposed that can be summarized in 4 main typologies such as: contact electrification, separation electrification, piezo-electrification (Ogawa et al., 1985;Zlotnicki and Cornet 1986) and atmospheric electricity generated by radon exhalation (Pierce, 1976). The electric field emissions generated in granite samples under bending or impact/shock have been studied by Ogawa et al. (1985) that has reconciled the observed phenomenology with the electric dipole momentum due to contact/separation electrification or piezo-electrification that would induce a charge separation. At a distance r from the dipole momentum (p), the near-, induced-and radiated-field (associated with increasing frequency range) are proportional to p/r 3 , dp/dt/cr 2 and d 2 p/dt 2 /c 2 r respectively.

Most Recent Hypotheses on the Generation of Electromagnetic Fields
Some authors such as Liperovsky et al. (2008a) suggested that the variety of variables and coupling mechanisms involved in generating precursors and the large phenomenology of observed anomalies could be reconciled claiming that different physical mechanisms can explain quite the same precursors. It is not necessary to quote Occam's razor to argue that the simplest explanation, with fewer parameters, is usually the correct one. However, a consistent unified model would certainly be preferable for more stringent checks of the variables involved. By summarizing results from literature (including cited works and also Kondo, 1968;Vershinin et al., 1999;Hao et al., 2000) the amplitude of on ground preseismic quasi-static electric field disturbances-detected on long spatial (hundreds of kilometers) and temporal scales (from 10 days up to few hours)-can be estimated never exceeding 100 V/m (even though sudden and more intense increases fluctuations have been reported on shorter distances). In space, these values reduce to about 10 mV/m (Zhang et al., 2014). The proposed hypotheses of an electrostatic origin of the lithosphere-atmosphere ionosphere coupling mechanism, through a direct propagation, have been ruled out by several authors such as Sorokin et al. (2020). Moreover, Denisenko (2015) and Denisenko et al. (2018) have demonstrated that the penetrations of both a lithospheric electric field and a current in the ionosphere is negligible. In order to explain preseismic electric field anomalies (Freund, 2010) proposed that they could be originated by tectonic stress applied to lithospheric rocks. In fact, laboratory experiments have pointed out that igneous rocks under stress are able to generate currents carried by both electrons and positive holes that can cause positive electric potential, ionization of air molecules and corona discharge on the rocks surface. This model has been criticized for two different points of view. The anomalies registered in space (of for example ionospheric plasma parameters, electromagnetic values, etc.) consist of both increase and decrease of the studied variable: some time the same variable can show an increase or a drop along the temporal series about a given seismic event. According with Pulinets and Ouzounov (2011) the "unipolar character" electric field suggested by Freund (2010) would be in contradiction with the experimental results of bidirectional behaviour of ionospheric preseismic effects and could be applicable only to earthquakes occurring on ground because the release of charges in atmosphere would need of solid rocks on the surface. Moreover, even though the mechanism proposed by Freund (2010) could contribute to explain short time electric field fluctuations (on the scale of about 10 min), Sorokin et al. (2020) argued that it would not be able to explain quasi-static electric field in space on a longer temporal scale. Pulinets and Ouzounov (2011) proposed the so-called LAIC model that can be summarized as follows. Radioactive gases (mainly radon) exhalations would be responsible of the preseismic chain. In the first step, the higher radon concentration (through alpha decay) would increase the local ionization in lower atmosphere and consequently the attachment centers for water vapor condensation causing a higher release of latent heat (Garavaglia et al., 2000;Dey and Singh, 2003;Silva and Claro, 2005). The resulting anomalies in the gradient of temperature would: 1) propagate from ground up to the high atmosphere (as detected by the OLR measurements), and 2) generate a flow of clusters/condensation nuclei (Svensmark and Friis-Christensen, 1997;Svensmark et al., 2007) that will affect the clouds formation rate and morphology (Ondoh, 2004). Simultaneously, the conductivity of atmospheric layers would change affecting the Global Electric Circuit inducing fluctuations in the ionospheric parameters such as: height scale, temperature and density of electron and ions, ion composition, etc. According with the authors, the model results into: the formation of ionospheric irregularities, electromagnetic disturbances, and the precipitation of charged trapped magnetospheric particles. Finally, the model forecast a feedback process when precipitating electrons would increase the D layer ionization affecting: VLF reflection height, whistler dispersion and more in general the e. m. propagation conditions in the Earth-ionosphere. The "bipolar character" of Pulinets and Ouzounov (2011) model aims at reconciling in a consistent scheme the extended phenomenology of the preseismic observations: the positive or negative direction of anomalous electric field on the ground can induce both positive and negative variations observed for several physical parameters (such as the electric air conductivity over the earthquake, the ionospheric parameters, the cloud formation, etc.).
More recently, Sorokin and Hayakawa (2014) proposed that preseismic processes would originate a current source near to ground that modify the lower atmosphere properties by generating an EMF (electro-motive force) that will cause ionospheric electromagnetic and plasma disturbances. The model would allow explaining several aspects of the phenomenology such as: the existence of quasi-static electric field on ground and in space over the seismic regions, the increase of conductivity with altitude and a mechanism for limiting the vertical component of the electric field on ground. Associated with EMF enhancement, ionospheric instability can generate AGW, field-aligned currents and plasma irregularities along the magnetic tubes, spectral broadening of VLF transmitter signals, enhancement in space of ELF (due to scattering of with conductivity irregularities) and changes in Schuman resonance harmonics.
A common feature to the models [such as Pulinets and Ouzounov (2011) and Sorokin et al. (2020)] based on radon release is that they could be applicable both for in land and undersea earthquakes and potentially independently from the seismo-tectonic characteristics of the involved faults, as long as a sufficient radon exhalation takes place.

Ionospheric Disturbances
The ionosphere is a highly dynamic environment characterized by global structures (evolving under the solar periodic and irregular activity) as well as by local irregularities on different scales and times. The list of ionospheric parameters routinely monitored by ionosondes and GPS that have been studied in correlation with seismic activity includes total electron content (TEC); F2-layer critical frequency (foF2); electron temperature at F2-layer heights; LF radio signals etc. Some references can be for example: (Strakhov and Liperovsky, 1999;Ondoh, 2003;Ondoh, 2009;Trigunait et al., 2004;Hobara and Parrot, 2005;Liu et al., 2006;Maekawa et al., 2006;Ondoh and Hayakawa, 2006;Dabas et al., 2007;Chum et al., 2016). Hayakawa et al. (2011) reported an increase of the amplitude and frequency shift of the fourth Schuman resonance before seismic events and proposed that this could be due to changes in the height of reflection and absorption of electromagnetic waves induced by variation of the ionization degree in the D layer.
Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 676766 Ondoh (2003) reported that changes in the altitudinal profile of ionospheric plasma density over a seismically active region could be accompanied by the formation of sporadic layers in the lower ionosphere. In addition, Pfaff et al. (2005) tried to reconcile the occurrence of sporadic E layer with seismic events, proposing that this would be induced by differential charging of top and bottom sides of clouds due to the migration of radon. Silina et al. (2001), Ondoh (2004), Ondoh and Hayakawa (2006) and Korsunova and Khegai (2006) and Korsunova and Khegai (2008) have identified ionospheric earthquake precursors by using sporadic Es layer. The method applied to ionosonde data has shown an increase of Es occurrence in conjunction with strong earthquakes (M > 7) in Japan and moderate earthquakes (with 5.5 ≤ M) in Italy (Perrone et al., 2010). Similarly, Perrone et al. (2018) have identified anomalous variations of the sporadic E-layer parameters (h′Es, foEs) and foF2 in occasion of earthquakes with magnitude M ≥ 6.0 in Greece during the 2003-2015 period. By studying the foF2 variations, through ionosondes measurements, Hobara and Parrot (2005) have detected a decrease near the epicenter of the Hachinohe earthquake (magnitude 8.3) of 1968. A decrease of foF2 (greater than 25%) has been reported also by Liu et al. (2006) within 5 days before of more than 150 seismic events (of magnitude greater or equal than 5) occurred in a period of 5 years. The authors pointed out a correlation between the entity of the decrease and the earthquake magnitude, whereas the effect decreases far from the epicenter and is not more evident more than 150 km from the epicenter. In addition, Kandalyan and Alquran (2010) have investigated the existence of a possible correlation between the occurrence of strong earthquakes and ionospheric scintillations. Even though, many TEC anomalies, detected by GPS networks, have been correlated to seismic events (see also Tsai et al., 2006;Singh et al., 2009;Jhuang et al., 2010;Hasbi et al., 2011;Kon et al., 2011;Ma and Wu, 2012;Contadakis et al., 2012 and references therein) nevertheless, the debate is still intense because some authors (e.g., Rishbeth, 2006) highlighted that generally TEC perturbations could be caused by the natural ionospheric and geomagnetic variability and that the precursors analyses are generally carried out after the earthquake occurrence (e.g., Mulargia and Geller, 2003;Afraimovich and Astafyeva, 2008). For this reason, statistical analyses on long time series are particularly valuable. By analyzing changes in total electron content from the global ionosphere map (GIM) within 2.5°l atitude and 5.0°longitude around the earthquake epicenter, Thomas et al. (2017) found no statistically significant TEC changes before the 1279 M ≥ 6.0 earthquakes for the period 2000-2014 (after declustering for aftershocks). The Wenchuan earthquake on May 12th, 2008, 14:28 LT, occurred in the eastern edge of the Tibetian plateau, is remembered not only for its devastating consequences-with tens of thousands of death and the large-scale disruptions in the western Sichuan-but also for having given a further input to Chinese scientific community, involved in studying earthquake precursors, to develop the satellite CSES mission designed for remote sensing the Earth electromagnetic environment. Zhao et al. (2008) have studied TEC data from 58 GPS receivers installed nearby the epicenter of the Wenchuan earthquake. Whereas the enhancement observed on May 3rd have been mainly reconciled with geomagnetic storm occurred in the 00:00-10:00 UTC, the analysis has pointed out a TEC anomalous disturbances detected on May 9th over the southern China and its conjugate point in the southern hemisphere could be related to the Wenchuan earthquake. For specific techniques that allow discriminating between real seismo-ionospheric anomalies and artifacts see also (He et al., 2014). By analyzing VTEC measurements, during all of the 20 M ≥ 6.0 earthquakes in the Taiwan area from September 1999 to December 2002, Liu et al. (2004) reported preseismic ionospheric anomalies (with respect a 15-days running average) before of the 80% of the analyzed earthquakes. More precisely, for the seismic events in the Taiwan, ionospheric VTEC remarkably decreased during 18:00-22:00 LT within 1-5 days before the earthquakes.
Observations are in agreement with the 93% reported by Liu et al. (2000). In occasion of the Chi-Chi earthquake Liu et al. (2010) reported that the correlations between the co-located NmF2 and VTEC about one week before the event was extremely high (about 0.95) suggesting that TEC measurements via GPS receiver would allow monitoring preseismic ionospheric anomalies. Liu et al. (2004) suggested that the anomalies appear first near the Earth's surface and then extend to higher altitudes through an unknown mechanism such as diffusion, atmospheric gravity waves or due to some vertical preseismic electric field. Since a large amount of contribution comes from higher altitudes, the appearance of the diurnal VTEC features may be somewhat later than that of the NmF2 (or foF2). By investigating TEC spatial distribution, Liu J. Y. (2009) and Zakharenkova et al. (2008) reported that preseismic TEC fluctuations seems accompanied by simultaneous changes of the electric field in the ionosphere [of the order of 1-10 mV/ m, also in agreement with Klimenko et al. (2012)], but without variations on ground; the authors also estimated that an increase of several times of the aerosols concentration on ground could cause TEC variations of several tens of percent (see also Ruzhin et al., 2014).
According with Freund (2010), the rocks under stress would constitute a dynamo able to generate preseismic currents of electrons and positive holes. The author of Kuo et al. (2011) estimated that current density of about (10 -7 -10 -6 ) A/m 2 would be able to generate a daytime TEC variations from 2% up to 25%. On the contrary, Sorokin et al. (2020) highlighted that this hypothesis would be incompatible with the measurements because: 1) the duration of the Freund currents (only several minutes) doesn't match with that of the TEC variations; for achieving the claimed current intensity it would be needed a vertical electric field of about (10 7 -10 8 ) V/m that is higher than the values observed on ground in the seismic areas.
The investigation of seismo associated disturbances on a long temporal scale can be carried out by the study of the variations of amplitude and phase of VLF signals (between ground based transmitters and receivers) traveling over the preparation zone of seismic events (Biagi et al., 2004;Rozhnoi et al., 2004;Hayakawa et al., 2011). Amplitude and phase of such signalsgenerally very stable as a function of time-have shown variations of the scale of tens of minutes possibly induced by variations the height of the Earth-ionosphere waveguide due to local ionization or secondary effects induced by seismic sequence.

Multi-Parametric Analyses Looking for Earthquake Precursors
The common feature of the most recent models for investigating earthquake precursors is to assume a multi-parameter coupling mechanism (which also includes several feedback processes, see for exampleİnan Pulinets and Ouzounov 2011;Freund 2013;Hayakawa et al., 2018;Ouzounov et al., 2018) according to which the exchange of energy and particles between the lithosphere and the lower and upper atmospheric layers of the ionosphere, during the earthquake preparation phase, causes a variety of precursor phenomena that should be considered as a whole [in the so called holistic approach of De Santis et al. (2019b)]. Therefore, one of the updated trends in precursor research is the simultaneous measurement of several physical variables [such as radon and other gas release, air temperature and humidity (e.g., Akhoondzadeh et al., 2018), thermal infrared emissions (e.g., Natarajan and Philipoff, 2018), electromagnetic fields and ionospheric parameters, including electron density (e.g., Zhang et al., 2011;De Santis et al., 2015) etc.] for a joint analysis of such complex and interrelated system . Because data needed for these investigations are collected both on ground and by satellite, we address the reader to the companion article Picozza et al. (2021) and to therein references, where this topic is further discussed for space-based observations of precursors. The possible correlation between seismic events and fluctuation of water column content has been studied for example by Dey et al. (2004) for the 2001 M7.8 Gujarat earthquake; by Ma et al. (2010) for the Hengchun (Taiwan) M 7.2 event of 2006 and by Wu et al. (2016) for the 2009 M6.2 L'Aquila earthquake. Anomalies of ozone have been investigated for example by Akselevich and Tertyshnikov (1995) and by Tronin (2002a). In this framework, Piscini et al. (2017) proposed the CAPRI (Climatological Analysis for seismic PRecursor Identification) algorithm searching for anomalies of the time series of climatological parameters. The authors analyzed skin temperature, total water vapor column, and total ozone column for the period from two months prior to the entire Amatrice-Norcia seismic sequence (central Italy)-which began on August 24, 2016 with an M6 earthquake and included two other large quakes (i.e., an M5.9 earthquake on October 26 and then an M6.5 on October 30)-comparing the measurements with time series of data from the previous 37 years in order to remove the possible effect of global warming. The method was also extended to analyze the M 6.2 L'Aquila (2009) earthquake. Piscini et al. (2017) reported persistent anomalies that emerged simultaneously in all parameters analyzed, conclusions that would be confirmed by comparing the data with those of the same months in other seismically quiet years. On the other hand, the authors have commendably pointed out that a single anomaly, if present, of an individual climatological parameter would have low statistical significance as it could be caused by several sources unrelated to the earthquake. Moreover, weather phenomena could move and mix the substances released from underground into the atmosphere could mask precursors, if any, lowering their signal-tonoise ratio and shortening their persistence in the atmosphere (Marchetti et al., 2019). Similar results have been obtained by De Santis et al. (2020) that have applied a multi-parametric study to the 2019 M7.1 Ridgecrest earthquake, by analysing furthermore methane exhalations [as suggested by Cui et al. (2019)], electron density fluctuations, and magnetic field measurements from Swarm satellites (Friis-Christensen et al., 2006;Zhu et al., 2019). The results can be summarized as follows: 1) precursor times would be much longer than those identified by other papers (especially about ionospheric precursors which seem to occur only a few hours to days before large seismic events) (see for example Heki, 2011;He and Heki, 2017;Yan R. et al., 2017); 2) the preparation phase would be much longer than few days [as also suggested by Liu et al. (2020), Marchetti et al. (2019), Marchetti et al. (2020), Sugan et al. (2014), Giovambattista and Tyupkin (2004)]; claimed precursors would follow the empirical Rikitake (1987) law, recently confirmed for ionospheric precursors from the satellite by De Santis et al. (2019a). In any case, it is a complex task to monitor several or all atmospheric effects due to the LAIC mechanism that also include linear cloud structures (the "earthquake clouds" that would repeat the shape of the tectonic structure in the sky (Jones and Stewart, 1997;Nissen et al., 2012); OLR [Ongoing Longwave Radiation-infrared emission at 10-13 μm recorded above the clouds ]; jet streams (Wu, 2007); etc. Even in a multi-parametric approach, the benefit of analysing the OLR-measured in the upper atmosphere -is that the OLR can take into account the cumulative effect of all thermal contributions possibly due to the earthquake, between the ground surface and the tropopause, according with the so-called "synergy of precursors" (Pulinets, 2011). Moreover, thanks to the long time series of data collected by the NASA Aqua and NOAA/AVHRR satellites, it is possible to compare the local and temporal variations of the OLR with a wellestimated reference background that would improve the statistical significance and reliability of pre-earthquake fluctuations of this parameter. Ouzounov et al. (2007) reported that a few days before the M 9.1 Sumatra earthquake of December 26, 2004, the OLR value was above 80 W/m 2 compatible with the estimate given by Kafatos et al. (2007). In general, the results of the multi-parameter analyses-also made possible by recent big data analysis techniques and large computational capacity-reinforce the idea of considering such an integrated anomaly recognition system as an effective tool for systematically finding earthquake precursors.

POSSIBLE APPROACHES FOR RECOGNIZING NATURAL AND ARTIFICIAL ELECTROMAGNETIC BACKGROUND
In order to identify electromagnetic earthquake precursors it is essential to reject carefully the background due to natural nonseismic sources and industrial electromagnetic noise. Although, the main driver of the electromagnetic Earth environment is the solar wind that shapes the magnetosphere and affects its dynamics, several other natural and artificial sources of electromagnetic emissions take place in the geomagnetic cavity. The terrestrial electromagnetic background noise of not-ionizing radiations (distributed-although not uniformly-in a broad range of power levels and frequencies from about mHz up to 300 GHz), is generated by a wide variety of natural (atmospheric thunderstorms (e.g., Neubert et al., 2008) or cosmic rays) and man-made (power line harmonic radiation, industrial and communications equipment, etc.) sources (Hayakawa, 1994). These emissions constitute a significant background with the respect to the elusive electromagnetic emissions and related ionospheric disturbances possibly caused by natural geophysical activities, such as earthquakes and volcanic eruptions .

Natural Non-Earthquake Electromagnetic Noise on Ground and in the Near Terrestrial Environment
Depending on the frequency, the natural electromagnetic noise in the magnetosphere arises from wave-particle interaction phenomena, while within the ionosphere originates from atmospheric electric discharges through several propagation mechanisms Simões et al. (2012). On ground, in some areas, the artificial noise generated by technological devices exceeded the natural one, while in the VLF-HF range, the atmospheric electromagnetic emissions exceed the artificial ones in rural areas (in the order of some tens of dB) and are comparable to them in industrialized zones. At higher frequencies, the electromagnetic noise induced by thunderstorms becomes less important, and that of cosmic origin (up to millimeter wavelength) prevails (Bianchi and Meloni 2007).
In the ULF band (from 1 mHz up to 1 Hz, i.e., from the lowest frequency that the magnetospheric cavity can support, up to the ions gyro-frequencies) are observed the geomagnetic pulsations and the "incoherent noise" (Lanzerotti et al., 1990;Jacobs et al., 1964). The low-frequency pulsations are generally originated by the Kelvin-Helmholtz instability in the magnetopause, through to the interaction of the solar wind with the magnetosphere, or by the waves "upstream" in the "foreshock" region. Wave-particle interactions in the magnetosphere give rise to the so-called "chorus" and "hiss" phenomena. The chorus are among the most intense electromagnetic emissions generated in the external magnetosphere and propagating up to the Earth surface where are observed at intermediate latitudes. The spectral characteristics of the chorus (from about 500 Hz to about 1.2 kHz) consist in a succession of tonespredominantly growing-that resemble to the birds chirping, that constitutes the origin of their name. Hiss are intense electromagnetic emissions, occurring mainly in the auroral zone, in a broad frequency band (from a few hundreds of Hz up to several tens of kHz). The most important electromagnetic phenomenon in the ELF band is the Schumann resonance, at about 7 Hz and higher harmonics, arising at the proper oscillation frequencies of the natural Earth-ionosphere waveguide (Polk, 1983). The tropospheric layers (with variable electrical conductivity) between the Earth crust and the ionosphere (both of them schematized as perfectly conductors) constitute an electromagnetic cavity in which electromagnetic radiation is trapped and waves can propagate. Waves constructive interference can excite resonances in the frequency band of about 6 ÷ 60 Hz in the above-mentioned Earth-ionosphere waveguide. Lightning are the main source of electromagnetic noise background in the ionosphere, where they generate emissions from ELF (about few Hz) up to VHF (about hundreds of MHz) although most of the energy is concentrated in the VLF band (from 0.1 to 10 kHz) with a typical power spectrum slope. Some thousands of storms are estimated to occur daily on the Earth (Shvets and Hayakawa, 2011), generating about 100 lightning per second, with discharges up to 10 kA per each and releasing an amount of energy from a few units up to a few tents of GJ, i.e., powers of the order of 0.1 ÷ 1 TW, for a total of 10 19 J released yearly around the world. The natural electromagnetic phenomena relevant in the ELF-VLF frequency bands are the so-called sferics, tweeks and whistlers, generated by lighting electromagnetic pulsed signals (of a few ms), that travel with a low attenuation in the Earth-ionosphere waveguide for thousands of kilometres (Helliwell, 1965;Hitchcock and Patterson 1995). The propagation of the spherical waves is determined by the variable ionospheric conditions and the sferics can be observed with directional antennas and AM receivers as typical disordered sounds called "statics". The tweeks (generally detected in the evening after the sunset) are sferics showing a spectral dispersion during their propagation and that, in acoustic spectra, sounds similar to birds singing with frequencies of about 1 ÷ 7 kHz. The higher harmonics penetrate in the ionosphere more in depth than the lower components that, being less attenuated, cover longer distances. Different paths imply different arrival times and the spectrograms show descending tones with a duration from about 25 ms up to about 150 ms. The whistlers are intense circularpolarized electromagnetic waves propagating at a frequency below the plasma frequency (from about 6 kHz up to a few hundred Hz), generated by lightning and perceived in radio receivers as characteristic decreasing tones. Whistler waves rotate clockwise propagating along the geomagnetic field lines through the ionosphere up to the magnetosphere, bouncing back and forth, and showing a significant spectral dispersion as function of the path length and the characteristics of the crossed medium. Also in the LF/MF/HF bands, up to the plasma frequency, the natural electromagnetic noise is mainly generated by atmospheric phenomena (with an amplitude decreasing at higher frequencies) and its propagation is affected by local ionospheric conditions and geometry of the paths. For frequencies higher than 15 ÷ 30 MHz the cosmic noise of astrophysical origin appears (Kraus, 1988;Erickson, 1990).

Anthropogenic Sources of Electromagnetic Emissions in the Near-Earth Space
The natural electromagnetic background is accompanied by the emissions due to human activities. The artificial noise originated by industrial technologies (power lines, radio and TV broadcasting stations, communications facilities, etc.) strongly depends on the distance from the sources, can vary also by many orders of magnitude, frequency and power, and show distinctive features such as continuous or impulsive regime, modulation and polarization. In the ELF band, the most powerful source of artificial noise (except Antarctica where the natural ELF Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 676766 emissions can be detected with few anthropogenic interference) is represented by power lines harmonic radiation (PLHR) which ideally operate at a single frequency of 50 Hz (60 Hz in the U.S.). At these frequencies, the electric and magnetic fields are virtually disconnected and (due to the multipolar configuration of the electric lines) decrease dramatically by increasing distance from the lines. Satellite observations show that PLHR can contribute-through non-linear interactions-to precipitation of Van Allen particles from the slot region in the radiation belts. The principal part of the PLHR energy dissipates in the lower ionosphere modifying ionospheric currents. Ground-based radio-navigation and communications VLF transmitters in the 10 ÷ 20 kHz frequency band can trigger new waves, ionospheric heating, wave-electron interactions, and particle precipitation via the cyclotron resonance. The wave and the particles interact when the Doppler-shifted wave frequency seen by the particles is close to the electron gyro-frequency. The trajectory of particles follows the magnetic field lines but the ray path of the injected waves is only field aligned if the wave is propagating in a whistler mode duct. The strongest interaction region is around the geomagnetic equator. Triggered emissions by coherent waves are related with non-linear wave growth caused by resonant particle trapping in a non-uniform magnetic field. At HF frequencies, powerful broadcasting stations can provoke ionospheric Joule heating by changing plasma temperature and density Erickson (1990). The dissipation in the ionosphere of the above-mentioned waves may contribute to the global warming of the Earth, since the change in global temperature increases the number of natural lightning discharges in the atmosphere and this produce more magnetospheric whistlers that may provoke heating and ionization in the lower ionosphere. This is a feedback mechanism since lightning are sources of NOx that influence the ozone concentration in the atmosphere contributing to the greenhouse effect. Moreover, precipitation of energetic electrons by anthropogenic waves may trigger other lightning discharges. This explains the importance in studying such anthropogenic electromagnetic emissions. For VHF or higher frequencies, that are relevant for radio astronomy observations Erickson (1990), the anthropogenic electromagnetic emissions are due to radio and television broadcasting stations, mobile communications, car ignition systems and industrial equipment, while radar and satellite devices, as well as highly directional SHF/EHF emitters, do not give a significant contribution to the electromagnetic noise background.

Rejecting Non-Earthquake Associated Effects From Data Analysis
The described features could allow distinguishing the effects of natural non-seismic, seismic and man-made electromagnetic emissions. In some frequency range (i.e., at 50 or 60 Hz or for signals from radio transmitters), it is easier to reject man-made electromagnetic emissions. On the other hand, in many cases, it is not straight to remove the background emissions due to magnetospheric disturbances and tropospheric sources aimed at identifying seismo-electromagnetic disturbances. This is particularly true for example in the statistical analyses (such as those based on the epochs overlap method) where data from different geographic regions (naturally exposed to different geomagnetic processes) are grouped and studied together. In this framework, the analyses published in literature have adopted different strategies. Unfortunately, in many cases, the studies have not carefully excluded data collected in geomagnetically disturbed periods. A reduction of the background can be achieved by using magnetospheric and ionospheric indices. Some examples are the well-assessed (global or aerial) indices Kp, Ap, AE, Dst, etc. released with relatively low time resolution (hour or multi-hour scale). Other ones are the more detailed indices able to describe the geomagnetic perturbations at mid-latitudes as a function of longitudinally asymmetric (ASY) and symmetric (SYM) disturbances for both H and D components (respectively parallel and perpendicular to the geomagnetic dipole axis). A further filtering, more difficult to be applied but that can be investigated, is to take care of the metrological conditions and of the lightning activity in the area of the earthquake and in the conjugated one. Frequencies range around the band of the known VLF transmitters should be also rejected in the analyses of seismo-associated disturbance because the power of the artificial signals can overwhelm the natural and faint emissions possibly associated to earthquakes. On the other side the high stability in power and frequency of these transmitters, as well as of the radio broadcasting stations, have been and are investigated in order to point out fluctuation in the transmission parameters possibly due to seismo-electromagnetic disturbances occurred along the transmission path from the transmitter and a fixed receiver station. In addition, the polarization of eventually directional signals can give hints on the internal or external origin of detected measurements candidates to be associated to lithospheric processes. In general, a strong interdisciplinary approach, between different research fields, is mandatory to clean data, to reject background and to better distinguish such complex phenomena. Hattori (2004) applied the principal component analysis (PCA) for decomposing and filtering time-domain series of observations in order to identify background signals (generated by DC noise, natural emitters and man-made devices) that could be removed in order to point out earthquake related anomalies (if any). The authors have applied the procedure to ground measurements of magnetic field (also including data gathered during strong geomagnetic disturbances) and reported that the first two principal components would be correlated to geomagnetic variation and anthropogenic sources, respectively, while the third component has been investigated in order to point out seismo-associated disturbances. By applying principal component analysis to magnetic data from six observatories (4 near Napa, California, together with two remote reference stations), Kappler et al. (2017) were able to identify and distinguish global geomagnetic signals (such as solar-generated noise) and anthropogenic signals.
Recently some authors have explored the application of machine learning methods for automatically classifying and recognizing earthquake precursors on ground and in space (see for example Rouet-Leduc et al., 2017;Li et al., 2018;Akyol et al., 2020;Johnson et al., 2021;Xiong et al., 2020;Xiong et al., 2021). The field of research is very interesting, but asks for some caution. Even though the amount of available observations is huge, the supervised methods could suffer the difficulty of learning from a limited number of tagged measurement: because even for scientists the signature of possible precursors is still unclear, it is hard to define a clear "learning path" for an automatic recognition. On the other side, the unsupervised methods could help in extracting/clarifying the signature of precursors, but the artificial intelligence systems should be trained also in distinguishing spurious phenomena for the needed background rejection.

REPORTS OF EARTHQUAKE LIGHTS FROM WITNESSES AND ANOMALIES IN ANIMAL BEHAVIOUR
One of the most controversial debates about earthquake precursors is the presumed correlation between the occurrence of seismic events and: 1) observations of earthquake lights; 2) reports of anomalous animal behaviour. We highlight that these two groups of alleged precursors are significantly inhomogeneous and very different from each other. Nevertheless, we present them in the same paragraph not with the purpose of mixing such a broad phenomenology, but because they often share a noninstrumental method of observation (with some recent, but still rare and unclear exceptions of photo/video recording) and suffer from the high variability of sensory perceptions by biological organisms with the consequent inherent difficulty of assessing the statistical significance of the reports.

Reports on Earthquake Lights
Pre-and co-seismic visible luminescence phenomena (so-called earthquake lights, EQLs for short) have been reported by several authors such as Galli (1910), Terada (1930), Richter (1958) and Yasui (1968) that published the first photograph of EQL; Derr (1973), Derr (1986), Tsukuda (1997) for the Kobe earthquake; Papadopoulos (1999), Stothers (2004), St-Laurent (2000 and Omori et al. (2007). The reports refer primarily to shallow earthquakes of high magnitude, but in some cases, observations have also been collected on medium events with deeper epicentres. EQLs have been reported from a few up to several kilometers from the epicenter, across the visible spectrum, with durations ranging from a fraction of a second to several seconds. EQLs have been sighted both before (from several weeks to a few seconds, near the epicenter) and during earthquakes (even far from the epicentre). More recently, in his valuable catalog of evidence of eyewitnesses -ranging from 9 months before up to five months after the main shock of the 2009 L'Aquila earthquake-Fidani (2010) compiled a careful classification (and an attempt to locate) observations of lights, flames and other bright observations-reported in a wide area up to 50 km far the city of L'Aquila and along the Aterno valley-looking for a correlation with the occurrence of the earthquakes of the L'Aquila seismic sequence. Heraud and Lira (2011) reported some EQLs that would have occurred during the 2007 Pisco earthquake-also supported by a video recording, and apparently not induced by thunderstorms or electrical faults-suggesting that they would have been generated by some electrical phenomena not yet understood in the atmosphere above the epicentral area. In general, some of the EQLs reported in the literature seem less convincing, however, some characteristics of the EQLs (shapes, colours, flames, etc.) are recurrent in the evidences (such as in the L'Aquila, Saguenay, and Pisco events) and the quantity of observations ask for a careful consideration of this phenomenon . Persinger, (1983) suggested that co-seismic EQLs would be related to energy release during fracturing of rocks; (Brady and Rowell, 1986;Kato et al., 2010;Martinelli et al., 2020b) and other authors showed that under high-stress conditions rocks could emit electromagnetic radiation prior to fracture. Additional hypotheses have been proposed in order to explain the origin of EQLs (Mizutani et al., 1976;St-Laurent et al., 2006), including a generation due to gas and aerosol exhalation Liperovsky et al. (2005) and a high-frequency electromagnetic origin Liperovsky et al. (2008b). Analysing the observation of claimed EQLs potentially associated to about 65 earthquakes (38 from Europe and 27 from the Americas) in different geotectonic settings, Thériault et al. (2014) proposed an updated model for the origin of EQLs associated with both intraplate and interplate earthquakes, which is based on the generation of electron charge carriers under high-voltage conditions (see also Freund et al., 2009;Freund 2010). Their thesis is that EQLs may be predominantly associated with intraplaque earthquakes within or nearby rift-related structures. EQLs do not always appear to occur before all strong earthquakes, and they vary in rate. Freund et al. (2021) suggested a solid-state physics mechanism that could take into account the variety of phenomenology. The peroxy defects in igneous rocks-mainly in gabbroic rocks that fill the sub vertical dykes and would be preferentially located along boundaries or between adjacent mineral grains-would make them highly susceptible to ever so slight displacements of mineral grains. The propagation of seismic waves would activate peroxy bonds generating charges displacement. The co-seismic EQL would be caused by the rupture of peroxy bonds, a discharge from the top of the dyke, removing some of the charge carriers. Less evident is how the other processes (corona discharges, thermal infrared emissions, air ionization as well as ion and electron fluctuation in the ionosphere, and electro-magnetic anomalies) would be generated by such mechanism before the propagation of seismic mechanical waves. Within the already highly debated topic of earthquake precursors Hough (2016), the complexity of the observations of pre-and co-earthquake EQLs is even more controversial because analysing eyewitness observations is much more difficult than studying instrumental measurements. Moreover, the possibility of misinterpretations is not negligible, also because the survey of witnesses is generally based on questionnaires filled out and collected after tragic events in very difficult and stressful environmental conditions. Therefore, the statistical completeness of anecdotic or sparse datasets could strongly influence the interpretation of the collected data, the over-or under-estimation of outliers and the rejection of spurious cases. In addition, at least in some cases, data are collected and analysed much later than the seismic events, when eyewitness accounts may no longer be independent. Photographic and video evidences are certainly a valuable tool to support the reliability and objectivity of claimed observations, but published testimonies leave room for uncertainty and interpretation. Despite the large amount of testimony collected in published about EQLs (see also Lockner et al., 1983;Johnston, 1991;Whitehead and Ulusoy, 2015 and reference therein), it appears that the rate of observation by witnesses for EQLs is significantly lower than even the most conservative estimates of hallucination prevalence in the normative population. This consideration is neither a bias due to the authors' skepticism nor a preclusion against the investigation and existence of EQLs, but something that should be cautiously taken into account when assessing the reliability of this intriguing phenomenon. Indeed, it is worth noting that the prevalence of hallucinations is by no means negligible even in the general population (which may include psychotics, whether they know they are psychotic or not. Some estimates of prevalence in nonclinical samples are about 10% (Sidgwick, 1891;West, 1948;Posey and Losch, 1983;Asaad and Shapiro, 1986;Sidgwick et al., 1994). Tien (1991) suggested 10-30 cases per 1,000 people per year while for Ohayon (2000) (although the methodology adopted is open to criticism) up to 40% of the 13,000 subjects included in the study would have experienced daytime hallucinations, demonstrating how hallucinations can occur sporadically even in healthy subjects in normative populations (including all possible types of dis-or mis-perception). More recently, Temmingh et al. (2011) suggest a lifetime prevalence of 10-15% for vivid sensory hallucinations. Because of the difficult conditions of a post-earthquake survey, reports published in EQLs include limited information on some key witness parameters (such as education, occupation, alcohol and drug use, and cognitive status) that play a role in estimating the status of participants in psychological studies (e.g., Badcock et al., 2017;Eaton and Kessler, 1985), and that should, at least, not be ignored when collecting and evaluating data from witnesses of earthquake observations. Memories of witnesses and therefore a precise time location of the observations are key points in distinguishing EQL occurred during main shocks or before them (but possibly during foreshocks). If the co-seismic intense flashes of light bursting out of the ground could be "more easily" reconciled with the propagation of seismic waves (Freund, 2019), the preseismic luminous events are more debated, even though the supporters of EQLs claim a common origin of pre-and co-earthquake visible emissions. Therefore, filtering psychological uncertainty and the future instrumental proofs could help in cross-checking eyewitness testimonials and clarifying EQLs existence and origin. The need to consider all possible biases, including possible "human illusions," in the difficult task of distinguishing EQLs from background noise is relevant not only to reject spurious events, but also to estimate the possible rate of EQLs. According to Papadopoulos (1999), De Ballore (1913 and Mallet (1855) and Thériault et al. (2014), EQLs would occur in about 10% of earthquakes (5-6% at night, while they would be hardly visible during the day) but for Persenger and Derr (1984) this percentage would represent a lower limit because many observed EQLs would never be reported and published in the scientific literature, or are interpreted as unidentified flying objects (UFOs) (Devereux et al., 1983). All these aspects should be carefully evaluated in surveys (including those aimed at estimating EQLs) in populations exposed to long and highly stressful conditions-such as the seismic sequence that often precedes main earthquakes-which obviously have an impact on the psychological and psychiatric response to environmental stimuli (individual and social) possibly varying the threshold of reaction as well as the sensitivity of witnesses, even if sincere and in good faith! Can Anomalies in Animal Behaviour be Reliably Correlated With Impending Earthquakes?
In the letter of 1949, recently published Dyer et al. (2021), in reference to the work of Karl von Frisch (Nobel Prize in Physiology in 1973) and the sensory perception of animals, A. Einstein wrote "It is conceivable that the investigation of the behaviour of migratory birds and homing pigeons may one day lead to the understanding of some physical process that is not yet known". A few words that show an insight into animal ethology and physiology that preceded by more than half a century the investigations and discoveries of our years about some animals capabilities (Wu and Dickman, 2012;Lambert et al., 2013;Mouritsen, 2018). But, can birds be sensitive to earthquakes as Yosef (1997) has suggested? It is well know that birds can use the Sun, stars, and magnetic field to orient themselves; young and adult turtles as well as salmon navigate using the geomagnetic field as a reference system (Lohmann, 2007;Lohmann and Lohmann, 2019); glass eels would have a magnetic compass linked to the tidal cycle (Cresci et al., 2017); cellular autofluorescence is sensitive to the magnetic field (Ikeya and Woodward, 2021); etc. Although, such research could also shed new light on the long-standing claim that animals are sensitive to earthquake precursors (see for example Ikeya, 2004;Freund and Stolc, 2013;Yamauchi et al., 2014;Grant et al., 2015 and therein references) however, this topic is highly debated to the point of overt skepticism even within the most open and convinced scientific community to investigate precursors. In their famous paper, Woith et al. (2018) analysed more than 700 reports of claimed correlation between earthquakes and "anomalies" in animal behaviour. Observations have been claimed for organisms belonging to more than one hundred different species-ranging from deep-sea fish (Orihara et al., 2019) to catfish (Musha, 1957), from mice (Ikeya et al., 1996) to elephants (Garstang, 2009), from pets to snakes (Tributsch, 1982), from cows (Fidani et al., 2014;Yamauchi et al., 2017) to mouse (Yokoi et al., 2003), etc.-published in nearly 200 papers. Hypotheses advanced to explain the claimed anomalous animal behaviours range from a high sensitivity (postulated but not demonstrated) of animals to earthquakerelated mechanical oscillations to their claimed ability to sense magnetic field fluctuations of a few nanotesla ) (which should be cautiously verified, pre-earthquake magnetic anomalies being detected even below this threshold). The claimed observation distance varies from a few up to Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 676766 hundreds of kilometers, but most of them were observed within 100 km of the epicentre and about 70% within 50 km. Although the advance in time varies from a few minutes to months before the earthquake, the number of claimed observations increases close to the seismic event and almost 60% of the cases fall in the last 5 min. Observations of claimed abnormal animal behaviour have been reported for several earthquakes worldwide, but more than 50% of the reports relate to only 3 earthquakes out of the 160 analysed by Woith et al. (2018). Beyond the specific reports, the objective difficulties in analysing data (and screening publications) on putative earthquake precursors based on animal behaviour can be summarized as follows: 1) the common bias of categorizing animal behavior as anomalous events only ex post seismic events; 2) the often unclear distinction of normal/anomalous behavior when a threshold and quantitative criteria are missing (preventing a clear definition of the signal-to-noise ratio); 3) the frequent cases of (too) short or partially published time series [focusing only on the "relevant" (claimed) portion of the "anomalous" data]; 4) the partial or missing monitoring of environmental variables (such as meteorological parameters, moon phase illumination, ethological constraints/variabilities, human or predatory conditioning, animal health conditions, etc.) that may influence animal behavior even under normal conditions; 5) statistical uncertainty in looking for recurrence/exceptionality when comparing too short time series of biological cycles with external sporadic events. Some attempts have been made in the past to evaluate (systematically and independently) the reliability of some earthquake prediction approaches and methods, such as the IASPEI initiative (Wyss, 1997;Wyss and Dmowska, 1997), or the international Collaboratory for the Study of Earthquake Predictability (CSEP) (Jordan, 2006;Michael and Werner 2018), etc. However, regardless of the initiatives, these approaches cannot be fruitfully applied to control animal observations because of the incomplete description of the original observations, the lack of data, and the difficulty of repeating observations under similar and controlled seismic and environmental conditions (Fidani, 2013;Fidani et al., 2014). The combined effect of too short time series, the periodicity of biological cycles, phase shifts, and the temporal distribution of earthquakes in long seismic sequences (with the difficulty of distinguishing foreshocks and aftershocks) can lead to misinterpreting random coincidences for positive correlations. Indeed, for example, the 80% success rate of catfish in alerting for an incoming earthquake claimed by Hatai and Abe (1932) should be interpreted more carefully because on 85% of the observation days an earthquake occurred by chance. Similarly, the claimed anomalies in toad behaviour (Grant and Halliday, 2010) were observed in too short a time period compared to the overlap between their life cycle and the long sequence of the L'Aquila earthquake. The correlation between earthquake and anomalies in ant behavior suggested by Berberich et al. (2013) was not confirmed in the longer study, up to two years, by Apostol et al. (2014). And the list could be longer. Furthermore, modeling the reasons that might have driven the emergence of animal sensitivity to earthquake precursors is even more complex because it is unclear what ethological stressor might have driven the natural selection of individuals capable of recognizing earthquake-related phenomena. Even less clear is why and how this sensitivity appeared in many species (of different phylum, class, order, genus, and species) and living in such diverse environments (on land, in the sea, and in the air). Indeed, although earthquakes are disruptive to human constructions, their occurrence is relatively infrequent (even in the zone of maximum seismicity) and generally safe for animals, so it is difficult to justify that mechanical effects in the destruction of burrows or eggs can explain an adaptive response and (generalized) modification of the genome and/or ethology of animals. The hypotheses advanced can be summarized into two main groups. Adaptive evolution may have developed/increased the animals' sensitivity to early mechanical oscillations of P-and S-waves (Wikelski et al., 2020). This hypothesis has the advantage that the presumed increased ability would not necessarily be a specific response to seismic stimuli, but could simply be an increased sensitivity to noise and infrasonic waves, naturally evolved as a reaction to predator pressure (Kirschvink, 2000). On the other hand, it has been suggested that animals would be highly sensitive to fluctuations in earthquake-related parameters [such as magnetic field, humidity (Tichy and Loftus, 1996), temperature, etc.]. Although, by altering the magnetic field it has been experimentally demonstrated that it is possible to change the direction of flight of birds within a Faraday cage, it is not correct to infer the statement that all birds or animals in all circumstances are equally sensitive to the magnetic field (Kirschvink, et al., 2010). For example, even in the pigeon-which is not a migratory bird but tends to return to the dovecote-the orientation system is more complex and not yet fully understood. Mora et al. (2004) suggested that iron-rich cells in the beaks of pigeons were nerve cells containing magnetite and thus able to aid navigation through the Earth's magnetic field. However, Treiber et al. (2012) and Treiber et al. (2013) realized that the iron-rich cells are actually immune system cells (macrophages) and not neurons. Magrophages could probably have a function in orientation toward the dove, but not in orientation with respect to the magnetic field. Furthermore, while sensitivity to the quasi-static geomagnetic field, which allows animals to map and recognize their environment, may support orientation as well as migration, the postulated sensitivity of animals to even the highest frequency content of electromagnetic spectra (which could be relevant to earthquakes) is less convincing. More recently, some complex scenarios have been suggested about a putative mechanism of magnetoreception by electromagnetic induction in the inner ear of some birds (Nimpf et al., 2019). From this perspective, nonstatic magnetic field sensitivity would not be an evolutionary capacity driven by seismic events, but a collateral/derived capacity consisting of a lowered response threshold. To explain the biological effects of weak magnetic fields, some molecular transduction mechanisms have been proposed (Binhi and Prato 2018;Bialas et al., 2019). While for animal navigation/orientation, the main hypothesis is a specialized magnetic sense associated with pairs of radicals located in the Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 676766 retina of the eye, nonspecific effects could occur due to the interaction of magnetic fields with the magnetic moments of rotating molecules dispersed in the organism. Indeed, Binhi and Prato (2018) have shown that the precession of the magnetic moments of these rotating molecules can be slowed due to a mixing of the quantum levels of magnetic moments (LMMs) inducing a magnetic field dependence that is in good agreement with experiments in which biological effects arise in response to the reversal of magnetic field orientation. Although these studies would suggest a (potentially "common") sensitivity of biological organisms to the magnetic field, even under non-static conditions, nevertheless, at present, the purported response of animals to earthquake precursors seems a bit of a blanket statement. Indeed, such sensitivity-even greater than the instrumental sensitivity of measurements made in various test campaigns to study earthquake precursors-would cut across so many different species that they do not share other important sensory characteristics. Systematic monitoring campaigns with continuous bio-logging of animal collectives, including movements and physiological parameters -such as the experiments about ultra-sensitive measurement of micromovement of cows in the stable and on pasture (see Brown et al., 2013;Wikelski, et al., 2020), and therein references-could yield valuable insights on this topic.

CONCLUSION
The hypothesis that before earthquake it could be possible to detect surface deformations is supported by the observation of dilatancy effects in labs experiments before rocks rupture. GPS and satellitebased SAR interferometry have given powerful tools for such a worldwide investigation both on local scale and on the continental one through the measurements carried out along the plate boundary. However the growth of the observed effects that would prelude to the earthquakes would not allow defining time and place of the occurrence. All the physical models claimed to explain precursors share the hypothesis that fast and non-linear processes in the rocks along the seismic fault (such as deformation, dilatancy, fluids flow changes, pores volume variation, etc.) could originate the anomalous variations of observed parameters. Even though experiments and theoretical models give some hints about the possibility to reconcile the observed phenomenology with physical mechanisms operative in the active faults-and nonlinear effects have been detected in laboratory tests-the scalability of results obtained in small-size labs experiments up to the large scale of real seismic fault dynamics is highly debated mainly for the large uncertainty about the values of many key parameters that are still barely known in the depth conditions of the faults.
The proposed models of a direct propagation of an electrostatic disturbance from ground to the upper atmosphere and then in the ionosphere have been ruled out by theoretical calculations and simulations. In this framework, it has been suggested that the lithosphere-atmosphere-ionosphere coupling could be due to atmospheric processes including acoustic gravity waves. It worth noting that the complexity of the observed phenomenology-and mainly the high variability of the spatial and temporal scales of the ionospheric preseismic fluctuations-could be hardly reconciled with only one single coupling mechanism prevailing over all the other proposed models. On this basis, the hypothesis that the lithosphereatmosphere-ionosphere coupling is implemented by multiple physical mechanisms together seems reasonable. For example, the variation of geochemical species (mainly radon exhalations) proposed in order to explain atmospheric thermal anomalies could also generate atmospheric oscillations that can trigger AGW.
Even though the spatial distribution of preseismic geochemical fluids variations (groundwater level, vapour emissions, gases releases and radioactivity fluctuations) would be correlated with the geographical system of seismic faults, the timing of the anomalies could still be basically random as well as the earthquakes themselves and their foreshocks, generating a <<mosaic of precursors>> (Sorokin et al., 2020) function of depth and magnitude.
Up to now, in several proposed models, a key role seems played by radon emission that would vary atmospheric parameters (such as conductivity) inducing a reaction in the global electric circuits. From one side, this prospective avoids the difficulties of a direct propagation of electromagnetic signal. However, from another side, probably, these hypotheses, if confirmed, expose the research of earthquake precursors to further conceptual and practical difficulties due such as: the complexity of transport and diffusion processes; the effect of convection and turbulence in the atmosphere; the difficulty in distinguishing seismo-associated signals from the large electromagnetic noise due to thunderstorm electricity, artificial signals and geomagnetic activity. However, achieving a better understanding of the physics of earthquakes (before, during, and after the seismic event) deserves the efforts being made by the involved scientific community worldwide.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

ACKNOWLEDGEMENT
This work was supported by the Italian Space Agency in the framework of the "Accordo Attuativo 2020-32.HH.0 Limadou Scienza+" (CUP F19C20000110005).