First Results on the Systematic Search of Land Surface Temperature Anomalies as Earthquakes Precursors

Abstract

Every year, earthquakes cause thousands of casualties and high economic losses. For example, in the time frame from 1998 to 2018, the total number of casualties due to earthquakes was larger than 846 thousand people, and the recorded economic losses were about USD 661 billion. At present, there are no earthquake precursors that can be used to trigger a warning. However, some studies have analyzed land surface temperature (LST) anomalies as a potential earthquake precursor. In this study, a large database of global LST data from the Geostationary Operational Environmental Satellite (GOES) and AQUA satellites during the whole year 2020 has been used to study the LST anomalies in the areas affected by earthquakes. A total of 1350 earthquakes with a magnitude larger than M4 were analyzed. Two methods widely used in the literature have been used to detect LST anomalies in the detrended LST time series: the interquartile (IQT) method and the standard deviation (STD). To the authors’ knowledge, it is the first time that the confusion matrix (CM), the receiver operating characteristic curve (ROC), and some other figures of merit (FoM) are used to assess and optimize the performance of the methods, and to select the optimum combination that could be used as a proxy for their occurrence. A positive anomaly was found a few days before the studied earthquakes, followed by the LST decrease after the event. Further studies over larger regions and more extended periods will be needed to consolidate these encouraging results.

1. Introduction

Earthquakes are one of the most destructive natural disasters, causing numerous casualties and economic losses. Between 1998 and 2018, the total number of casualties due to earthquakes was larger than 846 thousand people [1], and the reported material losses were roughly USD 661 billion [2].

Seismic events are unavoidable, but their effects could be mitigated if a warning is given in advance [3]. Using satellite data, remote sensing techniques can be exploited to make maps of zones exposed to specific hazards, including fire, hurricanes, floods, etc. However, earthquakes are unpredictable, and more efforts are required to prevent and monitor later seismic activities using earth observation (EO) data [4]. Compared to seismograph methods for earthquake monitoring, EO techniques have many advantages. For example, ground stations devoted to seismic event monitoring are spatially sparse; therefore, their ability to monitor crustal movements is limited. However, high spatio-temporal resolution data provided by satellite instruments can be used to overcome this limitation.

Several approaches have attempted to detect these land surface temperature (LST) anomalies using thermal infrared (TIR) datasets acquired from polar and geostationary satellites [5,6,7]. TIR data provide continuity over long periods, allowing a more precise calculation of the LST temporal anomalies over large regions. These anomalies are computed from the detrended LST values that are larger than a coefficient ‘C’ multiplied by the annual STD or IQT in the same position. More detailed explanations are provided in Section 2.3. Furthermore, satellite TIR data can depict extensive structures and short (from days to weeks) variations of LST anomalies over active faults [8]. Moreover, the movement of faults in earthquake-prone areas frequently affects micro-surface fracture growth, followed by variations in the Earth’s electric and magnetic fields [9]. These variations lead to pre-seismic LST anomalies in earthquake zones.

The main innovations of this work are: (1) the LST data pre-processing and analysis, especially the LST data extraction, quality and nighttime filtering, and LST data fusing and detrending over 1-year time series including seasonal effects as explained in Section 2.1 (Step 5), and (2) the use of techniques from detection theory such as the CM, the calculation of the ROC to optimize the detection threshold. Moreover, 1350 recent earthquakes of magnitude 4 or larger were studied for 2020 with consistent results.

Observations of thermal anomalies have been reported to precede earthquakes worldwide since the 1980s [5,10]. However, the physical relationship between the detected LST anomalies and the earthquakes has not yet been established [11]. Small earthquakes grow into strong ones until they spread along the fault line surface. Distortions are associated with the physical preparation of earthquakes; then, complex premonitory phenomena can occur in seismic dynamic regions under influential tectonic forces. The growing pressure energy establishes many changes in the Earth’s physical-chemical processes, appearing before the earthquake [6]. Unfortunately, most LST anomalies may have origins different from earthquakes and must be discarded [12].

In this study, after promising results presented in other previous studies, the LST variations in nearby areas affected by earthquakes have been analyzed [4] using the well-known interquartile (IQT) and standard deviation (STD) methods. The LST time series anomalies have been computed from the ABI/GOES and MODIS/AQUA data. Positive anomalies, i.e., detrended LST larger than usual, occurred a few days before the earthquakes, and a negative one, i.e., detrended LST less than normal, on 18 February 2018, and 2 April 2018, earthquakes in Mexico and Bolivia, respectively.

The Jiuzhaigou earthquake was studied in [8]. The TIR images and total electron content (TEC) fluctuations around the epicenter were recorded, encountering a high correspondence between the earthquake occurrence and the TEC anomalies.

A method was applied in [13] to find localized spatio-temporal fluctuations of LST using geostationary satellite data. The statistical results of the same regions over several years were also compared, and a distinct and repeated LST pattern was found in the case of an earthquake.

Significant pre-seismic TIR anomalies were also found 2 to 22 days before the earthquakes [14]. Thermal anomalies before large earthquakes are reported in [15], notably the Wenchuan earthquake in 2008 with M_w = 8.

