Potential Earthquake Proxies from Remote Sensing Data

Abstract

At present, there is no clear scientific evidence of reliable earthquake precursors that can be used as an early warning system. However, many studies have also reported the existence of faint signatures that appear to be coupled to the occurrence of earthquakes. These anomalies have traditionally been detected using data from in-situ sensors near high-seismicity regions. On the other hand, remote sensors offer the potential of large spatial coverage and frequent revisit time, allowing the observation of remote areas such as deserts, mountains, polar caps, or the ocean. This chapter revises the state-of-the-art of the understanding of lithosphere–atmosphere–ionosphere coupling. It also presents recent studies by the authors’ ongoing investigation on short-to-midterm earthquake precursors. The Earth observation variables discussed are (1) surface temperature anomalies from thermal infrared or microwave radiometer measurements, (2) atmospheric signatures, (3) ionospheric total electron density fluctuations or scintillation measured from GNSS signals, and (4) other geophysical variables, including geomagnetic field fluctuations, changes in the Schumann resonance frequency, or low-frequency electromagnetic radiation. However, despite the seismic hazard risk models that exist and the results shown by these studies, it is still very difficult to predict the occurrence of earthquakes.

1. Introduction

Earthquakes are natural disasters that cause devastating damages and significant human casualties. Over the past 50 years, earthquakes have been the first cause of death from natural hazards [1]. Despite the thousands of earthquakes occurring annually, only a tiny fraction is observable. Between 1998 and 2018, earthquakes caused 846 thousand deaths and approximately US$661 billion in economic losses [2]. While earthquakes are unavoidable, efforts have been made to mitigate and restrict their impact [3].

Forecasting earthquakes, especially the significant ones that cause severe destruction, has garnered widespread attention [4]. Numerous studies have explored potential earthquake precursors, aiming to understand the physical characteristics of seismic events and their interaction with the environment [5]. Despite the continued efforts, earthquake forecasting remains challenging, as reliable indicators of future earthquakes are infrequently apparent in seismic analysis. A clear physical link between the observed anomalies and the earthquakes is needed to make the correlation analysis more robust [6].

Remote Sensing (RS) has been widely employed in earthquake studies due to its capability to provide short revisit times and large coverage [7], continuity, and more consistent data than point measurements [8].

1.1 Remote sensing observables

There are different physical mechanisms involved in potential seismic precursors. First, the frictional forces involved in the preparation and occurrence of an earthquake may increase the surface temperature locally, which can be detected by satellite instruments. Therefore, despite cloud cover masks optical imagery, infrared data has already shown potential results for detecting temperature irregularities over active faults [9, 10].

Additionally, some atmospheric changes can also be observed before the occurrence of an earthquake [11], such as air ionization by radioactive radon emissions, thermal convection, anomalous linear cloud formation [12], or latent heat anomalies [13] are some of the means used as earthquake precursors.

The ionosphere is also a possible proxy of seismic activity. Anomalous TEC changes [14], and ionospheric scintillation (IS) perturbations related to earthquakes have been published in the last decades [15]. Different techniques to measure these parameters are summarized in this chapter. Traditional methods to measure these parameters include ionograms from ionosondes and GNSS satellite constellations to measure the effects produced in the signals when they are in view from a ground station. This chapter will review novel studies using GNSS RS techniques such as GNSS Reflectometry (GNSS-R) [16] and GNSS-Radio Occultations (GNSS-RO) to detect possible ionospheric precursors of earthquakes.

Other effects not easily categorized in previous environments have also shown the potential to serve as earthquake precursors. Those are, for example, geomagnetic perturbations [17], outgoing longwave radiation [18] and extremely low-frequency wave generation [19], or small perturbations in the Schumann resonance frequency [20].

1.2 Methodologies

Past studies have used mathematical models and statistical tools to study the correlation between these effects and the occurrence of earthquakes [21]. Typical statistical parameters used to detect the anomalies include the calculation of the standard deviation (STD), the interquartile range (IQT), and the Z-score. The influence of external phenomena on the studied parameters can hide or fake the correlation between the physical mechanism studied and the occurrence of the earthquake. Therefore, filtering techniques are very important to improve the data quality and reduce false alarms. Some previous observables can be affected by external factors, such as space weather, which is also largely influenced by solar activity and the interplanetary environment. In this sense, variables such as the planetary index (K_p) [22], the disturbance storm time index (Dst) [23], and the solar flux [24] should be considered to flag the affected observables.

1.2.1 Confusion matrix and receiver operating characteristic (ROC) curve

This section explains the methodologies used to compute the anomalies, highlighting explicitly the STD, the IQT, and the Z-Score techniques [17]. The Precursor deviation is determined from (Eq.(1)), wherein the Precursor mean corresponds to the average of the detrended Precursor over the entire period under consideration. The Dobrovolsky or strain radius (SR) estimates the size of the earthquake preparation area [25], and it is related to the Moment M_w through (Eq. (5)).

$${\text{Precursor}}_{\text{deviation}}=\text{Precursor}–{\text{Precursor}}_{\text{mean}}$$E1

$$\left|{\text{Precursor}}_{\text{deviation}}\right|\ge \text{C}·\text{STD}$$E2

$$\left|{\text{Precursor}}_{\text{deviation}}\right|\ge \text{C}·\text{IQT}$$E3

$$\text{ZScorePrecursorPrecursormean}/\text{STD}$$E4

$$S{R}_{\left[\text{km}\right]}={10}^{0.43Mw}$$E5

The test performance is assessed from six essential parameters used to compute the confusion matrix (CM): condition positive (CP), condition negative (CN), true positives (TP), false negatives (FN), false positives (FP), and true negatives (TN) [26]. The explanation for each of them is as follows:

• CP: count of samples with actual earthquakes.

• CN: count of samples without real earthquakes.

• TP: count of correct links where (Earthquake [+] and Anomaly [+]).

• FN: count of wrong links where (Earthquake [−] and Anomaly [−]).

• FP: count of wrong links where (Earthquake [−] and Anomaly [+]).

• TN: count of correct links where (Earthquake [+] and Anomaly [−]).

Additional metrics can be derived from them, such as the specificity, the precision, the recall, and the accuracy. Table 1 summarizes the different formulas used to compute these parameters, including the true positive rate (TPR), the false positive rate (FPR), the false negative rate (FNR), the true negative rate (TNR), the likelihood ratio (LR), the positive predictive value (PPV), the false omission rate (FOR), the negative predictive value, (NPV), the false discovery rate (FDR), the area under the curve (AUC), and the diagnostic odds ratio (DOR) [26].

