RST Analysis of Anomalous TIR Sequences in Relation with Earthquakes Occurred in Turkey in the Period 2004–2015

: The paper provides, for the ﬁrst time, a long-term (>10 years) analysis of anomalous transients in Earth’s emitted radiation over Turkey and neighbouring regions. The RST (Robust Satellite Techniques) approach is used to identify Signiﬁcant Sequences of Thermal Anomalies (SSTAs) over about 12 years (May 2004 to October 2015) of night-time MSG-SEVIRI satellite images. The correlation analysis is performed with earthquakes with M ≥ 4, which occurred in the investigated period/region within a pre-deﬁned space-time volume around SSTA occurrences. It conﬁrms, also for Turkey, the possibility to qualify SSTAs among the candidate parameters of a multi-parametric system for time-Dependent Assessment of Seismic Hazard (t-DASH). After analysing about 4000 images (about 400 million of single satellite records), just 155 SSTAs (about 4 every 100 images) were isolated; 115 (74% out of the total) resulted in earthquake-related (false-positive rate 26%). Results of the error diagram conﬁrms a non-casual correlation between RST-based SSTAs and earthquake occurrences, with probability gain values up to 2.2 in comparison with the random guess. The analysis, separately performed on Turkish areas characterized by different faults and earthquakes densities, demonstrates the SSTA correlation with a dynamic seismicity more than with static tectonic settings.


Introduction
The optimism following the successful prediction of the large Chinese Haicheng earthquake (4 February 1975, M = 7.3) made the scientific community confident in achieving 'earthquake prediction' within about a decade. The evacuation of the citizens living in that densely populated area was ordered after hundreds of foreshocks, which had been registering for three days before the mainshock, and other geophysical precursors. However, this first success was followed by a series of failures and missed events [1], such as the strong earthquake that occurred 18 months later at about 450 km from Haicheng (M7.8 Tangshan), as well as other devastating seismic events all around the world (e.g., Northridge, 1994;Kobe, 1995;Wenchuan, 2008). These unsuccessful situations caused the progressive damping of the initial enthusiasm. This declining feeling was also emphasised by the failure of the probabilistic seismic hazard analysis (PSHA), used for decades in the field of seismic hazard. In fact, many disastrous events, such as Wenchuan (2008), Haiti (2010), and Tohoku (2011), shook areas much more than the levels foreseen in hazard maps [2]. However, in the last decade, a renewed interest in earthquake forecasting has spread in the form of a "well-disposed scepticism" as opposed to the "wishful data mining" of early researchers [3]. The scientific community has begun to turn its attention to the new research lines based on alternative data analysis methods hoping to improve our present capability to assess the seismic hazards in the short-medium term. During decades 1.
The continuously increasing stress field determines an extensive process of microcrack formation with a consequent increase in degassing activity together with deepwater and convective heat flow rising toward the surface; 2.
When the stress field becomes high enough locally to close the cracks and the earthquake occurrence is approaching, all the above processes are expected to reduce up to the time of the earthquake occurrence; 3.
At the time of the earthquake occurrence, because of a major crack opening in the rupture zone, a new increase in degassing activity (and related phenomena) is expected before a gradual return to normality.
All the phenomena occurring during phases (a) and (c) contribute to increasing TIR emissions because of the above-mentioned local greenhouse effect (increasing near-surface temperature) and/or the increase in TIR spectral emissivity consequent to the increase in soil humidity due to deep-water rising. For the same reason, during phase (b), a reduction (up to a complete disappearance) of TIR anomalies is expected because of the reduction (up to a complete interruption) of gas release and deep-water rise.
To this aim, a 12-year-long time series (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) of satellite images, collected by the Spinning Enhanced Visible and Infrared Imager (SEVIRI) sensor on board Meteosat Second Generation (MSG) geostationary satellites, is analysed. For the first time, the analysis is performed not only on the whole investigated area (Turkey and neighbouring areas) but also considers separately different fault systems in the Anatolian Peninsula and the Eastern Mediterranean region. In addition, we use a customised Molchan diagram to compare the achieved results to the ones achievable by a random space-time distribution of Significant Sequences of Thermal Anomalies (SSTAs, such as in [49,54]) identified by RST.