The variations of in situ surface and atmospheric data from the Defense Meteorological Space Program (DMSP) Special Sensor Microwave Imager/Sounder (SSMIS) F17 satellite were studied in [16] for 20 days before two strong earthquakes of magnitude M_w > 7. A correlation among earthquakes, LST anomalies, and the air temperature was found during the seismogenic period.

The IQT, wavelet, and Kalman filter methods were used in [17] to detect LST anomalies before strong earthquakes (M_w > 6). Results show that the IQT was able to detect high anomaly values. In addition, the wavelet was sensitive to sudden changes, and the Kalman filter significantly detected minor LST anomalies.

A lithosphere–atmosphere–ionosphere–magnetosphere model was proposed in [18] to detect the pre-seismic signal and explain the process of air ionization before the earthquake.

Despite these evidences, pre-seismic LST anomalies are still not considered a reliable precursor of earthquakes; detection methods need to be improved, data should be made available to investigators for reproductivity purposes, and techniques to compute earthquake-related precursors need to be well standardized. For example, topography, soil moisture, vegetation cover, and other atmospheric phenomena might also impact pre-seismic thermal anomalies. This manuscript represents an attempt to improve the LST anomalies detection methods, and to select the optimum detection parameters. This fact led us to formulate the following research questions:

• How much can we trust the LST anomalies as earthquake precursors?
• What is the receiver operating curve of the detector, its optimum parameters, and the optimum probability of detection, false alarm, and other figures of merit?

This study aims at calculating and detecting earthquake LST anomalies using the STD and IQT methods in time series data [13] from the advanced baseline imager (ABI) [14] on GOES and from MODIS [15] on the AQUA satellite to inspect if any occurrence of LST anomalies before large, land-based earthquakes are detectable [19]. Only the data measured between midnight and 6:00 am local time are used to mitigate the diurnal cycle. The essential advantage of using multi-satellite data sources is the capability to generate accurate global and gap-free LST maps by data fusion of ABI/GOES 16 and 17, as explained in Section 2.1 (Step 4). Furthermore, earthquake events with a shock magnitude (M_w) larger than 4 have been studied globally during the year 2020. A few days or weeks before the seismic shocks, the LST anomaly of the pixels around the epicenters can be up to 1 °C to 4 °C [18]. This result is understood as the Earth’s crust’s thermal flux in seismically dynamic parts [20].

The rest of this manuscript is organized as follows. Section 2 describes the materials and methods used to process and calculate the LST anomalies. Section 3 shows the results, and Section 4 provides a discussion. Finally, Section 5 presents the main conclusions of the paper.

2. Materials and Methods

This section explains the methodology used to acquire, pre-process, and calculate the LST anomalies, including the EO datasets.

2.1. General Architecture for Data Acquisition and Processing

Figure 1 shows the LST data processing architecture from the data acquisition to the visualization of the predicted earthquake risk maps.

Step 1: Selection of relevant LST data sources [21], including MODIS/AQUA and ABI/GOES (East and West). These instruments have different spatial (SPR), spectral (STR), and temporal resolutions (TMR) [22], as detailed in

Table 1. LST Satellites Data Sources.

Satellite	Sensor	Spatial Resolution	Product	Input/Day	Rate/Day	Output/Day
GOES	ABI	4 km	LSTF	100 MB	60 files	95 MB
AQUA	MODIS	1 km	MYD11	1 GB	300 files	4 GB

.

Step 2: Data acquisition from the Earth Observing System Data and Information System (EOSDIS) of the National Aeronautics and Space Administration (NASA) [23]. This operation takes approximately 15 to 30 min, depending on the data transfer rate and the internet bandwidth [24], depending on the satellite sensor’s data latency, and the availability of the data center and data transmission [21].

Step 3: Pre-processing of the acquired LST data in 2020. First, the LST data is extracted from the Network Common Data Form (NetCDF, [25]) files using the SAT-ETL-Integrator software [26], and data filtering is conducted for land-only LST measurements. Erroneous, cloudy, and poor-quality data are filtered out, as described in

Table 2. LST Data Quality Products.

LSTF-M6	DQF	0	2	4	8	16
		Good retrieval valid input data, type, and LST clear conditions	Wrong or missing input data	Cloudy conditions	Degraded pixels	Invalid-water surface type
MYD11_L2	QC	0	1	2	3	-
		Pixel produced, good quality, cloud-free pixel	Pixel produced, unreliable quality, missing pixel	Pixel not produced due to cloud effects, fairly calibrated	Cloudy pixel not produced, poor calibration, processing skipped

[22]. Then, LST is converted from Kelvin (K) to degrees Celsius (°C), stored, and serialized into comma-separated values (CSV) files for each day [26]. The yearly total size of the processed data sumps up to 2 TB. Parallel and distributed tools were implemented to speed up the processing time using SAT-Hadoop-Processor, the distributed version of the SAT-ETL-Integrator. It is a new distributed software for remote sensing data pre-processing and ingestion using cloud computing technology, specifically OpenStack. The developed software discarded unneeded daily files and removed erroneous and inaccurate datasets. This parallel processing optimized the total execution time by 90%.