CP	TP	FN	TPR = TP/P	FNR = FN/P
CN	FP	TN	FPR = FP/N	TNR = TN/N
Prevalence = CP/CP + CN	PPV = TP/PP	FOR = FN/PN	(LR+) = TPR/FPR	(LR−) = FNR/TNR
Accuracy = TP + TN/CP + CN	FDR = FP/PP	NPV = TN/PN	DOR = LR+/LR−
	F1 = 2·TP/(TP + FP + FN + TN)		MCC = TP·TN-FP·FN/((TP + FP) (TP + FN) (TN + FP) (TN + FN))

Table 1.

The details of the confusion matrix and formulas.

The receiver operating characteristic (ROC) curve is a graphical plot illustrating a binary classification method’s performance that depends on a threshold value. By sweeping this threshold, values of TPR against FPR are depicted on an x-y graph, leading to a curve which, ideally, goes from (0,0) to (1,1). The AUC parameter mentioned above measures the integrated area under this curve, which ranges from 0 to 1, being larger for better classifiers.

1.2.2 Machine learning tools for earthquake prediction

The confluence of all these elements constitutes what can be considered a Complex System with different sources of perturbations, interacting non-linearly and resulting in tiny observable changes [27]. These features make the problem suitable to be studied using machine-learning techniques. Neural networks can effectively extract the main features of a system like this, helping to understand and categorize the input sources involved in the problem and their relative impact.

Earthquake forecasting is commonly branded into two main types: model-based and precursor-based techniques. Model-based approaches involve hypothesizing a mathematical or a machine-learning model for reliable forecasting, incorporating an analysis of signals indicative of an approaching significant quake. Unfortunately, these signals are often undetectable before the earthquakes.

An innovative machine-learning algorithm was introduced in [28] for real-time earthquake precursor detection, leveraging the Support Vector Machine (SVM). This technique integrated daily GPS-TEC geomagnetic indices data and identified spatio-temporal anomalies linking them to earthquakes. A case study in Italy from January 1, 2014, to September 30, 2016, utilized various training, validation, and testing periods. During the validation period (266 days), the algorithm successfully detected precursors in 17 out of 21 earthquakes, and during the test period (282 days), it identified 22 out of 24 earthquakes, with 7 and 13 false alarms, respectively.

In [29], the authors explored physical and dynamic changes in seismic data collected from ten infrared and hyperspectral measurements from 2006 to 2013. They introduced the Inverse Boosting Pruning Trees (IBPT) machine-learning method, using satellite data related to 1371 earthquakes with M_w > 6 for short-term forecasting. Compared to other machine-learning approaches, IBPT outperformed six selected baselines, enhancing earthquake forecasting across diverse databases.

After this introduction, the rest of this book chapter is structured as follows. Section 2 reviews the studies on LST proxies for seismic activity. Sections 3 and 4 analyze the atmospheric and ionospheric perturbations due to earthquakes, respectively. Section 5 summarizes additional observables to study earthquake precursors. Section 6 presents the main conclusions.

2. Land surface temperature proxies

Detecting thermal anomalies before significant earthquakes is crucial to understand and predicting seismic activities. This importance arises from the ability to recognize phenomena related to thermal radiation during seismic preparation phases [30]. Satellite data provides an effective means to monitor seismic regions globally, nearly in real-time. Various data sources have been incorporated into this field in recent decades, advancing progressive anomaly detection techniques. This section reviews the advances and developments in pre-seismic thermal anomaly detection technology over the past decade.

In Ref. [31], the authors explore the association between ground surface warming and earthquakes by analyzing thermal data from Meteosat. The objective is to mitigate the impact of cloud coverage and daily weather fluctuations, enhancing the resolution of thermal maps. Detected thermal anomalies show a potential link between LST anomalies and future seismic events, notably observed in the L’Aquila earthquake on April 6, 2009.

In Ref. [31], the authors evaluate LST anomalies in various earthquake events of M_w > 6 in Pakistan from 2000 to 2020. They use MODIS imagery data over the epicenter region to examine thermal anomalies associated with earthquake prediction. The findings reveal pre-seismic LST variations ranging from 30 to 54°C, indicating a promising correlation between thermal anomalies and earthquakes.

In Ref. [32], the study analyzes AQUA/MODIS and TERRA/MODIS LST products before and after the M_w 7.9 earthquake in Nepal on April 25, 2015. It explores the thermal information related to the earthquake, revealing a steady increase of LST before the event, reaching a maximum during the quake, and returning to normal afterward. This temporal pattern could serve as premonitory information associated with the Nepal earthquake.

In Ref. [33], the study focuses on the trend of LST based on MODIS daily LST products from 2000 to 2017. The study applies an annual temperature cycle model to derive parameters such as the mean annual surface temperature (MAST), yearly amplitude of surface temperature, and phase shift. The study offers insights into the impact of earthquakes on the mountain thermal environment and its diverse responses to changes in vegetation cover and climatic conditions.

In Ref. [2], three earthquakes of M_w > 6 were studied using AQUA/MODIS and TERRA/MODIS LST data. Anomalies are detected using the IQT and mean ± 2. STD methods, revealing positive daytime anomalies based on mean ± 2σ. However, negative (or positive) anomalies are observed at night before earthquakes in Mexico and Bolivia, and occasionally, a negative anomaly is detected during the 10 days following an earthquake.

The study in Ref. [21] employs a novel approach to detect localized spatio-temporal fluctuations in hyper-temporal, geostationary-based LST data. The study covers ten areas worldwide, analyzing 20 large and small land-based earthquakes (M_w > 5.5 and depth < 35 km). The comparison of times and locations with and without earthquakes shows no distinct repeated patterns, indicating that earthquakes do not significantly impact detected LST anomalies.

In Ref. [34], a segmented threshold method was anticipated to detect anomalous microwave brightness LST linked with strong earthquakes in Sichuan province, China. The Index of Microwave Radiation Anomaly (IMRA) computed by this method was found to be improved nearly 2 months before the three strong earthquakes in 2008 Wenchuan (M = 7.8), 2013 Lushan (M = 6.6), and 2017 Jiuzhaigou (M = 6.5). The study acknowledges the limitation in quantitatively evaluating IMRA due to the infrequent occurrence of earthquakes.