Tectonic Setting of the Investigated Area
Turkey and neighbouring regions in the eastern Mediterranean Sea are characterised by high seismicity and such rapid deformations that they are considered a natural laboratory, which is helping to better understand the complex processes in collision areas [58]. Turkey is in 10th place within the 50 most seismically active areas in the world [59].
The plate tectonic framework of the region, firstly described in the 1970s (see, for example, [60,61]), is governed by three major tectonic plates-African, Arabian, and Eurasiantogether with the Aegean and Anatolian minor plates. The tectonics of the area are strongly influenced by the northward subduction of the African plate below western Turkey and the Aegean region and the collision in the Caucasus and eastern Turkey [58].
Some authors [62], modifying an original scheme [63,64], identify the following five provinces showing different tectonic characteristics: North Anatolian Province (NAP), Eastern Anatolian Contractional Province (EACP), Central Anatolian Province (CAP), Western Anatolian Extensional Province (WAEP), and Southeast Anatolia Province (SEAP). EACP is under an S-N compressional tectonic regime, CAP is characterised by recent complex tectonic activity mainly with strike-slip faults, NAP, with very low microseismicity, is the area between the Black Sea and the North Anatolian Fault (NAF) system, SEAP, characterised by weak seismicity, represents the foreland of the Southeast Anatolian thrust zone, and WAEP is under an N-S extensional faulting mechanism. WAEP is further divided into smaller blocks such as the West Anatolia Graben (WAG) systems, the Isparta Angle (IA), and the Northwest Anatolia Transition Zone (NATZ).
The dextral NAF Zone, the sinistral East Anatolian Fault (EAF) Zone, the Hellenic (HA) and Cyprian Arcs (CA), as well as the sinistral Dead Sea Fault Zone (DSFZ), are the main plate boundaries, i.e., the most important tectonic structures [62,65] of the investigated area. Other relevant tectonic lineaments are the Middle Caspian Region-Kazakhstan fault system, the Eskişehir and Tuzgölü faults, which separate the CAP and WEAP, and the Caucasus thrust belt, representing the northern limit of the EACP. Figure 1 gives an overview of the above-mentioned tectonic provinces and main fault systems.
provinces showing different tectonic characteristics: North Anatolian Province (NAP), Eastern Anatolian Contractional Province (EACP), Central Anatolian Province (CAP), Western Anatolian Extensional Province (WAEP), and Southeast Anatolia Province (SEAP). EACP is under an S-N compressional tectonic regime, CAP is characterised by recent complex tectonic activity mainly with strike-slip faults, NAP, with very low microseismicity, is the area between the Black Sea and the North Anatolian Fault (NAF) system, SEAP, characterised by weak seismicity, represents the foreland of the Southeast Anatolian thrust zone, and WAEP is under an N-S extensional faulting mechanism. WAEP is further divided into smaller blocks such as the West Anatolia Graben (WAG) systems, the Isparta Angle (IA), and the Northwest Anatolia Transition Zone (NATZ).
The dextral NAF Zone, the sinistral East Anatolian Fault (EAF) Zone, the Hellenic (HA) and Cyprian Arcs (CA), as well as the sinistral Dead Sea Fault Zone (DSFZ), are the main plate boundaries, i.e., the most important tectonic structures [62,65] of the investigated area. Other relevant tectonic lineaments are the Middle Caspian Region-Kazakhstan fault system, the Eskişehir and Tuzgölü faults, which separate the CAP and WEAP, and the Caucasus thrust belt, representing the northern limit of the EACP. Figure 1 gives an overview of the above-mentioned tectonic provinces and main fault systems.

Data
The following sections give some details about the satellite data and seismic information analysed and used in this paper.

Satellite Data
The long-term analysis described in this paper is based on a 12-year-long series of night-time (time slot 00:00 GMT) MSG-SEVIRI satellite data. The time series consists of images acquired in the TIR spectral band centred at 10.8 µm (9.8-11.8 µm) from May 2004

Data
The following sections give some details about the satellite data and seismic information analysed and used in this paper.

Satellite Data
The long-term analysis described in this paper is based on a 12-year-long series of night-time (time slot 00:00 GMT) MSG-SEVIRI satellite data. The time series consists of images acquired in the TIR spectral band centred at 10.8 µm ( Figure 2). In this area, each SEVIRI pixel covers a surface of about 4.5 km (longitude) × 4.7 km (latitude).  Figure 2). In this area, each SEVIRI pixel covers a surface of about 4.5 km (longitude) × 4.7 km (latitude). The black line defines the area where the RST approach was applied to identify TIR anomalies on SEVIRI satellite images (SEVIRI cut area). The irregular shape of this area is due to the reprojection of the rectangular (215 rows × 476 columns) SEVIRI image portion from its satellite projection to a latitude-longitude grid. The red line (1 degree wider than the black contoured area) describes the area including all considered earthquakes (2484) before filtering. Background map: Wikimedia unlabelled map, ©OpenStreetMap contributors and available from [69,70].
A total of 3831 mid-night images represents the analysed data set (Ⴈ), which has been organised on a monthly basis.

Seismic and Tectonic Information
The performed analysis has taken account of seismic events with a magnitude greater than or equal to four (M ≥ 4) recorded in the catalogue of the United States Geological Survey [71]. The catalogue reports 2484 seismic events in the Turkish territory and neighbouring areas from April 2004 to October 2015. The actual area of the earthquake selection is one degree wider than the satellite cut area (see the red line in Figure 2). The distribution of seismic events shows a greater concentration along with the Hellenic Arc, the Isparta Angle, the East, and North Anatolian fault systems.
Since earthquake clusters could affect the results [72], following the approach described in [49] and before the correlation analysis, we filtered the events selecting among those that occurred on the same day within an epicentre-to-epicentre distance of 0.1 degree, the one with the highest magnitude. The resulting "filtered" catalogue, consisting of 2328 events, represents the reference seismic database used in the following analysis.