Step 4: Filtering the nighttime datasets between 00:00 and 06:00 am local time only to discard the diurnal cycle, i.e., the effect of variable sunlight and radiation on the LST values. Afterwards, data fusion (Only ABI/GOES East/West data, as they have the same spatial, radiometric, and spectral resolutions. MODIS/AQUA data is not fused with ABI/GOES data) aggregates measurements with the same location (longitude and latitude) and local time. Some figures, such as the average, standard deviation, and minimum and maximum values, are then evaluated. Finally, data are sorted and arranged chronologically in a time series.

Step 5: Computation of the average LST inside a region defined by the strain radius (SR), as explained in Section 4. The data processing showed that most LST anomalies occurred the week before the earthquake. Therefore, the average LST values within the SR region are detrended with breakpoints (BP) every 7 days to remove seasonality effects.

Step 6: Calculating and detecting the LST anomalies based on the STD and the IQT methods, as detailed in Section 3. Once the LST anomalies are calculated and detected, they are correlated and linked to the most probable earthquakes, as explained in Figure 3. Next, the CM is calculated to identify the optimum detection parameters. Finally, the correlation is visualized in charts and maps.

2.2. Input Data

LST is the radiative surface temperature of the land resulting from infrared radiation. Some factors can influence the derivation of LST, such as temperature variations within the pixel. Sub-pixel inhomogeneities in temperature and land covers, such as emissivity, atmospheric temperature, humidity, and clouds, can also affect the quality of LST measurement. LST is an important indicator of the energy flux between the atmosphere and ground, i.e., the Earth’s radiation budget.

Table 1 shows the satellites and sensors used in this study to provide the LST data. Their spatial resolution (SPR) ranges between 1 and 4 km, and their temporal resolution (TMP) resolution is between 1 h and 12 h, depending on the orbit [27]. In this study, global data is collected for 2020 [28].

2.2.1. LST Datasets from ABI/GOES

LST from ABI/GOES longwave infrared spectral channels, notably bands 14 (11.2 µm) and 15 (12.3 µm), with a NeDT of 0.1 °C [29]. Data were acquired from the Comprehensive Large Array-data Stewardship System (CLASS). LST dataset 2020 [30] was obtained from GOES 16 and 17 satellites. The product is available at a nominal SPR of 4 km, and a TMP of 24 images per day, although only seven images from midnight to 06:00 am were processed daily to avoid diurnal effects. Its spatial coverage encompasses the Americas. Quality information is also provided based on the total column moisture content, cloud coverage, and satellite positioning. Table 2 includes the ABI product’s data quality Flags (DQF). DQF ranges between 0 and 32; in our study, only data sets with a DQF equal to 0 are included, meaning that the data are good, have valid input, land type, and have cloud-free conditions.

2.2.2. LST Datasets from MODIS/AQUA

AQUA satellite is in a sun-synchronous orbit (SSO) at 1:30 pm local time of the ascending node. MODIS uses 36 spectral bands from 0.41 to 14.2 μm to sense the Earth’s surface [31] with a NeDT 0.25 °C in most channels [31]. The AQUA MYD11_L2 Global 1 km SIN Grid LST night has been used (Table 1) [32]. The LST values are also computed using a split-window algorithm. Table 2 includes the quality control (QC) of the MODIS product. In this study, only datasets with a QC equal to 0 were used, meaning that the pixel was produced with good quality and in clear conditions.

2.2.3. Earthquakes Datasets from the USGS Ground Stations

The United States Geological Survey (USGS) Earthquake Hazards Program (EHP) is a part of the National Earthquake Hazards Reduction Program (NEHRP) [33]. Earthquakes recorded within a 24 h window are clustered together. In this study, 1350 earthquakes with M_w ≥ 4 were studied, but only 3 with M_w > 6 (

Table 3. Details of earthquakes M_w > 6 in 2020 [33] that are presented in Figure 4. (From a total of 1350 earthquakes analyzed).

ID	Date	Location	Epicenter			Mw	SR (km)	Land Cover
			Lat (Deg)	Lon (Deg)	Depth (km)
ID977	2020-01-24	Ankara, Turkey	38.43	39.06	10	6.7	760	Grassland
ID5086	2020-05-15	San Jose, CA, USA	38.17	−117.85	2.7	6.5	623	Shrub
ID9629	2020-09-06	Coquimbo, Chile	−30.34	−71.49	30	6.3	511	Shrub

) are illustrated as examples in this paper to explain the methodology. Figure 2 shows the locations of the sample earthquakes listed in Table 3.

2.3. LST Anomaly Calculation

This section explains the two most widely used statistical methods to calculate LST anomalies in time series, notably the STD and the IQT methods.

2.3.1. STD Method