In Ref. [35], the authors aimed to examine the anomalies in LST from GOES/ABI 16, 17 [26], and AQUA/MODIS, potentially associated with 1350 global earthquakes with a magnitude M_w ≥ 4 during 2020. In the same direction, the authors tried to detect and find a correlation between earthquakes and LST data in a long time series from Fengyun-2F/VISSRS [36], Himawari-8/AHI, and MSG/SEVIRI. In other studies, data acquired from numerous satellites [37] was ingested [38] and processed [39] to identify unusual anomalies before earthquakes. The precursor value tends to be high near the earthquake’s epicenter 1 to 7 days before the seismic event [8]. Traditionally, this outcome has been attributed to the thermal flux emitted from the Earth’s crust in regions prone to seismic activity. To improve the data quality, filtering is applied exclusively to nighttime datasets between 00:00 and 06:00 am, thereby mitigating the influence of sunlight and radiation on precursor values. Subsequently, the data undergo fusion (averaging, standard deviation, minimum, and maximum values) for each location and day. The organized data is then chronologically sorted to form a time series. The subsequent phase involves computing the precursor average using the Strain Radius (SR) formula, delineated in Eq. (5). The final stage encompasses calculating anomalies based on prescribed formulas. Once the anomalies are computed, they are correlated and associated with potential earthquakes. A confusion matrix is generated to optimize detection parameters. Finally, the correlation is visually represented through charts and maps.

Most existing works apply RS data from geostationary and polar satellites equipped with thermal or microwave sensors covering temporal windows spanning a few months to years. Similarities in data processing methods, such as LST data ingestion filtering, time series analysis, and anomaly detection, are obvious. However, differences arise in the data processing and LST anomaly detection. Recent contributions stand out in several critical aspects: (1) satellite constellations by using multiple satellite constellations from geostationary satellites, ensuring comprehensive global LST coverage; (2) data filtering using diverse LST data filters, including quality flags, water bodies, and fire masks; (3) extensive time series from 2010 to 2021, facilitating correlation analyses across different seasons and temporal conditions; (4) a large number of land earthquakes, enabling analyses across varied land covers, altitudes, depths, and altitude regions for realistic statistical results; (5) statistical approaches not used before such as the Confusion Matrix and the ROC curves; (6) various thresholds analytically tested to determine the optimal configuration for correlation analysis; and finally (7) the classified correlation was conducted based on diverse parameters, including M_w, depth, land cover, regions.

Figure 1c1–c6 illustrates the processed AQUA/MODIS LST data near Ankara, Turkey, from January 20 to 25, 2020. The cyan star indicates the earthquake’s epicenter. A marginal rise in LST was observed 3 days before the earthquake. However, a notable positive anomaly in LST was documented on the day of the earthquake, followed by a decrease in LST the day after.

Figure 1.

Sample LST maps, original, detrended, and anomalies time series (°C) measured with ABI/GOES or MODIS/AQUA for epicenters ID977 in 2020 [35].

2.1 Confusion matrix and ROC curve

Figure 2 illustrates the outcome of the CM: the TPR and the FPR decrease as the “C” coefficient increases. Consequently, a small “C”, especially at C = 0.7, yields the optimal threshold to detect LST anomalies. However, this comes at the cost of a relatively high rate of false alarms. Conversely, larger “C” coefficients identify fewer LST anomalies. The higher the DOR, the better the test performance. In our study, a DOR of approximately 13.5 corresponds to the optimal values for C = 0.7. The shortest distance (d) from the receiver operating characteristic (ROC) curve to the point (TPR, FPR) (0,1) is approximately 0.49. Therefore, the optimal “C” for earthquake monitoring is 0.7, with an SR of 50 km.

Figure 2.

Optimum ROC curve for C = 0.7 of correlation between LST and earthquake occurrence.

2.2 Impact of Mw, depth, altitude, and land cover on the confusion matrix

Figure 3a shows a consistent TPR, AUC, and DOR increase as M_w increases. The LST anomalies detection is high where the M_w is large. TNR and accuracy remain relatively constant, around 0.8 across different M_w levels. Conversely, FNR and “D” escalate as M_w decreases. The peak TPR, AUC, DOR, minimum FPR, and “D” are observed for deep earthquakes. TNR, FPR, and ACC exhibit uniformity across depths.

Figure 3.

TPR and FPR separated by magnitude, depth, altitude, and land cover.

From Figure 3b, epicenters within 50 to 100 km and over 400 km strongly associate with prior LST anomalies, displaying high confidence and less noise, reaching the highest values of TPR, up to ~0.7.

At low altitudes (1001–2000 m), the highest TPR, AUC, and DOR, coupled with the lowest FPR and “D,” are achieved, highlighting numerous correlated LST anomalies. TNR, FPR, and ACC remain consistent across different altitudes, as shown in Figure 3c. Consequently, the land cover type and altitude choice are crucial in accurately detecting LST anomalies, as they significantly impact the measurement using RS techniques.

From Figure 3d, the maximum TPR, AUC, and DOR are noted for different land covers, while the lowest FPR and “D” occur in mosaic and snow ice covers, respectively. Savannahs, grassland, and cropland record an average detection, whereas forests and shrubs show the worst performance, suggesting that vegetation covers mask LST anomalies. TNR, FPR, and ACC demonstrate similarity across various land covers, emphasizing the significance of land cover type in LST anomaly detection using RS techniques.

Despite various promising approaches to identifying pre-seismic thermal anomalies, the LST anomaly detection method still needs to be more complete and requires higher precision. The detection of pre-seismic thermal anomalies may also be influenced by factors such as topography, soil moisture, vegetation, and other atmospheric conditions.

3. Atmospheric earthquake precursors

The atmosphere is the gaseous layer above Earth’s surface. It contains different layers characterized by different composition, density, and temperature. Up to 85% of the total mass of the atmosphere is confined in the lower layer, the troposphere. The overall composition of the Earth’s atmosphere is very stable, containing, on average, around 78.08% nitrogen, 20.95% oxygen, 0.93% argon, 0.04% carbon dioxide, and traces of hydrogen, helium, and other noble gases. A variable amount of water vapor is also present, about 1% on average at sea level. Several signatures of seismic activity have been reported within the atmosphere, which can be due to the proximity of these two layers and the energy interchanges between them. Among the affected physical parameters studied, the literature reports on radon emissions in the vicinity of seismic regions, the generation of gravity waves in the atmosphere, the formation of anomalous clouds in the troposphere, the so-called “earthquake clouds,” and anomalies in the surface latent heat flux [6].