The Robust Estimator of Thermal InfraRed Anomalies (RETIRA)
Details about the methodology used in this paper and its previous applications are summarized in [17,73]. In this section, we want just to recall that the RST approach rests on the Robust Estimator of TIR Anomalies (RETIRA) index, whose more general expression is the following: The black line defines the area where the RST approach was applied to identify TIR anomalies on SEVIRI satellite images (SEVIRI cut area). The irregular shape of this area is due to the reprojection of the rectangular (215 rows × 476 columns) SEVIRI image portion from its satellite projection to a latitude-longitude grid. The red line (1 degree wider than the black contoured area) describes the area including all considered earthquakes (2484) before filtering. Background map: Wikimedia unlabelled map, ©OpenStreetMap contributors and available from [69,70].
A total of 3831 mid-night images represents the analysed data set (  Figure 2). In this area, each SEVIRI pixel covers a surface of about 4.5 km (longitude) × 4.7 km (latitude). The black line defines the area where the RST approach was applied to identify TIR anomalies on SEVIRI satellite images (SEVIRI cut area). The irregular shape of this area is due to the reprojection of the rectangular (215 rows × 476 columns) SEVIRI image portion from its satellite projection to a latitude-longitude grid. The red line (1 degree wider than the black contoured area) describes the area including all considered earthquakes (2484) before filtering. Background map: Wikimedia unlabelled map, ©OpenStreetMap contributors and available from [69,70].
A total of 3831 mid-night images represents the analysed data set (Ⴈ), which has been organised on a monthly basis.

Seismic and Tectonic Information
The performed analysis has taken account of seismic events with a magnitude greater than or equal to four (M ≥ 4) recorded in the catalogue of the United States Geological Survey [71]. The catalogue reports 2484 seismic events in the Turkish territory and neighbouring areas from April 2004 to October 2015. The actual area of the earthquake selection is one degree wider than the satellite cut area (see the red line in Figure 2). The distribution of seismic events shows a greater concentration along with the Hellenic Arc, the Isparta Angle, the East, and North Anatolian fault systems.
Since earthquake clusters could affect the results [72], following the approach described in [49] and before the correlation analysis, we filtered the events selecting among those that occurred on the same day within an epicentre-to-epicentre distance of 0.1 degree, the one with the highest magnitude. The resulting "filtered" catalogue, consisting of 2328 events, represents the reference seismic database used in the following analysis.

The Robust Estimator of Thermal InfraRed Anomalies (RETIRA)
Details about the methodology used in this paper and its previous applications are summarized in [17,73]. In this section, we want just to recall that the RST approach rests on the Robust Estimator of TIR Anomalies (RETIRA) index, whose more general expression is the following: ), which has been organised on a monthly basis.

Seismic and Tectonic Information
The performed analysis has taken account of seismic events with a magnitude greater than or equal to four (M ≥ 4) recorded in the catalogue of the United States Geological Survey [71]. The catalogue reports 2484 seismic events in the Turkish territory and neighbouring areas from April 2004 to October 2015. The actual area of the earthquake selection is one degree wider than the satellite cut area (see the red line in Figure 2). The distribution of seismic events shows a greater concentration along with the Hellenic Arc, the Isparta Angle, the East, and North Anatolian fault systems.
Since earthquake clusters could affect the results [72], following the approach described in [49] and before the correlation analysis, we filtered the events selecting among those that occurred on the same day within an epicentre-to-epicentre distance of 0.1 degree, the one with the highest magnitude. The resulting "filtered" catalogue, consisting of 2328 events, represents the reference seismic database used in the following analysis.

The Robust Estimator of Thermal InfraRed Anomalies (RETIRA)
Details about the methodology used in this paper and its previous applications are summarized in [17,73]. In this section, we want just to recall that the RST approach rests on the Robust Estimator of TIR Anomalies (RETIRA) index, whose more general expression is the following: where: • r ≡ (x,y) indicates the geographical coordinates of the satellite pixel centre; • t is the satellite acquisition time with t 2004-October 2014 period. The black line defines the area where the RST approach was applied to identify TIR anomalies on SEVIRI satellite images (SEVIRI cut area). The irregular shape of this area is due to the reprojection of the rectangular (215 rows × 476 columns) SEVIRI image portion from its satellite projection to a latitude-longitude grid. The red line (1 degree wider than the black contoured area) describes the area including all considered earthquakes (2484) before filtering. Background map: Wikimedia unlabelled map, ©OpenStreetMap contributors and available from [69,70].
A total of 3831 mid-night images represents the analysed data set (Ⴈ), which has been organised on a monthly basis.

Seismic and Tectonic Information
The performed analysis has taken account of seismic events with a magnitude greater than or equal to four (M ≥ 4) recorded in the catalogue of the United States Geological Survey [71]. The catalogue reports 2484 seismic events in the Turkish territory and neighbouring areas from April 2004 to October 2015. The actual area of the earthquake selection is one degree wider than the satellite cut area (see the red line in Figure 2). The distribution of seismic events shows a greater concentration along with the Hellenic Arc, the Isparta Angle, the East, and North Anatolian fault systems.
Since earthquake clusters could affect the results [72], following the approach described in [49] and before the correlation analysis, we filtered the events selecting among those that occurred on the same day within an epicentre-to-epicentre distance of 0.1 degree, the one with the highest magnitude. The resulting "filtered" catalogue, consisting of 2328 events, represents the reference seismic database used in the following analysis.