LST values further away than ±1.5. σ from the mean were flagged [34,35], and ±1.5. STD and ±2. STD thresholds were used in [36]. The LST deviation is calculated as LST deviation = LST – LST mean, where the LST mean is the average of detrended LST inside the strain radius (SR). For the earthquake precursor, the following thresholding is considered. The LST anomaly is detected if the flowing condition is met:

| LST deviation | ≥ C × STD,

(1)

where “C” is a coefficient that multiplies the STD of the annual detrended LST to detect LST anomaly.

In this study, 9 values for the C coefficient have been considered [0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1, 1.5, 2]. The first step calculates the STD values for each coefficient during 2020. Values satisfying Equation (1) are flagged as LST anomalies in this time series. A sensitivity analysis of this approach showed that this threshold captures fluctuations between 1 °C and 3 °C, as compared to LST neighboring background. The link between LST anomalies and studied earthquakes will be explained in Section 2.4.

2.3.2. IQT Method

The IQT is a statistical parameter that measures the numerical spread and data variability by dividing the dataset into quartiles. Any dataset is ordered and splits into 4 distinct intervals, including [0–0.25], (0.25–0.50], (0.50–0.75], and (0.75–1.00]. In this method, 9 values of C have been considered: 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1, 1.3, 1.5, and LST anomalies are the values satisfying the following condition.

| LST deviation | ≥ C × IQT

(2)

2.3.3. LST and Earthquakes Aggregation

Previous studies detected LST anomalies even 500 km away from the epicenter [37]. However, in [38], it is claimed that because the latent alarm zones are too large, observed phenomena farther than 60 km from the earthquake epicenter are not useful as precursors, even if they were spatially correlated to the epicenters. The SR is the radius of the circle of the earthquake’s groundwork zone, and it is related to the moment magnitude (M_w) by [39].

SR (km)= 10 ^{0.43 × Mw}

(3)

For instance, on 24 January 2020, an earthquake occurred in Ankara, Turkey, with M_w = 6.7, and SR is ≈720 km. This study uses the SR (Equation (3)) to calculate the averaged LST for every pixel. The neighboring radius (NR) is used to search for near earthquakes to expand the correlation between LST anomalies and earthquakes. In this study, NR is computed as:

NR max (km)= 2 × SR

(4)

Figure 3a shows the aggregation process between the epicenters and the LST anomalies in those pixels. First, neighbor earthquakes are searched within the NR and associated with positive LST anomalies up to seven days before. As a result, multiple aggregations are possible; the one selected is based on the strongest earthquakes, the most recent one, the closest epicenter, and the largest LST anomaly, in this logical order. For instance, the San Jose, USA earthquake (ID5086) showed two positive LST anomalies, 1 °C and 2.8 °C, recorded 3 and 6 days before the earthquake, respectively. Accordingly, the first anomaly was associated with the earthquake because it is the most recent one.

In 2020, the same area was examined with or without an earthquake occurrence. This approach allows testing the effect of anomalies in periods without earthquakes. Prior studies have frequently tested various periods, but only before and after the earthquake. This leaves the possibility that similar irregularities may occur throughout the year(s) without any associated earthquake [40]. In our case, we have found more than 60,000 LST anomalies from 370,000 without any associated earthquake, as also found in [41,42]. Furthermore, it has been shown that anomalies in the neighborhood of earthquake epicenters and everywhere in the study areas happen during the year, regardless of the existence of an earthquake.

2.4. LST Data Pre-Processing and Normalization

The high-level procedure to detect the LST anomalies as seismic precursors were outlined in Section 3. Here, more details on the LST processing are provided. Among the 1350 earthquakes studied, three examples corresponding to Table 3 are also illustrated. First, the thermal conditions around the epicenter’s nighttime are visually analyzed seven days before and four days after the earthquake. A graphical analysis of the thermal images shows the approximate number of days of the appearance of an LST anomaly. Finally, the plots of time series of LST anomalies are presented for the whole year to show other examples in the same region but on different dates.

A spatial low pass filter (LPF) with a radius equal to the SR is applied to discard spatial patterns. This operation eliminates local inhomogeneities of pixels, humid conditions, or a heatwave spreading over large zones. Afterwards, only local anomalies remain highlighted in the subsequent time series. Before an LST anomaly can be designated as a consistent precursor, it should be considered whether or not it is an artificial anomaly associated with an investigated event or a random anomaly. Before strong and moderate seismic events, earthquake precursors appear at different distances. The zone of strain accumulation is estimated following the definition of earthquake preparation area as described in [42].

Figure 4 illustrates the above procedures with three examples of the studied epicenters (ID5086, ID9629, and ID977), as detailed in Table 3. On top, it shows the fused and LPF the LST 4 days before, and 1 day after; the cyan star represents the epicenter of the earthquakes, and the dashed circle is the region defined by the SR. The orange ellipse is the central studied epicenter. Before each earthquake, an LST increase around the epicenters occurs 2 to 4 days before the seismic shock. Then, the LST cools down after the earthquakes. Figure 4(a1–a4, b1–b4, and c1–c4) show the original, and the detrended LST for each subplot during 2020 to show all the LST anomalies in the studied area. Differences in seasonal detrended LST are due to the latitude (north vs. south, summer vs. winter), which is more stable and independent of seasonal effects.