Radon emissions are one of the most reported seismic precursors within the atmosphere. Radon (Rn) is a noble gas in the sixth period, in which all its isotopes are radioactive. Its most important one (Rn-222) belongs to the radium and uranium-238 decay chains, with a half-life of 3.8 days. Being a radioactive element, it can be easily tracked by the radiation it generates and its sub-products. Additionally, the short half-life limits its detection to short periods after its generation in the atmosphere. It can come from subsurface trapped gas after new faults and cracks appear in preparation for an earthquake.

A very extensive review of the methods using radon emissions, especially as earthquake precursors, is made in [40], where many past studies and techniques are discussed.

An especially large magnitude earthquake was reported in [41], where the sudden increase of the soil-gas Rn-222 was detected 22 days before a 7.1-magnitude earthquake in November 1994 in the Philippines. The peak reached a value of up to six times the background concentration of Rn. In the study, it was discarded that the close pass of typhoon Teresa was the possible origin of the peak because 1 year later, the larger typhoon Angela did not produce the same disturbance in the Rn concentration.

Surface latent heat flux (SLHF) and air and surface temperature anomalies have been reported in [42] for two earthquakes located in a seismic area of Iran in 2010 and 2011, which reached magnitudes 6.5 and 6, respectively. The SLHF is the thermal energy flux transmitted from the Earth’s surface to the atmosphere. Therefore, it is highly related to the surface temperature, which may increase during the earthquake preparation.

A study of the SLHF anomalies found before seven earthquakes in India, Taiwan, and Mexico from 1993 to 2003, with magnitudes ranging from 5.4 to 7.8, was presented in [43]. The SLHF anomalies were recorded within 2 to 7 days before the earthquakes, and they were attributed to the increase of infrared thermal temperature in the earthquake preparation region.

In Ref. [44], air temperature anomalies were also found before three earthquakes in Turkey, Canada, and Argentina, with magnitudes of 6.1, 6.2, and 6.4, respectively. This work studied the average and STD values of the atmospheric daily temperature for 5 years before the earthquake date, and the positive or negative anomalies were recorded for the earthquake year. As a result, mostly positive anomalies were found 1 month before the earthquake and some a few days before. Additionally, some anomalies were found after the earthquake.

Another possible atmospheric earthquake precursor reported in the literature is the anomaly observed in the formation and evolution of tropospheric clouds near a fault before the earthquake. In Ref. [15], three case studies are reported using geostationary satellite images in Italy and Iran. The authors claim that the linear stationary shape of clouds, which did not move for about 8 hours, is a clear signature of an earthquake precursor, estimating the magnitude and even the date. According to the study, they only use the cloud shape and size to predict, and the prediction could be improved by using other data, such as radon or groundwater data. Nevertheless, the authors do not describe the methodology used to derive the earthquake occurrence with the shape or size of the anomalies. Additionally, the formation and shape of clouds can be linked to the terrain orography, which needs to be discussed in the report. Two years later, the study was questioned in Ref. [45].

In addition to previous atmospheric signatures, in Ref. [46], atmospheric ozone disturbances were observed before two large earthquakes occurred in Indonesia (M = 7.3) and Fiji (M = 8.2) in 2019 and 2018, respectively. The study analyzed the ionospheric TEC parameter and the ozone data from AIRS, OMI, and TOMS-like satellites.

4. Ionospheric earthquake precursors

Along with the previously analyzed earthquake precursors, the ionosphere has been reported in many studies to be affected by seismic activity before and after earthquakes. This section discusses recent advances in the precursor signals observed in the ionosphere, particularly in the electron density (Ne) and TEC, as well as the IS phenomena.

4.1 Electron density and TEC anomalies

In Ref. [47], the daily variations of TEC in the ionosphere and the daily energy released near three studied earthquakes were examined and correlated, confirming a minor correlation. It was deduced that more parameters are essential to understand the cause of ionospheric TEC anomalies [48].

In Ref. [49], a study is presented using plasma density, electron temperature, and slant TEC from the ESA Swarm constellation. Two Mw ≥ 6.5 earthquakes were analyzed in this study. The first was a 6.7 magnitude earthquake in Turkey in March, and the second was in the USA in January 2022, with a magnitude of 6.5.

In Ref. [50], positive and negative VTEC anomalies were found before the December 2016 earthquake in Chile. Disturbances occurred 8 to 6 days before the earthquake within the Dobrovosky preparation region. The solar and geomagnetic conditions were also studied, showing that they remained quiet and stable.

In Ref. [51], the Ne was analyzed from the Swarm mission 1 month before and after 12 strong earthquakes of Mw > 6, between 2014 and 2017. The study used the Dobrovosky SR around epicenters to limit the region to correlate with Ne anomalies.

In Ref. [52], an innovative methodology based on TEC measurements was presented to estimate the moment of Co-seismic Ionospheric Disturbances (CID). The study focused on three large earthquakes: the 2015 Illapel, the 2014 Iquique, and the 2011 Sanriku-Oki. Results indicated that CIDs arrived 250–430 seconds after the seismic wave peak or 350–700 seconds after the earthquake onset time.

4.2 Ionospheric scintillation measured from GNSS signals

Ionospheric scintillations (IS) are the rapid fluctuations of the phase and intensity of an electromagnetic wave traversing the ionosphere. It is mainly due to the spatio-temporal changes in the local ionospheric Ne, and it is primarily driven by the solar influence, usually showing a maximum during the post-sunset hours.

Only a small fraction of the studies relating ionospheric perturbations to earthquakes use the IS as a potential proxy. These studies usually take GPS data from ground stations [53] or ground-based ionosondes [54] to measure the S4 intensity scintillation index, and correlate it to regional earthquakes.

This section presents IS signatures as possible earthquake precursors reported in the literature. They are studied using data recorded in ground stations and novel techniques to derive the IS from remotely sensed data from satellites, particularly from GNSS Reflectometry and GNSS radio-occultation.

4.2.1 GNSS ground stations

GNSS ground stations are infrastructures installed worldwide to control and monitor the signals received from GNSS constellations, offering important TEC and IS data. This information is used to make ionospheric corrections for navigation and to evaluate the quality of service. Over the past decades, these data have been used to monitor the state of the ionosphere. However, coverage is limited. GNSS stations are usually strategically placed at fixed ground locations, resulting in sparse installations that only offer data for their immediate local regions. Therefore, ground station data is a good approach when studying a local event or a high-seismicity region.