The Robust Estimator of Thermal InfraRed Anomalies (RETIRA)
Details about the methodology used in this paper and its previous applications are summarized in [17,73]. In this section, we want just to recall that the RST approach rests on the Robust Estimator of TIR Anomalies (RETIRA) index, whose more general expression is the following: , being 2004-October 2014 period. The black line defines the area where the RST approach was applied to identify TIR anomalies on SEVIRI satellite images (SEVIRI cut area). The irregular shape of this area is due to the reprojection of the rectangular (215 rows × 476 columns) SEVIRI image portion from its satellite projection to a latitude-longitude grid. The red line (1 degree wider than the black contoured area) describes the area including all considered earthquakes (2484) before filtering. Background map: Wikimedia unlabelled map, ©OpenStreetMap contributors and available from [69,70].
A total of 3831 mid-night images represents the analysed data set (Ⴈ), which has been organised on a monthly basis.

Seismic and Tectonic Information
The performed analysis has taken account of seismic events with a magnitude greater than or equal to four (M ≥ 4) recorded in the catalogue of the United States Geological Survey [71]. The catalogue reports 2484 seismic events in the Turkish territory and neighbouring areas from April 2004 to October 2015. The actual area of the earthquake selection is one degree wider than the satellite cut area (see the red line in Figure 2). The distribution of seismic events shows a greater concentration along with the Hellenic Arc, the Isparta Angle, the East, and North Anatolian fault systems.
Since earthquake clusters could affect the results [72], following the approach described in [49] and before the correlation analysis, we filtered the events selecting among those that occurred on the same day within an epicentre-to-epicentre distance of 0.1 degree, the one with the highest magnitude. The resulting "filtered" catalogue, consisting of 2328 events, represents the reference seismic database used in the following analysis.