Figure 4(a1–a4) show the calculated LST anomalies used for the San Jose earthquake and the STD and IQT methods. Each chart has a main epicenter and neighbor epicenters closer than the NR max. Thus, a positive LST anomaly can be noticed 1 to 4 days before most studied earthquakes, as shown in Figure 4(a5). Recent LST anomalies to earthquakes have been linked in this study because, logically, seismic events rapidly impact the LST. Thus, for instance, +1 °C and +3 °C were recorded 1 and 4 days before earthquake ID5086, as shown in Figure 4(a5). Accordingly, the closest one in time, +1 °C, was associated with the earthquake.

Figure 4(b1–b6) show the analyzed, fused, and LPF LST around Coquimbo (Chile) between 3 and 7 September 2020. The epicenter is shown with a cyan star. Abnormal LST increments of 1 to 4 °C are detected around the epicenter, increasing to a maximum build-up of a thermal anomaly 2–4 days before the seismic events. Then, the LST anomaly cools down immediately after the earthquakes.

Figure 4(c1–c6) illustrate the processed AQUA/MODIS LST data near Ankara, Turkey, between the 20 and 25 of January 2020. The cyan star represents the epicenter. We noticed a slight increase in LST 3 days before the earthquake. Still, a significant positive LST anomaly was reported on the earthquake day, and then the LST cooled down the day after.

2.5. LST and Earthquakes Correlation Analytics Using the Full Confusion Matrix

The CM characterizes the performance of a detector. It helps measure other parameters such as recall, precision, specificity, and accuracy. Each row of the matrix represents the cases in an actual class. Conversely, each column denotes the samples in a predicted class or vice versa [11]. The four components of CM in this study are defined as follows:

• True Positive (TP) is the number of correct predictions: An earthquake occurs, and an LST anomaly occurs.
• False Negative (FN) is the number of incorrect predictions: An earthquake occurs, and there is no LST anomaly.
• False Positive (FP) is the number of incorrect predictions: An earthquake does not happen, but there is an LST anomaly.
• True Negative (TN) is the number of correct predictions: An earthquake does not occur, and there is no LST anomaly.

Figure 3a shows the flowchart describing the CM calculation algorithm. Samples were divided into two parts: the CP (an earthquake actually occurs) and CN (no earthquake occurs). In CP, only samples with M_w ≥ 4 are retained, then every earthquake is associated with its best LST anomaly, as explained in Section 3. The LST anomalies could either be positive or neutral. If a positive LST anomaly and a neighboring earthquake are found for up to 7 days, the TP is incremented. In contrast, the FP is incremented. If a neutral LST anomaly and a neighboring earthquake are found for up to 7 days, the FN is incremented. In contrast, the TN is incremented. Finally, the algorithm calculates the other metrics using the formulas provided in

Table 4. Details of the Confusion Matrix and Formulas of FoM Parameters.

Condition Positive (CP)	TP	FN	True Positive Rate (TPR) = TP/P	False Negative Rate (FNR) = FN/P
Condition Negative (CN)	FP	TN	False Positive Rate (FPR) = FP/N	True Negative Rate (TNR) = TN/N
Prevalence (P) = CP/CP + CN	Accuracy (ACC) = TP + TN/CP + CN		Positive Likelihood Ratio (LR+) = TPR/FPR	Negative Likelihood Ratio (LR−) = FNR/TNR
G-mean $=\sqrt{\frac{\mathrm{T}\mathrm{P}\ast \mathrm{T}\mathrm{N}}{\left(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}\right)\ast \left(\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{P}\right)}}$		Kappa Coefficient $=\frac{2\ast \left(\mathrm{T}\mathrm{P}\ast \mathrm{T}\mathrm{N}\right)-\left(\mathrm{F}\mathrm{N}\ast \mathrm{F}\mathrm{P}\right)}{\left(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}\right)\ast \left(\mathrm{F}\mathrm{P}+\mathrm{T}\mathrm{N}\right)+\left(\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{N}\right)\ast \left(\mathrm{F}\mathrm{N}+\mathrm{T}\mathrm{N}\right)}$		Diagnostic Odds Ratio (DOR) = (LR+)/(LR−)

.

3. Results

This section describes the results obtained after the pre-processing, processing, and statistical calculation of the ABI/GOES and MODIS/AQUA data. First, the numerical correlation analysis is provided. Then, the results of the CM are presented, and finally, the outcome of the CM for each M_w.

3.1. Statistical Analysis of LST Anomalies and Earthquakes Correlations

Table 5. Average LST Anomalies (°C) for TP as a Function of Mw and FP.