In Ref. [55], GPS and low-frequency measurements from multiple stations were analyzed to identify ionospheric anomalies related to the South West Banten (Indonesia) earthquake with an Mw of 7.4. The analysis revealed increased IS 5 days before and after the earthquake, with the anomaly being more extensive after the main seismic shock. In [56], the authors examine the impact of the seismic activity associated with the 2021 La Palma (Canary Islands, Spain) volcanic eruption using IS anomalies detected from ground stations near the volcano. The ground station database was pre-processed using the algorithms described in [17]. The S4 data here are obtained using the 1 min averaged signal’s intensity received in ground stations from several constellations of GNSS satellites. Only GNSS satellites above 30° of elevation were used to avoid wrong high scintillation values. The results showed a small but detectable correlation between the IS and the earthquake energy released per unit of time. The best correlation between both variables was found to occur 18 h, and 7–8 days after earthquakes. The work also studied IS measured by GNSS-R and GNSS-RO satellites, in this last case where larger correlations were found before the earthquakes. A cross-correlation between the S4 measurements and the Kp and solar flux (F10.7) was also performed, showing smaller and not overlapping in time with the correlations found between earthquakes and S4.

4.2.2 Satellite GNSS-R data

In Ref. [57], the authors demonstrated the efficiency of high-rate GNSS data in detecting low-magnitude anthropogenic earthquakes. This study showed the potential impact of the infrastructure of GNSS ground stations and was successfully integrated into moment tensor calculations, focusing on two earthquakes with magnitudes greater than 3.5.

GNSS-R works as a multi-static radar that uses opportunistic navigation signals. The receiver calculates signal delay and Doppler frequency, enabling simultaneous tracking of multiple transmitters, and their specular reflection points [58]. Applications of GNSS-R include ocean altimetry, surface wind speed, cryosphere studies, and surface soil moisture. NASA’s CYGNSS represents the first operational mission based on GNSS-R, featuring eight micro-satellites equipped with SGR-ReSi receiver payloads from SSTL. The SGR-ReSi produces Integrated Delay Doppler Maps (DDMI) by capturing reflected GPS signal power, which depends on delay and Doppler shift [59]. In this study, CYGNSS data at level 1, version 2.0, has been used. The study explores GNSS signals traversing the ionosphere in their down-welling and up-welling paths, providing insights into IS at GNSS frequencies.

In previous contributions [18, 60, 61], the IS data derived from CYGNSS GNSS-R data were scrutinized and correlated to ocean-based earthquakes from 2017 to 2021. Again, a small, but detectable correlation was identified between S4 anomalies and earthquake occurrences a few days before.

After acquiring CYGNSS data, a nighttime filtering process is employed between 00:00 and 06:00 am local time to eliminate diurnal effects and post-sunset scintillation. While IS may persist beyond midnight on certain occasions, this local time window is chosen to ensure a substantial dataset for processing. The data is then aggregated to determine the mean, STD as shown in Figure 4a, IQT, and total count of pixels with identical geolocation (longitude and latitude) and local time. The spatial resolution in the figures is set to 0.5°.

Figure 4.

(a) STD of S₄ with 4 m/s sea wind speed filters between 2017 and 2021, (b) S₄ anomalies classification (TP, FP, and FN), (c) S₄ anomalies associated with the Solomon Islands earthquake, and (d) S₄ anomalies associated with the Maldives earthquake.

Subsequently, the dataset is converted into a raster format, with a filter applied to exclude pixels with high sea wind speed. This filter is important as a coherent reflection for observing IS in GNSS-R data, which can only occur when wind speeds are below 2–3 m/s. An elevated solar activity flag eliminates potential sources of S4 anomalies unrelated to earthquakes. It is worth noting that the STD, IQR, mean, and median are then computed over 7 years to derive daily S4 anomalies, as illustrated in Figure 4b.

In areas near earthquake epicenters, there is a general increase in IS, particularly when the Sea Wind Speed Filter (SWSF) is less than 2 m/s, and the Fixed Search Radius (FSR) is under 50 km. S4 anomalies, ranging from 0.08 to 0.23, can be observed up to 7 days before earthquakes. Despite a visual and numerical correlation, not all earthquakes connect to S4 anomalies, and vice-versa. The forecast of earthquakes remains complex due to the sporadic occurrence of minor S4 anomalies.

For example, in the Solomon Islands earthquake (Figure 4c) in February 2018, with an Mw of 4.4 at a depth of 10 km, a positive S4 anomaly of 0.11 was noticed 3 days before the event. Minor positive S4 anomalies, approximately 0.1, were detected 1 and 3 days before the Maldives earthquake (Figure 4d) occurred in April 2019.

Figure 5 illustrates that the ROC curves derived from CYGNSS mission data using the STD show remarkable similarity and near-identical patterns. As the parameter “C” increases, the TPR and FPR exhibit distinct trends: TPR rises while FPR decreases. Consequently, detecting S4 anomalies improves for smaller “C” coefficients. However, it is important to note that the false alarm rate is also notably elevated.

Figure 5.

Confusion matrix for S₄ anomalies using STD and IQT methods with 2 m/s sea wind speed filters.

For instance, when “C” is set to 0.1, TPR equals 1, and FPR equals 1. Conversely, at “C” = 2.5, TPR is 0, and FPR is 0. The optimal configuration is observed at “C” = 0.5, and a false alarm rate FSR of 50 km. The result is the highest area under the curve (AUC) and a straight distance (d) of 0.62. Therefore, the recommended threshold is “C” = 0.5, paired with a sliding window size factor SWSF of 2 m/s.

4.2.3 Satellite GNSS-RO data

GNSS-RO stands for GNSS Radio-Occultation, and it constitutes a technique employed for atmospheric and ionospheric sounding by using the fluctuations of the received signal from GNSS satellites, but in this case, in a radio-occultation geometry. In this configuration, rays do not reflect on the Earth’s surface but cross the ionosphere tangentially when the GNSS satellites fall over or rise above the horizon from another satellite, the GNSS-RO receiver. This technique has the advantage of not relying on surface properties, which may affect the properties of the wave, as was the case in GNSS-R.