The Robust Estimator of Thermal InfraRed Anomalies (RETIRA)
Details about the methodology used in this paper and its previous applications are summarized in [17,73]. In this section, we want just to recall that the RST approach rests on the Robust Estimator of TIR Anomalies (RETIRA) index, whose more general expression is the following: the temporal support [31] identifying the time series of homogeneous (same month of the year, same time of the day) collection of images; • ∆T(r,t) = T(r,t) − T(t) is the difference between the TIR brightness temperature T(r,t) and the spatial average T(t) of T(r,t) on the image at hand. It should be stressed that T(t) computation takes account only of cloud-free pixels, within the investigated region, which are part of the identical category (i.e., only sea or land pixels if r is on the sea or land, respectively); • µ ∆T (r,L) and σ ∆T (r,L) are, respectively, the temporal mean and standard deviation of ∆T(r,t,L) computed on cloud-free pixels belonging to the chosen dataset (t The black line defines the area where the RST approach was applied to identify TIR anomalies on SEVIRI satellite images (SEVIRI cut area). The irregular shape of this area is due to the reprojection of the rectangular (215 rows × 476 columns) SEVIRI image portion from its satellite projection to a latitude-longitude grid. The red line (1 degree wider than the black contoured area) describes the area including all considered earthquakes (2484) before filtering. Background map: Wikimedia unlabelled map, ©OpenStreetMap contributors and available from [69,70].
A total of 3831 mid-night images represents the analysed data set (Ⴈ), which has been organised on a monthly basis.

Seismic and Tectonic Information
The performed analysis has taken account of seismic events with a magnitude greater than or equal to four (M ≥ 4) recorded in the catalogue of the United States Geological Survey [71]. The catalogue reports 2484 seismic events in the Turkish territory and neighbouring areas from April 2004 to October 2015. The actual area of the earthquake selection is one degree wider than the satellite cut area (see the red line in Figure 2). The distribution of seismic events shows a greater concentration along with the Hellenic Arc, the Isparta Angle, the East, and North Anatolian fault systems.
Since earthquake clusters could affect the results [72], following the approach described in [49] and before the correlation analysis, we filtered the events selecting among those that occurred on the same day within an epicentre-to-epicentre distance of 0.1 degree, the one with the highest magnitude. The resulting "filtered" catalogue, consisting of 2328 events, represents the reference seismic database used in the following analysis.

The Robust Estimator of Thermal InfraRed Anomalies (RETIRA)
Details about the methodology used in this paper and its previous applications are summarized in [17,73]. In this section, we want just to recall that the RST approach rests on the Robust Estimator of TIR Anomalies (RETIRA) index, whose more general expression is the following: ). On a monthly basis, we generated two images (µ ∆T and σ ∆T images) used as 'reference images' for the calculation of the RETIRA index. They are representative of expected monthly thermal conditions. To reduce the possible negative impact of the massive presence and/or asymmetric spatial distribution of meteorological clouds on the computation of reference fields and the consequent proliferation of possible false positives (reported, for instance, in [26,35,46]), we adopted here the improved RST pre-processing phases firstly proposed by [49]; • L × L represents the dimension (in pixel units) of the elementary spatial unit centred at location r. L = 1 corresponds to the RETIRA classical configuration (used for Turkey already by [25,56]). For L > 1 (only odd numbers), the variable ∆T(r,t,L) is the spatial mean of the punctual cloud-free ∆T(r',t) values belonging to the L × L pixel box, centred at location r. In all computation phases, the box is considered cloudy when a threshold percentage CT (Cloud Threshold) of cloudy pixels within the L × L pixel box is overcome; • we define Thermal Anomaly (TA) as a (not-cloudy) location where ⊗ ∆T (r,t,L) ≥ K.
The configuration with L > 1 demonstrated (see [11,13] who chose L = 3) to be particularly useful in reducing sporadic and spatially localized TIR anomalies (e.g., due to forest fires, gas-flaring, or accidents related to industrial activities), saving or even highlighting, at the same time, spatially extended TIR anomalous transients, as expected (see next section) in the case of fluctuations related to seismic events. In this paper, K = 3 and, similarly in [11,13], L = 3 and CT = 55% were chosen. The OCA (One-channel Cloudy-radiancedetection Approach, [74]) algorithm was used everywhere to identify cloudy pixels.

Space-Time Persistence Criteria, Significant Thermal Anomalies (STAs) and Significant Sequences of Thermal Anomalies(SSTAs) Definitions
As explained in detail in previous papers [25,26,35,46], natural and observational particular conditions may generate spurious TAs, which have recognisable features (e.g., cloud shadows, coastline shape replication), often appearing isolated in the space (i.e., punctual) and time (i.e., sporadic) domain. Unlike them, the TAs possibly associated with earthquakes are expected (see [17,73] and references therein) to be persistent in the space (e.g., affected area) and time domains (e.g., at least N repetitions within D days). To this aim, specific selection criteria on TAs have been introduced to discard spurious TAs, saving the ones possibly associated with seismic activity and, for this reason, named Significant Thermal Anomalies (STAs). In this paper, in order to isolate STAs, the same criteria adopted in [11,13,49,54] were applied. This can be so summarised in the following: • Identification and removal of spurious TAs due to massive (more than 80% of pixels) presence of clouds on the scene and/or to the so-called 'cold spatial average effect' [26,46]. In both cases, land and sea pixels are separately considered (so that a land/sea TA is excluded if more than 80% of land/sea pixels in the scene are, respectively, cloudy). Similarly, a land/sea TA is removed if the following expression T(t') > µ T -2σ T is verified, being T(t') the spatial average of T(r,t) computed on the cloud-free, land/sea pixels of the image at hand (acquired at time t = t'), µ T and σ T are, respectively, the temporal average and standard deviation of T(t), computed using the homogeneous dataset of images belonging to the temporal domain  Figure 2). In this area, each SEVIRI pixel covers a surface of about 4.5 km (longitude) × 4.7 km (latitude). A total of 3831 mid-night images represents the analysed data set (Ⴈ), which has been organised on a monthly basis.

Seismic and Tectonic Information
The performed analysis has taken account of seismic events with a magnitude greater than or equal to four (M ≥ 4) recorded in the catalogue of the United States Geo- Navigation and co-location errors. TAs close to coastlines clearly replicating their shape [35] are easily identified and removed; 2.
Other known spurious effects. TAs associated with the effects reported in [26,46,49,54] are easily identified and removed.
Moreover, a single STA is considered part of a Significant Sequence of Thermal Anomalies (SSTA) if:

1.
Spatial persistence: it is not spatially isolated being part of a group of STAs (1-degree maximum away from each other) covering an area (affected area) ≥150 km 2 ; 2.
Temporal persistence: the same STAs reappears at least another time in the seven preceding/following days.

Data Analysis and Results
Applying the aforementioned rules to all SEVIRI TIR images from May 2004 to October 2015, we identified 155 SSTAs in the whole time series.
We performed a correlation analysis between SSTA appearances and the location, time, and magnitude in the seismic events of the filtered catalogue. To this aim, we applied empirical rules based on previous studies (see [17,73] and references therein) and physical models (e.g., [20]) proposed to explain TA appearances before and after large earthquakes. In detail, we considered each STA observed at a geographical site (x,y) at the time t and belonging to a specific SSTA, in correlation with an M ≥ 4 seismic event, if the earthquake takes place within the following space-time correlation volume:

•
Temporal window: up to 30 days after (pre-earthquake anomaly) the last or until 15 days before (postseismic/coseismic anomaly) the first appearance of TAs; • Spatial window: within a distance D ≤ R from whatever TAs belonging to the considered SSTA. The distance D is defined under the conditions of 150 km ≤ R ≤ 10 0.43M , the upper limit being the Dobrovolsky radius (in km) [75], corresponding to an earthquake of magnitude M.
The space-time correlation analysis based on the above-mentioned rules highlights that 115 SSTAs (74%) out of 155 could be related to seismic events, while the remaining 40 SSTAs (26%) turn out to be false positives. An example of an earthquake-related SSTA is in Figure 3. More detailed analysis of the earthquake-correlated SSTAs showed that 51 out of 115 (44%) precede the related seismic event, 9 out of 115 (8%) just follow the correlated earthquake (post-seismic SSTAs), and the remaining 55 (48%) appear between two or more earthquakes (pre & post-seismic SSTAs).
In addition, we also analysed the earthquake correlated SSTAs versus the magnitude of related events (Figure 4). It is worth noting that pre-seismic SSTAs always prevail over post-seismic ones, with a decreasing difference as far as the magnitude limit increases. More detailed analysis of the earthquake-correlated SSTAs showed that 51 o 115 (44%) precede the related seismic event, 9 out of 115 (8%) just follow the corre earthquake (post-seismic SSTAs), and the remaining 55 (48%) appear between tw more earthquakes (pre & post-seismic SSTAs).
In addition, we also analysed the earthquake correlated SSTAs versus the m tude of related events (Figure 4). It is worth noting that pre-seismic SSTAs always pr over post-seismic ones, with a decreasing difference as far as the magnitude lim creases. It should be noted that, from an OEF (Operational Earthquake Forecast) poi view, the pre & post-seismic SSTAs are of scarce importance. In fact, their appear between the times of occurrence of different earthquakes cannot be interpreted in rel to future seismic events, but rather to already occurred events. As the goal of this wo It should be noted that, from an OEF (Operational Earthquake Forecast) point of view, the pre & post-seismic SSTAs are of scarce importance. In fact, their appearance between the times of occurrence of different earthquakes cannot be interpreted in relation to future seismic events, but rather to already occurred events. As the goal of this work is still not to provide an OEF, in order to try to solve such an ambiguity, each pre & postseismic SSTA has been re-classified as pre-seismic (post-seismic) if it occurs before (after) the time of occurrence of the earthquake, with the higher magnitude among the ones included in its space-time correlation volume. In case more than one seismic event shares the same magnitude, the closest event in the space domain has been retained. After this re-classification, the number of pre-seismic SSTAs passes from 51 (48%) to 72 (63%), and post-seismic SSTAs passes from 9 (8%) to 43 (37%). Figure 5, reporting the results of such an analysis with reference to earthquake magnitudes, confirms the prevailing appearance of the SSTAs before earthquakes. In particular, although differences poorly change with the magnitudes, the trend of pre-seismic SSTAs seems to increase slightly with the magnitude. However, such a prevalence could be simply due to the fact that the considered pre-seismic temporal window is double the duration (30 days) of the post-seismic one (15 days). So that, after normalizing such a result from this effect, just a small prevailing of post-seismic SSTAs can be registered (from 6% up to 0% moving from the lowest to the highest magnitude class) in comparison with pre-seismic ones.
To better evaluate the significance of the achieved results, the Molchan error diagram ([76] and references therein) was considered. In particular, a customised version of such an approach, fully described in previous papers [11,49], allows us to verify the actual value of the SSTAs as compared with a random alarm function. It plots ν(M), i.e., the percentage of the earthquakes with a magnitude greater than (or equal to) M that are not preceded or followed by SSTAs (i.e., missed earthquakes), versus τ(M), i.e., the percentage of the warned space-time volume considering earthquakes with a magnitude greater than (or equal to) M. These two parameters, ν(M) and τ(M), are computed according to the following equations: (63%), and post-seismic SSTAs passes from 9 (8%) to 43 (37%). Figure 5, reportin results of such an analysis with reference to earthquake magnitudes, confirms the vailing appearance of the SSTAs before earthquakes. In particular, although differe poorly change with the magnitudes, the trend of pre-seismic SSTAs seems to inc slightly with the magnitude. However, such a prevalence could be simply due to th that the considered pre-seismic temporal window is double the duration (30 days) o post-seismic one (15 days). So that, after normalizing such a result from this effect, small prevailing of post-seismic SSTAs can be registered (from 6% up to 0% moving the lowest to the highest magnitude class) in comparison with pre-seismic ones. To better evaluate the significance of the achieved results, the Molchan error gram ( [76] and references therein) was considered. In particular, a customised versi such an approach, fully described in previous papers [11,49], allows us to verify th tual value of the SSTAs as compared with a random alarm function. It plots ν(M), i.e percentage of the earthquakes with a magnitude greater than (or equal to) M that ar preceded or followed by SSTAs (i.e., missed earthquakes), versus τ(M), i.e., the per age of the warned space-time volume considering earthquakes with a magnitude gr than (or equal to) M. These two parameters, ν(M) and τ(M), are computed accordi the following equations: In addition, following Aki [77], the probability gain lines were calculated a 1-Gτ (where G represents the probability gain) together with the confidence curves, based on the null-hypothesis (see [78]), at levels of 95% surrounding the diag of a random guess. Figure 6 shows the result of the error diagram. The events are analysed by m tude classes, referring to all earthquake-related SSTAs ('ALL') or just to the pre-se In addition, following Aki [77], the probability gain lines were calculated as ν = 1-Gτ (where G represents the probability gain) together with the confidence limit curves, based on the null-hypothesis (see [78]), at levels of 95% surrounding the diagonal of a random guess. Figure 6 shows the result of the error diagram. The events are analysed by magnitude classes, referring to all earthquake-related SSTAs ('ALL') or just to the pre-seismic ('PRE') ones. In both cases, the data fall within the optimal strategy zone (left side) with a "significant" deviation from random guessing (i.e., confidence levels of 95%) and, therefore, a non-casual correlation between the SSTAs and earthquakes may be assumed. Looking at the probability gain (Figure 6a), we find values between 1.3 (for M ≥ 4) and 2.2 (for M ≥ 5).
It should be noted that: (a) The non-casual correlation, found in the 'ALL' case, confirms those models (e.g., [20]), suggesting the possible appearance of thermal anomalies both before and/or after seismic events; (b) The non-casual correlation, found in the 'PRE' case, confirms, instead, the predictive capability of the considered parameter; (c) The gain factor seems to be greater for higher (and rarer) magnitude class events, reinforcing the idea that the correlation is driven by physical relations and not just by the high number of events.
Such results, which are in line with the ones achieved by the previous, already quoted, similar studies, also confirm a regional dependence that can be recognized by looking, for instance, at the achieved false-positive rates. Applying similar RST analyses to comparably long (≥10 years) time series over Greece [49], Italy [13], and Japan [11], report quite different false-positive rates (respectively 7%, 40%, 38% to be compared with the 26% of this study) depending on the considered region. The false-positive rate, better than other indicators, can qualify the predictive capabilities of the considered parameter (i.e., thermal anomalies based on RST analyses) within a multi-parametric t-DASH system, even if just in terms of the reliability of the forecast. In fact, the missing rate, another important indicator, would be strongly underestimated in this case due to the possible presence of meteorological clouds, preventing TIR radiation from reaching the satellite sensor.
Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 16 ('PRE') ones. In both cases, the data fall within the optimal strategy zone (left side) with a "significant" deviation from random guessing (i.e., confidence levels of 95%) and, therefore, a non-casual correlation between the SSTAs and earthquakes may be assumed. Looking at the probability gain (Figure 6a), we find values between 1.3 (for M ≥ 4) and 2.2 (for M ≥ 5).
It should be noted that: (a) The non-casual correlation, found in the 'ALL' case, confirms those models (e.g., [20]), suggesting the possible appearance of thermal anomalies both before and/or after seismic events; (b) The non-casual correlation, found in the 'PRE' case, confirms, instead, the predictive capability of the considered parameter; (c) The gain factor seems to be greater for higher (and rarer) magnitude class events, reinforcing the idea that the correlation is driven by physical relations and not just by the high number of events.
Such results, which are in line with the ones achieved by the previous, already quoted, similar studies, also confirm a regional dependence that can be recognized by looking, for instance, at the achieved false-positive rates. Applying similar RST analyses to comparably long (≥10 years) time series over Greece [49], Italy [13], and Japan [11], report quite different false-positive rates (respectively 7%, 40%, 38% to be compared with the 26% of this study) depending on the considered region. The false-positive rate, better than other indicators, can qualify the predictive capabilities of the considered parameter (i.e., thermal anomalies based on RST analyses) within a multi-parametric t-DASH system, even if just in terms of the reliability of the forecast. In fact, the missing rate, another important indicator, would be strongly underestimated in this case due to the possible presence of meteorological clouds, preventing TIR radiation from reaching the satellite sensor.
Quite extended space/time data gaps can be responsible for a number of "missed" events just due to the absence of the satellite measurements that are useful to identify possible SSTAs. To give an idea, 65% of TIR observations used by [11] over Japan were discarded due to the massive presence of clouds; 26% of images were discarded for the same reason in this study. Figure 6. Error Diagram. Full triangles refer to the events (with M belonging to a specific magnitude class) that occurred after SSTAs appearance (pre-seismic SSTAs), identified by RST on the whole studied period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). Empty triangles refer to the seismic events that occurred after or before SSTA appearances (i.e., all SSTA-related events). Each colour indicates a specific magnitude class. (a) Main gain lines are displayed; (b) Coloured lines define the 95% confidence limit curves.
Quite extended space/time data gaps can be responsible for a number of "missed" events just due to the absence of the satellite measurements that are useful to identify possible SSTAs. To give an idea, 65% of TIR observations used by [11] over Japan were discarded due to the massive presence of clouds; 26% of images were discarded for the same reason in this study.
Such variability of results, already reported by [11] for different regions of Japan (apparently in association with the latitudinal variation of the local meteorological regimes), has been investigated for Turkey (for regions that are mostly at the same latitudes) in relation to seismogenic fault distribution.
To this aim, an analysis concerning the SSTA spatial distribution was carried out in order to highlight SSTAs possible concentration in specific areas. As confirmation of the existing link between SSTA appearance and migration of optically active gases, such as CO 2 and CH 4 , alongside seismogenic faults (e.g., [23,79]), almost all (more than 99%) SSTAs were found near fault lines or tectonic plate boundaries, but not in all of them. This could be explained by considering the different seismic and/or degassing activity associated with the different fault lines. Figure 7 shows the SSTAs and corresponding fault systems where they are located. SSTAs over the different fault systems (coloured slices of the inner pie chart) are divided into earthquake-related SSTAs (blue slices of the outer pie chart) and earthquake-not-related (red slices of the outer pie chart) are presented. It is worth noting that: SSTAs were found near fault lines or tectonic plate boundaries, but not in all of them. This could be explained by considering the different seismic and/or degassing activity associated with the different fault lines. Figure 7 shows the SSTAs and corresponding fault systems where they are located. SSTAs over the different fault systems (coloured slices of the inner pie chart) are divided into earthquake-related SSTAs (blue slices of the outer pie chart) and earthquake-not-related (red slices of the outer pie chart) are presented. It is worth noting that:

•
The greatest number of SSTAs are located along with the East Anatolian Fault system (32%), with about 80% of earthquake-related ones (i.e., 20% of false positives); • The Isparta Angle and the North Anatolian Fault system register 15% and 14% of the total identified SSTAs, respectively; 80% of successes (i.e., 20% of false positives) are associated with the former, about 75% of successes (i.e., 25% of false positives) with the latter; Just focusing on the Turkish territory, results on NAF, EAF, IA, and TF were compared with the earthquake and fault densities by [80] (Figure 8), also in agreement with Figures 1  and 2. It is worth noting that the region of the Tuzgölü Fault has a medium-low value of fault density but a very low density of seismic activity. That is in line with the very low percentage of identified SSTAs along with this fault system (2%, corresponding to 3 SSTAs, one of them is earthquake-correlated). On the whole, there is a good correspondence between the areas with a greater number of identified SSTAs (NAF, EAF, IA) and higher values of earthquake density. Such a comparison suggests that SSTAs are not simply related to fault systems (indeed, they appear quite often there), but particularly to those that are seismically active. value of fault density but a very low density of seismic activity. That is in line with the very low percentage of identified SSTAs along with this fault system (2%, corresponding to 3 SSTAs, one of them is earthquake-correlated). On the whole, there is a good correspondence between the areas with a greater number of identified SSTAs (NAF, EAF, IA) and higher values of earthquake density. Such a comparison suggests that SSTAs are not simply related to fault systems (indeed, they appear quite often there), but particularly to those that are seismically active.  [80]); (c) SSTA density map together with the main fault systems (see Figure 1) and earthquakes (green points) considered in this study (see Figure 2).