C/Mw	TP (Mw = 7–7.9)		TP (Mw = 6–6.9)		TP (Mw = 5–5.9)		TP (Mw = 4–4.9)
Sensor	M	A	M	A	M	A	M	A
STD method: Mean LST anomalies (°C)
0.7	-	0.56	1.41	3.1	2.69	1.46	2.9	1.83
0.75	-	0.51	1.31	3	2.54	1.38	2.71	1.75
0.8	-	0.45	1.21	2.89	2.39	1.31	2.53	1.68
0.85	-	0.4	1.1	2.79	2.24	1.24	2.36	1.61
0.9	-	0.35	1	2.69	2.1	1.17	2.19	1.54
0.95	-	0.29	0.89	2.59	1.96	1.1	2.03	1.48
1	-	0.24	0.79	2.5	1.82	1.04	1.87	1.41
1.5	-	-	0.24	1.7	0.79	0.58	0.74	0.9
2	-	-	0.03	1.12	0.27	0.29	0.26	0.58
IQT method: Mean LST anomalies (°C)
0.7	-	0.42	0.83	3.02	2.16	1.44	2.02	1.81
0.75	-	0.36	0.7	2.92	1.98	1.37	1.82	1.74
0.8	-	0.29	0.59	2.82	1.81	1.3	1.64	1.67
0.85	-	0.23	0.52	2.72	1.65	1.23	1.47	1.59
0.9	-	0.16	0.46	2.62	1.5	1.16	1.32	1.52
0.95	-	0.1	0.4	2.52	1.35	1.1	1.18	1.46
1.1	-	-	0.23	2.26	1	0.93	0.84	1.27
1.3	-	-	0.07	1.98	0.65	0.75	0.54	1.06
1.5	-	-	0.01	1.71	0.42	0.59	0.36	0.89
Other metrics
Days before	-	−6	−1.86	−2.2	−3.01	−2.4	−2.81	−2.7
Count	1	1	22	11	238	93	2972	1255

Cells with (-) are missing data due to gaps or poor-quality measurement data. (A) means ABI and (M) means MODIS.

shows the average LST anomalies for each M_w. It can be noticed that the detected LST anomalies are almost similar using the STD and IQT methods with the same coefficient. In addition, the values of LST anomalies decrease when coefficient C increases. The ABI/GOES LST anomalies are higher than the MODIS/AQUA ones with high M_w and vice versa [43,44]. This could be explained because ABI/GOES and MODIS/AQUA have different coverages and spatial resolutions [45]. MODIS scans the whole globe, while ABI scans only America. The number of samples and pixel size are very different [46].

The ABI/GOES and MODIS/AQUA have different SPR and TMR, affecting the values of the calculated LST anomalies. Furthermore, the measured LST anomalies are higher with small M_w and lower with strong earthquakes. LST anomalies frequently appear earlier than smaller ones. Finally, the average LST anomalies in FP values recorded from ABI are higher than MODIS.

3.2. Output of the CM and Visualization of the ROC Curves

This section comprises the results obtained after the CM calculation. First, the numerical analyses of the CM acquired from ABI/GOES and MODIS/AQUA satellite instruments are provided. Then, the results of the partial CM for each M_w are presented. Finally, the ROC curves are shown.

3.2.1. CM of ABI/GOES for M_w > 3

Figure 5 shows the results of the CM measured from ABI/GOES. It can be noticed that all metrics are almost equal for both the STD and IQT methods using the same ‘C’ coefficient. The ACC, TPR, FNR, and TNR increase with increasing the ‘C’ coefficient, while the FPR decreases. Accordingly, the LST anomaly detection increased for the ‘C’ coefficients, but the rate of false alarms is also significantly high. Few LST anomalies are detected with large ‘C’ coefficients, with high confidence and less noise. The diagnostic odds ratio (DOR) measures the effectiveness of a diagnostic test. It ranges from zero to infinity. However, it is larger than one for practical tests, and the higher the DOR, the better the test performance is. The geometric mean (G-mean) indicates the balance between the classification performances of majority and minority classes. A low G-mean means a low prediction of the TP cases, even if the TN cases are correctly classified, and vice-versa. The kappa coefficient calculates the agreement between predicted and valid values. A kappa coefficient value of 1 denotes perfect agreement, while a value of 0 represents no agreement at all. In our case, the DOR (≈13), G-mean (≈0.78), and kappa coefficient (≈0.56) reach the maximum value for C = 0.7, which confirms the selection of this value as the optimum detection parameter for both methods.

3.2.2. CM of MODIS/AQUA for M_w > 3

From Figure 6, it can be noted that for the STD method, using C = 0.7, the TPR, G-mean, kappa coefficient, and the DOR reach the maximum values: 0.81, 0.79, 0.57, and 13.56, respectively. However, the FPR is also high at 0.24, meaning a high prediction sensitivity at the expense of a high rate of false alarms with small coefficients. The IQT method also uses C = 0.7, the TPR, G-mean, kappa coefficient, and the DOR reach only a maximum value, correspondingly, 0.67, 0.76, 0.57, and 11.51. In contrast, the lower FPR recorded was 0.16. Accordingly, for M_w ≥ 4, a C = 0.7 is the optimum method parameter.

3.2.3. CM as a Function of the Earthquake Magnitude (M_w = 7, 6, 5, and 4)