In Ref. [62], a novel method was presented for automatically categorizing disturbances in the lower ionosphere. This approach utilized GNSS-RO data from Spire’s CubeSats constellation, integrating signal processing methods and semi-supervised machine learning, employing spectral clustering within a metric space of wavelet spectra. This innovative approach established an automated system for global monitoring of ionospheric disturbances.

As shown in Ref. [63], a similar analysis was conducted, but using GNSS-RO data to examine the S₄ intensity scintillation index throughout 2022 in the Coral Sea region, encompassing Papua New Guinea, Solomon Islands, Vanuatu, and New Caledonia. The selection of this area is based on the substantial occurrence of earthquakes, aimed at enhancing the reliability of the statistical analysis. Within this timeframe and geographical scope, 201 earthquakes with magnitudes M_w ≥ 5 and 22 earthquakes with M_w ≥ 6 have been documented.

The IS data in this study was sourced from the GNSS-RO COSMIC-2 constellation [64]. Applying GNSS-RO for IS also enables the investigation of terrestrial and oceanic regions. Earthquake data was obtained from the United States Geological Survey (USGS) database, encompassing earthquake magnitude, date, location, and hypocenter depth [61].

COSMIC-2 provides GPS and LEO coordinates for each occultation but lacks a reference coordinate for determining the S₄ location. The tangent point is calculated using the GPS and LEO satellite positions from COSMIC-2 using the same methods described in [65], which can be seen in [66]. This is later used to filter data based on the distance around earthquake epicenters or latitude-longitude pairs.

Before data processing, various filters are applied to ensure accurate S₄ estimation. Initially, it is mandated that the pierce points in the profile data of the occultation path account for a Slant TEC value ≥80% of the maximum STEC value in that occultation profile as described in [56] and [63]. This filter enhances the likelihood of isolating points within the ionosphere susceptible to scintillation. Additionally, the elevation must be equal to or lower than 0° to confirm the occurrence of an occultation. Another filter involves a nighttime criterion, as explained in previous sections.

To compare the S₄ scintillation index with seismic activity, the mean and the 95th percentile of S₄ values were computed in 1-day intervals throughout the entire study period. The spatial resolution is set at 1° of latitude and longitude. Throughout the statistical analysis, every latitude-longitude pair within the study area undergoes systematic examination. For each pair, S₄ values within a search radius of 100, 200, and 300 km from the pair’s center are gathered, and the mean and 95th percentile of these values were calculated. Epicenter’s coordinates are retrieved from the USGS database. To reconcile the spatial resolution disparity between S₄ and earthquake data, the latitudes and longitudes of earthquakes are approximated to their nearest latitude-longitude pair.

This study spans 2022, repeating the analysis for all latitude-longitude pairs within the defined study area. Consequently, there are instances where no earthquakes occur on a specific day and latitude-longitude pair. This aspect of the research aided in evaluating non-detection reliability while simultaneously investigating IS anomalies around earthquake epicenters.

Next, the results of two case studies in this area are shown. These two case studies presented occur 8 days apart, with epicenters positioned approximately 250 km from each other. Consequently, there is a possibility of mutual influence between the earthquakes. For each case study, a temporal analysis and spatial representation are conducted. In the temporal analysis, we examine the daily changes in S₄ values near the epicenter within a temporal window of ±15 days surrounding the event. Concerning the spatial representation, the distribution of S₄ values is illustrated geographically in pixel maps computed daily near the earthquake. The first earthquake under consideration occurred on September 2nd, 2022, with a magnitude of M_w = 6.1 and ID “us7000i4st”.

Figure 6 illustrates the daily averages of the S₄ scintillation index within search radii of 100, 200, and 300 km from the epicenter. Notably, ionospheric intensity scintillation anomalies are discernible several days before the earthquake, particularly between 5 and 6 days preceding the event. To quantify this deviation, the median value of the S₄ averages is computed within a 200 km radius, resulting in 0.079. In contrast, values 5 and 6 days before the event are 0.30 and 0.26, respectively—approximately 3 to 4 times larger than the median of the averages within the ±15 days surrounding the event. The same trend is observed if the same analysis using the 95th percentile is represented, as seen in [63].

Figure 6.

Average values of the S₄ indexes within circular areas of 100, 200, and 300 km from the earthquake’s epicenter “us7000i4st”.

In addition to the temporal analysis, a spatial representation is examined. This analysis involves 0.2° × 0.2° pixels in latitude-longitude for both case studies, each assigned a distinct color based on S₄ average values. Additionally, a gray circle with a radius of 200 km is delineated, centered on the epicenter of the analyzed earthquake, symbolized by a black star. The spatial representation examines the proximity of high S₄ values to the earthquake’s epicenter, providing a visual portrayal of the scene.

For this case study, the spatial representation is presented in Figure 7. Various days are depicted, including those with peak S₄ values (Figure 7a and b), as well as the day of the earthquake occurrence (Figure 7c). These figures reveal that the days with peak S₄ values are 5 and 6 days before the earthquake, specifically on the 28 and 27 August, respectively. These days, certain pixels inside the drawn radii exhibit high values of S₄. Conversely, the maximum values are lower on the event day, September 2.

Figure 7.

Average S₄ maps for selected days centered on the earthquake “us7000i4st”.

The second case study that is shown involves an earthquake on September 10, 2022, with a magnitude of M_w = 7.6 and ID “us6000iitd”. Figure 8 is created using the same method described in Figure 6. In this instance, two notable peaks are observed again, occurring 13 and 14 days and 3 and 2 days before the seismic event. The first set of days, corresponding to 5 and 6 days from the “us7000i4st” earthquake (earthquake M6.1 in Figure 7), is due to the proximity between the two events. Consequently, the ionosphere is affected by anomalies 3 and 2 days before the earthquake. The S₄ representative values within a search radius of 200 km for the 3 and 2 days before the earthquake are 0.37 and 0.30, respectively. Comparing these values to the dataset’s median, which is 0.15, the anomalies are approximately 1.97 to 2.43 times higher than the median.

Figure 8.

Average values of the S₄ indexes within circular areas of 100, 200, and 300 km from the earthquake’s epicenter “us6000iitd”.

The spatial representation maps in Figure 9 correspond to the peaks occurring 2 and 3 days before the earthquake. Figure 9a and b correspond to the 8th and 7th of October, respectively, and the day of the earthquake itself, Figure 9c. It is evident, again, that only in Figure 9a and b there are values of S₄ around or higher than 0.4. In contrast, there are only values below 0.1 on the day of the earthquake.

Figure 9.