Conclusions
In this paper, a long-term correlation analysis was performed for the first time over Turkey and neighbouring areas in order to fully characterise the predictive capabilities of the parameter TIR anomaly identified by an RST-RETIRA analysis to be a "candidate" in the framework of a time-Dependent Assessment of Seismic Hazard (t-DASH) system [13,73].
To this aim, a 12-year-long series (about 4000 images) of SEVIRI-MSG images collected in the period of April 2004-October 2015 was analysed using the RST methodology to first identify the TIR anomalies and then SSTAs. A correlation analysis was carried out in order to evaluate the possible correlation of a SSTA appearance with earthquakes (with M ≥ 4) that occurred within a predefined space-time volume in the investigated region/period. In order to evaluate a possible dependence on sub-regional seismic as-

Earthquake density map
Fault density map  [80]); (c) SSTA density map together with the main fault systems (see Figure 1) and earthquakes (green points) considered in this study (see Figure 2).

Conclusions
In this paper, a long-term correlation analysis was performed for the first time over Turkey and neighbouring areas in order to fully characterise the predictive capabilities of the parameter TIR anomaly identified by an RST-RETIRA analysis to be a "candidate" in the framework of a time-Dependent Assessment of Seismic Hazard (t-DASH) system [13,73].
To this aim, a 12-year-long series (about 4000 images) of SEVIRI-MSG images collected in the period of April 2004-October 2015 was analysed using the RST methodology to first identify the TIR anomalies and then SSTAs. A correlation analysis was carried out in order to evaluate the possible correlation of a SSTA appearance with earthquakes (with M ≥ 4) that occurred within a predefined space-time volume in the investigated region/period. In order to evaluate a possible dependence on sub-regional seismic aspects (namely fault density and seismic activity), the analysis was also performed for the sub-regions corresponding to different fault systems.
The analysis allowed us to identify 155 SSTAs on the whole for the period and the whole area of interest. Of which, 74% of them appear earthquake-related (i.e., the falsepositive rate is 26%), occurring within the pre-defined space-time correlation volume around the time and location of the occurrence of M ≥ 4 earthquakes. The random error test performed using a customised Molchan diagram proved that the appearance of SSTAs in spatial/temporal correspondence with seismic events was not casual. Compared with the random guess, we found probability gain values from 1.3 to 2.2 (for M ≥ 4 and M ≥ 5 earthquakes, respectively).
Finally, an additional analysis was carried out considering the SSTA distribution over the main fault systems that affect the Anatolian Provinces and neighbouring areas. It is noticeable that SSTAs seem to appear much more where the seismic density is higher than where the fault density is higher (for instance, along with the East Anatolian Fault system and not in the Tuzgölü Fault area, Figure 8), confirming the dynamic relationship of SSTAs with actual seismicity more than with static tectonic settings. The fact that results are different in the different considered regions (such as in the Hellenic Arc where the false-positive rate reaches the zero value on 19 identified SSTAs) confirms their dependence on specific local/regional tectonic, seismic, and probably also geochemical settings, as shown by the variety of results achieved in previous similar studies [11,13,49,54].
Author Contributions: Methodology, N.G. and V.T.; formal analysis, A.C., R.C. and M.L.; writing, C.F. and V.T.; review and editing, N.P.; supervision, C.F. and V.T. All authors have read and agreed to the published version of the manuscript.