Table 6. Optimum Coefficients of The CM for Each M_w.

	TPR		FPR		FNR		TNR		DOR
Sensor	A	M	A	M	A	M	A	M	A	M
Mw = 7 (GOES CP = 1 and AQUA CP = 1)
0.7-STD	1.00	1.00	0.18	0.24	0.00	0.00	0.82	0.76	-	-
0.7-IQT	1.00	1.00	0.18	0.17	0.00	0.00	0.82	0.83	-	-
Mw = 6 (GOES CP = 42 and AQUA CP = 79)
0.7-STD	0.88	0.94	0.18	0.24	0.12	0.06	0.82	0.76	34.7	48.0
0.7-IQT	0.86	0.85	0.18	0.17	0.14	0.15	0.82	0.83	26.4	26.5
Mw = 5 (GOES CP = 185 and AQUA CP = 461)
0.7-STD	0.8	0.80	0.18	0.24	0.2	0.20	0.82	0.76	18.7	13.3
0.7-IQT	0.8	0.75	0.18	0.17	0.2	0.25	0.82	0.83	17.6	14.1
Mw = 4 (GOES CP = 2465 and AQUA CP = 5338)
0.7-STD	0.73	0.80	0.18	0.24	0.27	0.20	0.82	0.76	12.7	13.4
0.7-IQT	0.74	0.70	0.18	0.17	0.26	0.30	0.82	0.83	12.3	11.2

Cells with (-) are missing data due to gaps or bad quality measurement data. (A) means ABI and (M) means MODIS.

includes the CM measured from ABI/GOES and MODIS/AQUA data for earthquake magnitudes from 4 to 4.9, 5 to 5.9, 6 to 6.9, and 7 to 7.9. It can be noticed that the optimum thresholds do not depend on the earthquake magnitude. Furthermore, as expected, the TPR, TNR, and DOR increase as earthquake magnitude increases, indicating more detectability. Furthermore, the performance is dominated mainly by M_w = 4 earthquakes, which are much more abundant. The DOR improves because the FNR decreases with increasing magnitude, while TPR, FPR, and TNR are roughly constant. There is no monotonic behavior based on the magnitudes because the number of studied earthquakes by magnitude is very different for different magnitudes. For instance, there was only one earthquake with M7, hundreds of M5, and thousands with M5 and M4. Thus, we could not deduce a such a rule.

3.2.4. ROC Curves of ABI/GOES and MODIS/AQUA for M_w ≥ 4

A ROC curve is a plot illustrating the performance of a binary classifier system as a function of the threshold. The probability of a false alarm (FPR) is on the x-axis, while the probability of detection (TPR) is on the y-axis. It can be used to find the optimum threshold that minimizes the distance between the perfect classifier (Pfa, Pdet) = (0,1) and the ROC.

From Figure 7, it can be noticed that the four ROC curves calculated from ABI/GOES and MODIS/AQUA data using the STD and IQT methods are very similar and have almost the same shape. With C = 0.1, the TPR = 1, and the FPR = 1, while with C = 3, the TPR = 0, and the FPR = 0.

When C = 0.7, the shortest distance from the ROC to the point (0,1) is achieved, ranging between 0.3 and 0.36 depending on the data used and method; as a result, the optimum C for earthquake prediction is 0.7.

4. Discussion: Correlation between LST and Earthquakes for Strong Earthquakes (M_w> 6)

The proposed approach correlates earthquake occurrence and previous LST anomalies. As can be seen, each method (STD and IQT) leads to a minor different LST anomaly, despite using the same data source. The current research used nighttime ABI/GOES and MODIS/AQUA LST data. During the nighttime, small but detectable, positive LST anomalies around 1 to 4 °C were observed before the earthquakes. A negative LST anomaly is expected during 1 to 3 days after the earthquake.

For the studied earthquakes, LST generally increases along the fault zone of the epicenter region, and the location of LST anomalies is within the SR of each studied epicenter.

Table 7. Measured LST Anomalies (°C) from GOES and AQUA M_w > 6 Studied Earthquakes In 2020.

C/ID	ID6704		ID3436		ID5086		ID334	ID9629		ID977
Satellite	A	G	A	G	A	G	A	A	G	A
0.7-STD	-	0.56	2.4	0.73	-	1.04	0.85	-	0.27	1.05
0.7-IQT	-	0.42	2.11	0.64	-	0.94	0.35	-	-	0.67
Days before	-	6	1	1	-	3	5	-	2	4

Cells with (-) are missing data due to gaps or bad quality measurement data. (A) means AQUA, and (G) means GOES.

summarizes the LST anomalies for the six strongest earthquakes during 2020, with small LST anomalies from 0.4 to 2.4 °C occurring from 1 to 6 days before the earthquakes. Still, earthquake prediction is very complex and debatable since the observed changes are tiny, even around the epicenter and surrounding regions.