Average S₄ maps for selected days centered on the earthquake “us6000iitd”.

An analysis of all earthquakes within the studied region and period is conducted to assess the presence of a positive correlation and to determine the optimal detection method through confusion matrices. As detailed in Table 2, various parameters are tested to generate diverse confusion matrices and ROC curves. The outcomes of the confusion matrices, along with their associated metrics and ROC curves, contribute to finding the optimal approach.

Parameter swept	Values swept
The radius of the search area	100, 200, and 300 km
Type of S4 representative value	Computed with mean or 95th percentile
Minimum magnitude of earthquakes studied	Mw = 5, 5.5, 6
Depth of the earthquake	[0, 100) km, [100, 200) km, [200, ∞) km
Magnetic field activity	No filter or KP ≤ 5
S4 Anomaly occurring	[−6, −3], [−3, 0], [0,+3] days before/after the earthquake
Detection threshold values	1, 1.25, 1.5, 1.75, 2 … 4

Table 2.

Description of optimization parameters to generate different confusion matrices.

The investigation of these parameters aims to identify the range from which IS effects can be discerned. This serves as a filtering mechanism, determining the effectiveness of the detection method for earthquakes with specific characteristics. For example, a magnetic field activity filter is implemented to mitigate anomalies from solar or geomagnetic storms. The K_P is utilized to ascertain the presence of a solar or geomagnetic storm, describing the intensity of magnetic field activity.

Ratios are established to determine the likelihood of an S₄ anomaly occurrence numerically. Two types of ratios are defined, leading to two additional iterations to select the optimal detection method. The first ratio (r_{1, x}) involves the ratio between a specific day’s S₄ mean value and the median of the S₄ representative value within ±15 days from the computed day. The second ratio (r_{2, x}) is an IQR-based threshold, calculated as the difference between the median of the S₄ representative values and the daily value of the S₄ representative value, divided by IQR/2. These are also described mathematically in [63], and are computed for all days of the year, not solely around the days of earthquake occurrences. Once both ratios are calculated, they are tested with threshold values ranging from 1 to 4 in increments of 0.25. An anomaly is classified as true or false based on whether r_{1, x,} or r_{2, x} exceeds the tested threshold.

After the optimization involving all the parameters, it is shown that the most favorable results are obtained with the window of occurrence day [−6, −3] and the 200 km radius, the S₄ representative value calculated using the mean and the r_{1, x} type of ratio as the method of detection. As for other filters, as shown in Figure 10, the best results are obtained when all the filters proposed in Table 2 are used altogether, obtaining the highest values of AUC.

Figure 10.

ROC curve comparison for S₄ average within the searching area of 200 km, ratio r_1,x, and anomaly occurrence [−6, −3] days with different filters.

The final step in this study was to determine the optimal threshold. This involves examining the ROC curve in Figure 10 and considering metrics such as F₁, MCC, and DOR. After considering those metrics, the best results were obtained using a threshold of 3.75. Out of 14 positive cases, 5 were correctly detected, yielding a TPR of 35.7%. Additionally, an FPR of 7.1% is observed, with 21,027 false positive cases among 295,636 negative cases. The ROC curve suggests a positive correlation, given that the TPR value is approximately five times higher than the FPR, and the AUC is 0.69.

5. Other earthquake precursors

In addition to the previously reported possible earthquake precursors in the LST, the atmosphere, and the ionosphere, this section reviews other related phenomena that are difficult to characterize in any of the previous sections. In this section, geomagnetic field perturbations, electromagnetic emissions, and disturbances in the Schumann resonance as possible earthquake precursors are revised.

5.1 Geomagnetic field perturbations

Geomagnetic disturbances have also been detected in association with earthquakes. Some studies have related post-seismic geomagnetic disturbances. For example, a geomagnetic pulsation lasting for about 3.6 minutes was observed in Thailand 12 minutes after the Sumatra 9.3 magnitude earthquake on December 26, 2004, as reported in [66].

Moreover, in [47], efforts were made to detect perturbations in the geomagnetic field in the ionosphere potentially caused by earthquakes in the lithosphere. A magnetic field perturbation was observed on the orbit track crossing the studied epicenters. The study found 12 satellite tracks with potential co-seismic disturbances related to 10 earthquakes ranging from M_w 5.6 to 6.9, confirming the probability that seismic events produced the detected disturbances.

In Ref. [67], magnetic field data from the Swarm Alpha satellite were analyzed using non-negative matrix factorization (NMF) before the earthquake in Ecuador that occurred in 2016, with M_w = 7.8. The study included data collected under quiet and disturbed geomagnetic conditions, minimizing false alarms. Observations suggested a potential association described by one of the Lithosphere–Atmosphere–Ionosphere-Coupling (LAIC) models.

In Ref. [68], authors deployed a magnetic monitoring network in earthquake-prone areas of Sichuan, China, conducting long-term observations and real-time analysis of geomagnetic data. Results showed a high correlation between pre-earthquake geomagnetic anomalies and earthquakes. A new algorithm using geomagnetic anomaly energy was introduced to predict the time of an earthquake’s occurrence.

In Ref. [69], the authors inspected magnetic field variations from earthquake-associated satellites, utilizing the China Seismo-Electromagnetic Satellite (CSES) and ESA’s three spacecraft Swarm fleet. The study accessed links between specific lithospheric or near-surface sources and ionospheric magnetic field measurements, revealing increased fluctuations in higher-frequency components.

Swarm satellites represent an interesting way to study possible earthquake precursors in the geomagnetic field, as reported in several studies in the last few years.

In Ref. [70], the 2015 Nepal earthquake, with a magnitude of 7.8, is under study. Time series of the magnetic field vector and intensity were analyzed. Only data from night orbits were used to avoid noisy magnetic periods. The time series were derived and detrended to show fast fluctuations. The number of anomalies within the Dobrowolski circular area around the earthquake epicenter was recorded, and it was revealed that they increased in the days before the earthquake. The perturbations are observed mainly in the Y component (East), smaller in the magnetic field intensity and other elements.

5.2 Electromagnetic radiation

Some studies, including some of the previously revised in the chapter, report that earthquakes can generate low to very low-frequency electromagnetic waves, which can travel from the preparation earthquake region. These waves fall into the extremely low frequency (ELF) or very low frequency (VLF) ranges, 3–30 Hz and 3–30 kHz, respectively. Some satellite missions have been particularly designed to try to detect and correlate this electromagnetic emission with earthquakes. The first one was the French CNES Demeter mission, launched in 2004. It could detect changes in the magnetic field and electromagnetic waves; it also had an ion analyzer, a Langmuir probe, and an energetic particle detector.