The Oaxaca earthquake, Mexico (ID6704), was the strongest land earthquake in 2020. It occurred on 23 June 2020, with an Mw = 7.4 at a depth = 20 km in a forested area. A slight positive anomaly of just 0.6 °C was detected 6 days prior to the seismic event. The LST increased up to a maximum of 24.4 °C, then decreased gradually to 23 °C the day after the earthquake. The Portland earthquake, USA (ID3436) occurred on 15 May 2020, with M_w = 6.5 at a depth =10 km in a shrubland area. The calculated SR is 626 km. A positive LST anomaly of approximately 2 °C was detected 3 days before the seismic event, with a maximum recorded LST of 18.3 °C one day before the earthquake. The LST then cooled down to 15.6 °C 1 day after. The same observation was detected in the San Jose earthquake, USA (ID5086). However, a positive LST anomaly was observed only one day before the seismic shock with a maximum LST equal to 4 °C, and it decreased progressively to 2.5 °C the day after the earthquake. Regarding the Chile (ID9629) earthquake, a small prior positive anomaly equal to 0.3 °C was detected 2 days before the earthquake.

In the MODIS/AQUA LST night, positive anomalies around 1 °C were detected 4 and 5 days before the Anadyr, Russia (ID344) and Ankara, Turkey (ID977) earthquakes. The maximum LST recorded were −31.8 °C and −4.1 °C, respectively, then the LST cooled down steadily to −38 °C and −5.3 °C, respectively, the day after the main earthquake.

The anomalies are small, within the range of 0~2 °C. In addition, other examples are in the background (not shown on paper because we do not have enough space) that have been associated with high LST anomalies (>2 °C).

Many researchers have observed abnormal changes in LST anomalies before strong earthquakes, which can be considered potential earthquake precursors. However, many studies argue that the changes in these parameters are not related to earthquake activities, since even stronger LST anomalies also appear in the absence of earthquakes.

The relationship between longtime LST data series and dates of earthquakes required more investigations to determine whether the anomalies were related to the earthquake. Our results confirm that the STD and IQT methods with proper selection of the C parameter (Copt = 0.7), the proper spatial filtering, and the temporal detrending are capable of detecting LST anomalies linked to the occurrence of earthquakes with a probability of detection (Pdet) = 0.67–0.81 and probability of false alarm (Pfa) = 0.16–0.24 depending on the method and sensor. This performance can be improved by repeating the test on consecutive days.

5. Conclusions

This study analyzed the LST anomalies for 1350 earthquakes that occurred in 2020 using time series data of ABI/GOES [47] and MODIS/AQUA. To gain a wider perspective, this period could be extended to a few years to enhance the quality of correlation and to cover global extents with more different land covers and regions. Only three earthquakes with M_w > 6 have been visualized because they clearly show the correlation. The pre-processing of this massive volume of data was accomplished efficiently using the SAT-ETL-Integrator [26] and SAT-Hadoop-Processor software [27].

LST anomalies were identified using the IQT and STD methods. However, the efficiency of integrating the above-mentioned methods to detect anomalies in LST variations in the strong earthquakes studied has been presented. This has been illustrated by the San Jose, Chile, and Ankara earthquakes (Figure 4). As a result, the detected LST anomalies can be trusted as a precursor for earthquake events. The IQT and STD methods were applied with different threshold coefficients, and the CM and the ROC curves were generated for each C to select the optimum threshold, concluding that C_opt = 0.7. However, it is acknowledged that some anomalies noticed near the earthquake may not be related to the earthquake but may be masked by vegetation cover, topography effects, or other climatological conditions.

The results of the pre-earthquake LST anomalies detected before strong and moderate earthquakes (M_w ≥ 4) studied in this work may be used in the future in conjunction with other proxies for short-term earthquake risk prediction. In addition, we can reuse the same model for other parameters, notably ionospheric scintillation parameters [48,49].

At the time of writing the original manuscript, a sequence of earthquakes occurred between 12 and 19 September 2021, in La Palma Island (Canary Islands), preceding the volcanic eruption on 19 September. As illustrated in Figure 8a, b, there were hundreds of small earthquakes with M_w ≤ 3, located primarily in the southwest part, that progressively shifted towards the east after the volcano eruption. The analysis of AQUA/MODIS and TERRA/MODIS nighttime, high-quality, cloud-free conditions data showed an increase of the LST of up to ~7 °C from 2 September to 11 September, in the partially active part of the island (see Figure 9a). On 12 September, after the earthquakes started, the LST suddenly dropped by 4 °C (see Figure 9b). The rest of the time series uses data not corrupted (lava, clouds, etc.), and the LST remains approximately constant. These results illustrate the power of LEO satellites for measuring and detecting the LST anomalies on islands or small land bodies because of the higher spatial resolution of MODIS (1 km) compared to GEO satellites.

Other factors can produce thermal anomalies, e.g., human activities, but LST anomalies induced by earthquakes are larger in magnitude and impact vast areas (ten to a hundred kilometers). As a result, the percentage of false alarms due to these other sources is low and does not affect the result. Future work will address vegetation, climatological, human activities, and topographic effects to minimize the probability of false alarms and missed detections, and enlarge the data record’s temporal length.