Another satellite investigating electromagnetic radiation related to earthquake occurrence is the one already mentioned. China Seismo-Electromagnetic Satellite (CSES). This satellite orbits in a 507 km Sun-synchronous orbit with a descending node at 14:00 LT. Its eight instruments include magnetometers, electric field detectors, a GNSS-RO instrument, a plasma analyzer, a Langmuir probe, and high-energy particle detectors. The mission is intended to study the correlation between ionospheric-electromagnetic signatures and seismic activity.

In a study from 1992 [71], seismo-electromagnetic waves were observed by the COSMOS-1809 satellite orbiting above seismic regions in Armenia. The earthquake under study occurred in Spitak, Armenia, in December 1988, with a magnitude of 6.7. Intense EM waves were detected in frequencies under 450 Hz at the L-shells of the earthquake during 12 out of the 13 orbits that passed within 6° in longitude from the epicenter.

In Ref. [72], the VLF frequency band was analyzed as a possible source of earthquake precursors in the shape of seismic radio wave signals. They found pulse-like signals traveling long distances up to 500 km from the epicenter of preparing earthquakes, indicating waveguide mode propagation between the ionosphere and the ground. The great 1994 Kuril Islands Earthquake of magnitude 8.1 occurred East of Hokkaido, and the number of pulses detected in the 5 hours before the earthquake increased about nine times above the background value, with its maximum around 20 minutes before the earthquake.

In Ref. [73], seismic anomalies before two strong earthquakes with M_w > 7 in western China were detected using ground-based observation of VLF signals from monitoring stations. A full wave model was used to research possible factors inducing seismic anomalies on the received VLF electric field. Results demonstrated that anomalies could be caused by ascending/descending the ionosphere’s bottom height.

In Ref. [74], the research inspected the correlation between earthquakes with M_w ≥ 6 globally and fluctuations in ionospheric Ne between August 2018 and March 2023. Using data from the Langmuir probe onboard the CSES during nighttime, statistical analysis revealed that larger earthquakes were associated with more anomalous phenomena. Anomalies appeared at specific intervals, including co-seismic effects and possible precursor candidates.

In Ref. [75], Ne disturbances were detected around the Maerkang earthquake (M_w = 6) on June 9th, 2022, using multi-data from the global ionospheric map, ionospheric TEC inverted from GPS observations, and the critical frequency of the F2 layer from the ionosonde. Seismo-ionospheric disturbances were observed 7 to 5 days before the earthquake, serving as significant signals of the upcoming main shock.

5.3 Schumann resonance frequency

The Schumann resonance occurs due to the resonant cavity formed between the Earth’s surface and the ionosphere, acting as an electric conductor. The global extension of the cavity adjusts some resonance frequencies that fall in the ELF band, with its fundamental frequency at 7.83 Hz, and higher modes at 14.1 Hz and 20.3 Hz. The excitation of the cavity comes from natural sources such as electric discharges from lightning.

A possible link between earthquakes and the Schumann resonance has been recently theorized, and several studies have been conducted on the detection and correlation of these phenomena. A possible physical explanation of these anomalies is that ELF radiation can be generated during the preparation of earthquakes, which can rise from underground and interact with the Schumann resonance [76].

In Ref. [55], two anomalies on the Schumann resonance are discussed before large earthquake occurrences in Taiwan and Japan. The first phenomenon is the enhancement of the Schumann resonance fourth harmonic in Japan, possibly associated with the September 1999 Chi-chi earthquake in Taiwan. This earthquake had an magnitude of 7.6, and the hypocenter was 30 km deep. The second phenomenon was associated with two huge earthquakes in Japan. The resonance frequency was shifted by about 2 Hz from the typical Schumann resonance frequency: ~16 Hz and 18 Hz in two different modes. In this case, the earthquakes producing these anomalies were the 2004 Niigata-chutes earthquake, and the 2007 Noto-Hanto earthquake, with magnitudes of 6.8 and 6.9, respectively. The anomalies started around 20 h before the earthquakes, were maximum around 2 h before, and decreased after the earthquake.

In Ref. [77], a theoretical study is performed on how earthquakes affect the Schumann resonance spectra. The authors utilized a numerical model of the ELF radio wave scattering generated by a localized non-uniformity of the air conductivity profile likely to appear before an earthquake. The results showed an overall elevation of the Schumann resonance spectra, with some dependence from the distance to the earthquake epicenter.

6. Conclusions

This chapter has summarized the potential application of remotely sensed data as proxies for earthquake prediction. Numerous subtle indicators linked with earthquakes have been detected, such as land displacements, accelerations, gas emissions, temperature anomalies, and ionospheric effects. RS data offers benefits over in-situ data collection, such as extensive spatial coverage in remote regions. This chapter highlights investigations into various Earth observation observables, including anomalies in surface temperatures, atmospheric and ionosphere precursors, and additional physical mechanisms like fluctuations in the geomagnetic field. Despite the presence of seismic hazard risk models and the discoveries from these investigations, accurate earthquake forecast remains a significant challenge.

Acknowledgments

This work was supported in part by the grant GENESIS PID2021-126436OB-C21 from the Programa Estatal para Impulsar la Investigación Científico- Técnica y su Transferencia, del Plan Estatal de Investigación Científica, Técnica y de Innovación 2021-2023 (Spain), and in part by the European Social Fund (ESF), Grant RYC-2016-2018 financed by MCIN/AEI /10.13039/501100011033. Badr-Eddine Boudriki Semlali received support from an FI grant: 2021 FI_B 00471 from FI AGAUR 2021.

Data availability statement

The GOES/ABI, MSG/SEVIRI, and Himawari-8/AHI data used in this study comes from the Copernicus global land service https://land.copernicus.eu/global/products/lst and earthquake data comes from the open USGS earthquake database available at https://earthquake.usgs.gov/earthquakes/search/.

CYGNSS GNSS-R data is freely available in the NASA PODAAC database. https://doi.org/10.5067/CYGNS-L1X30.

GNSS-RO data from the COSMIC-2 satellite is freely available in the UCAR/NCAR CDAAC database. https://doi.org/10.5065/t353-c093

Other data are the property of a third-party organization, such as the Spire GNSS-RO data.