Long-Term RST Analysis of Anomalous TIR Sequences in Relation with Earthquakes Occurred in Greece in the Period 2004–2013

Real-time integration of multi-parametric observations is expected to accelerate the process toward improved, and operationally more effective, systems for time-Dependent Assessment of Seismic Hazard (t-DASH) and earthquake short-term (from days to weeks) forecast. However, a very preliminary step in this direction is the identification of those parameters (chemical, physical, biological, etc.) whose anomalous variations can be, to some extent, associated with the complex process of preparation for major earthquakes. In this paper one of these parameters (the Earth’s emitted radiation in the Thermal InfraRed spectral region) is considered for its possible correlation with M ≥ 4 earthquakes occurred in Greece in between 2004 and 2013. The Robust Satellite Technique (RST) data analysis approach and Robust Estimator of TIR Anomalies (RETIRA) index were used to preliminarily define, and then to identify, significant sequences of TIR anomalies (SSTAs) in 10 years (2004–2013) of daily TIR images acquired by the Spinning Enhanced Visible and Infrared Imager on board the Meteosat Second Generation satellite. Taking into account the physical models proposed for justifying the existence of a correlation among TIR anomalies and earthquake occurrences, specific validation rules (in line with the ones used by the Collaboratory for the Study of Earthquake Predictability—CSEP—Project) have been defined to drive a retrospective correlation analysis process. The analysis shows that more than 93 % of all identified SSTAs occur in the prefixed space–time window around (M ≥ 4) earthquake's time and location of occurrence with a false positive rate smaller than 7 %. Molchan error diagram analysis shows that such a correlation is far to be achievable by chance notwithstanding the huge amount of missed events due to frequent space/time data gaps produced by the presence of clouds over the scene. Achieved results, and particularly the very low rate of false positives registered on a so long testing period, seems already sufficient (at least) to qualify TIR anomalies (identified by RST approach and RETIRA index) among the parameters to be considered in the framework of a multi-parametric approach to t-DASH.


Introduction
A renewed interest on the study of preparatory phases of earthquakes has been solicited in recent years by the, everyday more evident, weakness of traditional approaches to seismic hazard assessment as well as from the significant consequences of their failures in terms of human and economic losses (e.g., WYSS et al. 2012). For instance GELLER (2011) reports that ''…since 1979, earthquakes that caused 10 or more fatalities in Japan actually occurred in places assigned a relatively low probability''. KOSSOBOKOV and NEKRASOVA (2012) measured quantitatively, the errors of the Global Seismic Hazard Project (GSHAP) maps by the difference between observed and expected earthquake intensities. They found that GSHAP accelerations significantly underestimated the observed ones. No better results were achieved in other regions where seismic hazard assessment is based mostly (if not exclusively) on the study of earthquakes catalogs. This is for instance the case of Italy where, on the basis of probabilistic seismic hazard analysis (PSHA) methods (e.g., CORNELL 1968), all the events that caused 10 or more fatalities in the past 20 years were largely unexpected and/or underestimated in terms of magnitude or peak ground acceleration (PGA). For instance the San Giuliano earthquake (31 October 2002, M L = 5.4) that killed 27 kids in a school, occurred in an area that was previously considered of minor concern. Particularly enlightening is the explanation given by CHIARABBA et al. (2005): ''…Seismic hazard for the region had not been previously retained high and the earthquake was mostly unexpected by seismologists. The reason was that neither historical or instrumental events had been previously reported in seismic catalogues for that area''.
In the case of recent Emilia earthquake (20 May 2012, M W = 5.8), the observed PGA ([0.25 g) in the epicentral zone (PANZA et al. 2014) was significantly higher than the one (\0.175 g) predicted by the PSHA map assumed as reference by the Italian Ministry of Infrastructures (see ZUCCOLO et al. 2011 and reference herein) in defining the new rules for building in seismic areas.
Also for these reasons, an everyday increasing interest of scientific community, has been addressed to alternative observational techniques and data analysis methods suitable for improving our present capability to assess seismic hazard in the shortmedium term. In this context a renewed role could be played by the research on earthquake precursors if it is addressed to develop/improve systems for time-Dependent Assessment of Seismic Hazard (t-DASH, TRAMUTOLI et al. 2014) instead to the deterministic earthquake predictions. Several geophysical parameters (see for instance TRONIN 2006 andCICERONE et al. 2009, and reference herein) have been proposed, since decades, as possible earthquake precursors. Although a large scientific documentation exists about the occurrence (in apparent relationship with earthquake preparation phases) of anomalous spacetime transients of their measurements (see for instance TRONIN 2006 and reference herein;CICERONE et al. 2009), no one single measurable parameter and no one observational methodology have demonstrated, until now, to be sufficiently reliable and effective for the implementation of an operational earthquake prediction system (see also GELLER 1997).
A multi-parametric approach seems, instead, to be the most promising approach in order to increase reliability and precision (e.g., HUANG 2011a, b; OU- ZOUNOV et al. 2012;TRAMUTOLI et al. 2012a) of shortterm seismic hazard forecast.
To this aim several research programs have been initiated-like the EU-FP7 project PRE-EARTH- QUAKES (www.pre-earthquakes.org;TRAMUTOLI et al. 2012b) and the Italian INGV-S3 project, https:// sites.google.com/site/ingvdpc2012progettos3/)-or reiterated (e.g., the iSTEP-integrated Search for Taiwan Earthquake Precursors;TSAI et al. 2006) which apply a real-time integration of multi-parametric observations to develop improved systems for t-DASH and earthquake short-term (from days to weeks) forecast.
However, a very preliminary step of whatever multi-parametric approach is to identify those parameters (chemical, physical, biological, etc.) whose anomalous variations can be, to some extent, associated with the complex process of preparation for a big earthquake.
To this aim physical models (e.g., SCHOLZ et al. 1973;TRONIN 1996;FREUND 2007;PULINETS and OU-ZOUNOV 2011;HUANG 2011b;TRAMUTOLI 2013a) that are able to justify the reason why such parameter should/could exhibit significant variations in relation with the preparation phases of an earthquake surely can help but, in some case, just statistical correlation analyses have been proposed in order to establish if or not, and in which measure, such parameter anomalous variations (identified with some specific data analysis technique) can be or not related to an impending earthquake.
In both cases, a convincing demonstration that a not casual relationship exists, among anomalous variations of a candidate precursor and earthquake occurrence, should be provided before deciding to include it in a multi-parametric t-DASH scenario.
Such long-term correlation analyses have been already successfully performed for several parameters like plasma frequency at the ionospheric F2 peak foF2 (e.g., LIU et al. 2006), ionospheric ion density recorded by the detection of electro-magnetic emissions transmitted from earthquake region (DEMETER) satellite (e.g., LI and PARROT 2013), and ultra low frequency (ULF) geomagnetic signal (e.g., HAN et al. 2014;XU et al. 2013). Studies like these represent good exempla of what we can establish as the minimum starting point of whatever multi-parametric observational approach: to identify suitable candidate parameters by measuring their level of correlation with earthquake occurrence on a sufficiently long time series of measurements. This preliminary step will also serve to determine the 286 A. Eleftheriou et al. Pure Appl. Geophys. relative weight to be attributed to each parameter (in comparison with the other contributing parameters) in a multi-parametric system for dynamically updating seismic hazard estimates. However, such long-term correlation analyses, even if necessary, are not always possible due to the lack of long-enough time series of consistent observations. In this context satellite-based measurements could represent a unique reservoir of long-term (more than 30 years in some case), global, continuous dataset.
The fluctuations of Earth's thermally emitted radiation, as measured by satellite sensors operating in the thermal infrared (TIR, e.g., WANG and ZHU 1984;GORNY et al. 1988;QIANG et al. 1991;TRONIN 1996;TRAMUTOLI et al. 2001TRAMUTOLI et al. , 2015aOUZOUNOV and FREUND 2004 and reference herein) as well as in the longwaves (OLR, e.g., OUZOUNOV et al. 2007 and reference herein) spectral range, have been proposed since long time as potential earthquake precursors. However, very refined data analysis techniques are required in order to isolate residual TIR variations, potentially associated with earthquake occurrence, from the normal variability of TIR signal due to other causes (see for instance . More than 10 years (since 2001) of applications of the general Robust Satellite Techniques (RST;TRAMUTOLI 1998TRAMUTOLI , 2007 methodology to this issue, have shown the ability of this approach to discriminate anomalous TIR signals possibly associated with seismic activity (hereafter referred simply as TIR anomalies or TAs) from normal fluctuations of Earth's thermal emission related to other causes (e.g., meteorological) independent of the earthquake occurrences.
Being based on a statistical definition of TIR anomalies and on a suitable method for their identification even in very different local (e.g., related to atmosphere and/or surface) and observational (e.g., related to the time/season or satellite view angles) conditions, RST approach has been widely applied to tens of earthquakes, covering a wide range of magnitudes (from 4.0 to 7.9) and those occurred in very different geo-tectonic contexts (compressive, extensional, and transcurrent) in four different continents.
RST intrinsic exportability permitted its implementation on TIR images acquired by sensors on board of different polar (see DI BELLO et al. 2004;FILIZZOLA et al. 2004;LISI et al. 2010;PERGOLA et al. 2010;TRAMUTOLI et al. 2001) and geostationary satellites (see ALIANO et al. 2007aALIANO et al. , b, 2008aALIANO et al. , b, c, 2009CORRADO et al. 2005;GENZANO et al. 2007GENZANO et al. , 2009aGENZANO et al. , b, 2010GENZANO et al. , 2015. In this paper the level of correlation among TIR anomalies identified by the RST methodology (TRA-MUTOLI et al. 2005) and earthquake occurrence, is evaluated for the first time in a quite long time period (10 years). To this aim 10 years (from May 2004 to December 2013) of TIR satellite records, collected over Greece by the geostationary satellite sensor Meteosat Second Generation-Spinning Enhanced Visible and Infrared Imager (MSG-SEVIRI), have been retrospectively analyzed using the RST approach, and time and location of TIR anomalies occurrences were compared with the ones of earthquakes with M C 4. Achieved results will be discussed in the perspective of a multi-parametric approach for a t-DASH.

RST Methodology
The RST methodology is based on the general approach Robust AVHRR Technique (RAT; TRAMU-TOLI 1998), which, being exclusively based on satellite data at hand (do not require whatever ancillary data), is intrinsically exportable on different satellite packages, the reason why the original name RAT was changed in the more general Robust Satellite Techniques (RST, TRAMUTOLI 2005. The RST approach is based on a multi-temporal analysis of historical data set of satellite observations acquired in similar observational conditions (e.g., same month of the year, same hour of the day, same sensor, etc.).
Such preliminary analysis is devoted to characterize the measured signal (in terms of its expected value and variation range) for each pixel of the satellite image to be processed. In this way spacetime anomalies are identified always by comparison with a preliminarily computed signal behavior. On this basis, TIR anomalies possibly that are related to earthquake occurrences could be defined and isolated from those signal variations which are related to Vol. 173, (2016) Long Term RST Analysis of TIR Anomalies in Greece 287 known (but also unknown) natural and/or observational factors (see ) that could be responsible for false alarm proliferation.
Since the first RST application to the thermal monitoring of earthquake prone areas, TIR fluctuations were identified using the RETIRA index (Robust Estimator of TIR Anomalies, FILIZZOLA et al. 2004;, which can be computed as follows: where: x,y represent the coordinates of the center of the ground resolution cell corresponding to the pixel under consideration on a satellite image; t is the time of the measurement acquisition with t [ s, where s defines the homogeneous domain of multi-annual satellite imagery collected in the same time slot of the day and period (month) of the year; -DT(x,y,t) = T(x,y,t) -T(t) is the value of the difference between the punctual value of TIR brightness temperature T(x,y,t) measured at the location x, y, acquisition time t, and its spatial average T(t) computed on the investigated area considering only cloud-free locations, all belonging to the same, land or sea, class (i.e., considering only sea pixels if x, y is located on the sea and only land pixels if x, y is located on the land). Note that the choice of such a differential variable DT(x,y,t) instead of T(x,y,t) is expected to reduce possible contributions (e.g., occasional warming) due to day-to-day and/or year-to-year climatological changes and/or season time drifts; -l DT (x,y) time average value of DT(x,y,t) at the location x, y computed on cloud-free records belonging to the selected data set (t [ s); -r DT (x,y) standard deviation value of DT(x,y,t) at the location x, y computed on cloud-free records belonging to the selected data set (t [ s).
In this way DT (x,y,t) gives the local excess of the current DT(x,y,t) signal compared with its historical mean value and weighed by its historical variability at the considered location. Both, l DT (x,y) and r DT (x,y) are computed, once and for all, for each location x, y processing several years of historical satellite records acquired in similar observational conditions. They are two reference images describing the normal behavior of the signal and of its variability at each location x, y in observational conditions as similar as possible to the images at hand.
Excess DT(x,y,t) -lDT(x,y) then represents the Signal (S) to be investigated for its possible relation with seismic activity. It is always evaluated by comparison with the corresponding natural/observational Noise (N), represented by r DT (y,x) which describes the overall (local) variability of S including all (natural and observational, known and unknown) sources of its variability as historically observed at the same site in similar observational conditions (sensor, time of day, month, etc.). This way, the relative importance of the measured TIR signal (or the intensity of anomalous TIR transients) can naturally be evaluated in terms of S/N ratio by the RETIRA index. It should be noted that prescriptions on the temporal domain s, which select TIR images acquired in the same period of the year (e.g., month) and in the same hour of the day (e.g., midnight), guarantee the reduction of the signal variability due to the daily (diurnal variation of the temperature) and annual (seasonal variation of the temperature but also of emissivity, which is mainly related to the different vegetation coverage) solar cycles.

A Refined Implementation of RST Approach
In this work all (3151) TIR images acquired by MSG/SEVIRI satellite sensor in the IR10.8 channel (at 9.80-11.80 lm) in the first time slot of the day (00:00700:15 GMT; i.e., 02:00702:15 LT) since May 2004 up to December 2013 have been analyzed. Night-time TIR images were preferred, as usual, because less influenced by effects related to soil-air temperature differences (which are normally higher during other hours of the day) and less sensitive to local variations (due for instance to cloud cover or shadows) of solar illumination which could represent a further element of variability of TIR signal, independent from seismicity. bottom-left 33.5°N-16.7°E) which includes the whole Greek peninsula. Following RST prescriptions, the computation of reference fields was made considering only the cloud-free pixels. In fact thick meteorological clouds are not transparent to the passage of the Earth's emitted TIR radiation so that measured signal in those pixels refers to the cloud top temperature (usually very low) and not to the near surface conditions. Errors in the identification (and the consequent non-exclusion) of cloudy pixels, could heavily condition quality of reference fields which, in general, will result biased toward lower values of averages l DT (x,y) and higher values of standard deviations r DT (x,y). Even if the latter effect-increasing the denominator of expression (1)-could compensate the first one (increasing the numerator of the same expression)-avoiding, thanks to the robustness of RETIRA index, a proliferation of false positives-a significant reduction of the overall sensitivity can be however observed as a consequence of cloudy pixel identification errors.
In order to identify (and discard from reference field computation) cloudy affected pixels, the OCA (One-channel Cloudy-radiance-detection Approach; PIETRAPERTOSA et al. 2001;CUOMO et al. 2004) method (still RST based)-devoted to identify cloudy radiances (i.e., radiances which deviate significantly from the expected values for a specific place and time of observation)-was preferred to traditional cloud detection methods (devoted to identify pixels containing clouds) which are much more exposed to commission (i.e., classifying pixels as cloudy independently if clouds affect or don't affect the measured radiance in the considered spectral band) and omission (mostly because of the use of a fixed threshold approach) errors (see for instance PIETRAP- ERTOSA et al. 2001;TRAMUTOLI et al. 2000).
However, in this paper, due to their importance, particular attention has been paid to cloudy pixel handling, introducing the following refinement of the standard RST pre-processing phases.
(a) In order to be sure that only cloud-free radiances contribute to the computation of reference fields, after cloudy pixels have been identified by OCA, not only those pixels but (like in ENEVA et al. 2008) also the 24 ones in a 5 9 5 box around it (very often belonging to cloud edges) have been excluded by the following computations of reference fields. (b) As shown by ALIANO et al. (2008a) and GENZANO et al. (2009a), spatial distribution of clouds, over a thermally heterogeneous scene, can significantly change the value of the measured signal DT(x,y,t) = T(x,y,t) -T(t) in the remaining (cloud-free) pixels of the scene belonging to the same land/sea class. In fact the same T(x,y,t) value of measured TIR signal can be associated with an higher or lower DT(x,y,t) values depending on the spatial average T(t) computed on the remaining cloud-free pixels (belonging to the same land/sea class of the pixel centered at the x, y coordinates). It has been shown-firstly by ALIANO et al. (2008a) and then by GENZANO et al. (2009a) who named it cold spatial average effect-that, if clouds mostly cover the warmer part of the land (or sea) portion of the scene, the spatial average T(t) will result lower than expected in clear sky conditions. As a consequence anomalously higher values of the signal DT(x,y,t) = T(x,y,t) -T(t) can be measured over the remaining, cloud-free, land (or sea), portion of the scene which are due only to such an anisotropic distribution of clouds along the North-South direction. If not properly taken into account, such a pure meteorological phenomenon not only could (occasionally) introduce false positives in the interpretation phases but will also strongly affect l DT (x,y) and r DT (x,y) reference fields which will be both biased toward higher values with a strong reduction of the overall sensitivity of RETIRA index. To face the problem TIR images suffering from such a cold spatial average effect have been automatically identified and excluded from the computation of reference fields. In order to identify them the values of the temporal average l T of T(t) and its standard deviation r T have been computed (using all the dataset of SEVIRI images collected over Vol. 173, (2016) Long Term RST Analysis of TIR Anomalies in Greece Greece in between May 2004 and December 2013) and for those scenes having T(t) \ l T -2r T [being T(t), l T and r T all computed separately for land or sea pixels] the corresponding land (or sea) portions of image have been excluded from reference field computation. (c) Even if not producing a cold spatial average effect, an extended cloud coverage can determine values of T(t) and then of the considered signal DT(x,y,t) scarcely representative of the actual conditions of cloud-free pixels. So, when the cloudy fraction of land (or sea) portion of the scene was [80 %, then that portion (i.e., all pixels belonging to that land or sea class) have been excluded from the computation of the reference fields l DT (x,y) and r DT (x,y).
On this basis l DT (x,y) and r DT (x,y) reference fields, more reliable than the ones achievable using the standard RST pre-processing phases, have been computed (Fig. 1).
In Fig. 2 it is shown, separately for land and sea classes, the day-by-day results of the analysis performed on the whole data set of TIR images in order to apply tests (b) and (c).
RETIRA indexes DT (x,y,t) have been then computed for all MSG-SEVIRI TIR images belonging to the dataset producing one Thermal Anomaly Map (TAM) D T (x,y,t) for each day t in between May 1st 2004 and December 31st 2013. Locations with DT (x,y,t) C 4-i.e., with signal excess DT(x,y,t)l DT (x,y) C 4r DT (x,y)-will be particularly addressed in this paper and hereafter we will refer to them simply as Thermal Anomalies (TAs).

Identification of Significant Sequences of TIR Anomalies (SSTAs)
As already widely discussed in previous papers (see for example FILIZZOLA et al. 2004;, the RETIRA index, being based on timeaveraged quantities, is intrinsically not protected from the abrupt occurrence of signal outliers related to particular natural (see ALIANO et al. 2008a;GENZANO et al. 2009a) or observational (see FILIZZOLA et al. 2004;ALIANO et al. 2008b) conditions. The particular spatial distribution of this kind of TA and their transitory character in the temporal domain, normally allows to identify them and, in any case, to distinguish them from the spatially and temporally persistent ones possibly related to an impending earthquake, even in the case where they have similar intensity.
This is the reason why (together with relative intensity) spatial extension and persistence in time are requirements to be satisfied in order to preliminarily identify what we call significant thermal anomalies (STAs). Like in all the previous applications to thermal monitoring of earthquake prone area (ALIANO et al. 2007a(ALIANO et al. , b, 2008a(ALIANO et al. , b, c, 2009BONFANTI et al. 2012;CORRADO et al. 2005;DI BELLO et al. 2004;FILIZZOLA et al. 2004;GENZANO et al. 2007GENZANO et al. , 2009aGENZANO et al. , b, 2010GENZANO et al. , 2015LISI et al. 2010;PERGOLA et al. 2010;PULINETS et al. 2007;TRAMUTOLI et al. 2001TRAMUTOLI et al. , 2009TRAMUTOLI et al. , 2012aTRAMUTOLI et al. , b, 2013aTRAMUTOLI et al. , b, 2015b) TA highlighted by RST methodology have been subjected to such a preliminary space-time persistence analysis before it was qualified among STAs.
However, other well-known (see for instance FILIZZOLA et al. 2004;ALIANO et al. 2008a;GENZANO et al. 2009a) spurious effects exist that prevent to include among STAs some, even space-time persistent, sequence of TAs. The ones already identified are: -TAs due to meteorological effects (Fig. 3) These are all anomalous pixels appearing in the TIR scenes affected by a wide cloudy cover or in the TIR scenes affected by an asymmetrical distribution of clouds mainly over the warmest portions of a scene, which expose the remaining clear portions of the scene to the appearance of spurious anomalies (cold spatial average effect, ALIANO et al. 2008a;GENZANO et al. 2009a). Such a circumstance could appear in the portions of TIR scene having the daily spatial average ‹T(t)› B ‹l T › -‹r T › (being ‹l T › and ‹r T › the monthly average and corresponding standard deviation of T(t) computed for the same month of the image at hand using the whole historical dataset of TIR images) or having a cloudy coverage C80 % of total pixels of the same  classes (land/sea). Moreover, also TAs generated by local warming due to night-time cloud passages have been recognized as artifacts of the meteorological effects (see ALIANO et al. 2008a). -TAs due to errors in image navigation/co-location process (Fig. 3) Although, this artifact is not rare for polar platforms, also in the cases of geostationary platforms a wrong navigation may cause intense TAs where sea pixels turn out to be erroneously co-located over land portions (see FILIZZOLA et al. 2004;ALIANO et al. 2008b). -Space-time persistent TAs due to extreme events Usually, these TAs can be observed in relation with particularly rare (over decades) events increasing (for more than 1 day) measured TIR signal because of an increase of surface temperature (e.g., in the case of extremely extended forest fires) or emissivity (e.g., extremely extended floods).
By this way an operational definition of STA can be given by considering a location x, y affected by a STA at the time t if the following requirements are satisfied: (a) Relative intensity DT (x,y,t) [ K (with K = 4 in our case) (b) Control on spurious effects Absence of known sources of spurious TAs (see above) (c) Spatial persistence It is not isolated being part of a group of TAs covering at least 150 km 2 within an area of 1°9 1°( d) Temporal persistence Previous conditions (i.e., the existence of a group of TAs covering at least 150 km 2 within an area of 1°9 1°around x, y) are satisfied at least one more time in the 7 days preceding/following t.
After applying the above-mentioned rules to the whole data set of 3151 SEVIRI scenes over Greece 62 Significant Sequences of Thermal Anomalies (SSTAs) were identified where each one is composed by several STAs spanned on 2 or more TAMs.

Long-Term Correlation Analysis
In order to evaluate the possible correlations existing among the appearance of SSTAs and time, location and magnitude of earthquakes, empirical rules were applied which were mostly based on the long-term (more than 14 years) experience on TAM analyses (ALIANO et al. 2007a(ALIANO et al. , b, 2008a(ALIANO et al. , b, c, 2009BONFANTI et al. 2012;CORRADO et al. 2005;DI BELLO et al. 2004;FILIZZOLA et al. 2004;GENZANO et al. 2007GENZANO et al. , 2009aGENZANO et al. , b, 2010GENZANO et al. , 2015LISI et al. 2010;PERGOLA et al. 2010;PULINETS et al. 2007;TRAMUTOLI et al. 2001TRAMUTOLI et al. , 2009TRAMUTOLI et al. , 2012aTRAMUTOLI et al. , b, 2013aTRAMUTOLI et al. , 2015b) performed by authors in four different continents, different tectonic settings, for tens of earthquakes with magnitudes ranging from 4.0 to 7.9.
By this way each single STA observed at the time t in the location (x,y) will be considered possibly related to seismic activity if: -It belongs to a previously identified SSTA; -An earthquake of M C 4 occurs 30 days after its appearance or within 15 days before 1 (temporal window) -An earthquake with M C 4 occurs within a distance D, from the considered STA, so that 150 km B D B R D being R D = 10 0.43M the DO- BROVOLSKY et al. (1979) distance (spatial window).
By this way, starting from each STA belonging to an SSTA, different possibly affected areas can be built for different possible magnitudes of future/past c Figure 2 Results of the analysis performed to reduce effects related to cloud coverage extension and/or distribution across a scene. Each rhomb represents the spatial average T(t) computed on a TIR image collected over Greece on the day t at 00:00 UTC considering only cloud-free pixels. Blue rhombs correspond to scenes used for the computation of reference fields, the yellow ones to scenes removed from the used data sets. Red lines and light red bands represent, respectively the temporal averages (‹l T ›) and ±2 sigma bounds (2‹r T ›) computed considering the images collected during the same month in the past between May 2004 and December 2013. Vertical gray bars represent the percentage of cloudy pixels identified over each scene. The dashed horizontal blue line indicates the cloudiness limit of 80 % adopted to exclude cloudy scenes from reference field computation (see text). Top only over land; bottom only over sea 1 On the models which foresee the occurrence of similar anomalies also immediately after the quake; see for instance SCHOLZ et al. (1973) and, with reference to TIR anomalies, TRA-MUTOLI et al. (2005, 2013a.

Daily BT spaƟal average (‹T(t)›) of usable images
Monthly BT spaƟal average (‹μ T ›) VariaƟon range of ‹μ T › (±2‹σ T ›) Cloudy coverage (%) Daily BT spaƟal average (‹T(t)›) of discarded images Vol. 173, (2016) Long Term RST Analysis of TIR Anomalies in Greece 293 earthquakes. The convolution of the contours drawn for all the STAs belonging to the same SSTA, allow to draw the contours of the areas (different for different magnitudes) possibly affected by future/past earthquakes. The possible correlation among previously identified SSTAs and earthquake occurrence was investigated considering all earthquakes 2 with magnitude M C 4 occurred from April 1st 2004 to January 31st 2014 in: -The area (contoured in red in Fig. 3)  At first glance it is noticeable that most of STAs respecting the above-described correlation rules appear few days around the time of the earthquake occurrence being generally localized near main tectonic lineaments of the epicentral area.
Looking for instance at the example shown in   In addition to the magnitude of earthquakes, also the temporal gap from the first appearance of STAs is indicated by a number (N) in parentheses (±N means that the earthquake occurred N days after/before the first appearance of STAs). Contours in different colors correspond to different space/magnitude windows (see text). The red contoured box indicates the limits of analyzed SEVIRI TIR scenes (Thermal Anomaly Map area) days after the occurrence of several low magnitude earthquakes (i.e., from 4.2 to 4.5) in the Sea of Crete.
The 4 SSTAs apparently not associated with earthquakes (rows 8, 24, 31, 43) are bounded by a continuous black line 296 A. Eleftheriou et al. Pure Appl. Geophys. reported in Fig. 5, where each row corresponds to one SSTA.
Considering the examples reported in Fig. 4, which are summarized in rows 19 (i.e., Peloponnesus region) and 20 (i.e., Crete island) of Fig. 5, the cells corresponding to the day of appearance of the first STAs (which are part of the same SSTA) are colored in yellow, in that case 26 and 27 June 2007 respectively, and they are located at the day zero on the temporal (horizontal) axis. Days of further appearances of STAs, i.e., 28 June for Peloponnesus region and 29 June for Crete island, in the same area (i.e., persistence) are depicted in red in the same row.
Observational gaps (no data) due to data missing or to overcast conditions are indicated respectively with black and gray cells. Only earthquakes with M C 4 occurred within the prefixed correlation space/magnitude windows (bordered by colored lines in Fig. 4) are reported with numbers indicating their magnitude inside the corresponding cells (days of occurrence).
The analysis performed by applying previously established correlation rules to all the 62 SSTAs identified on the whole time series of SEVIRI TIR observations in the period May 2004-December 2013 highlighted 58 SSTAs (*93 % of the 62 previously identified) in apparent space-time relations with earthquake occurrence and only 4 SSTAs (*7 %) apparently not related to documented seismic activity. Looking at Fig. 6 it is possible moreover to note that: -As foreseen by general (e.g., SCHOLZ et al. 1973) and specific (TRAMUTOLI et al. 2013a) physical models, SSTAs appear mostly before but also after the occurrence of seismic events, showing however an increasing tendency to appear mostly before (more than 66 % of the total) in the cases of medium-high magnitude earthquakes (M C 5); -The presence of meteorological clouds (gray cells in Fig. 6) prevents to guarantee continuity to the observations and to fully appreciate possible space-time persistence of TAs. This is for instance is the case of STAs observed on 2nd and 3rd October 2013 in the Peloponnesus area (sequence 62 in Fig. 5) 10 and 9 days before a large event (M = 6.2) occurred in the Peloponnesus-Cretan Ridge on 12 October 2013 (Fig. 7) after 8 overcast days.

Comparison with a Random Alarm Function: Molchan Diagram
In order to better qualify the possible contribution of the use of SSTAs in the framework of a multiparametric system for a t-DASH, it is important to verify its actual added value in comparison with a random alarm function (see for instance ZECHAR and JORDAN 2008 and reference herein). It is, in fact, particularly important to understand to which extent the very low rate of false positives observed in the considered case of Greece (2004Greece ( -2013 is due to the high prognostic capability of SSTAs or to the very high seismicity of the considered region coupled with an eventually too large used correlation space/time window. To this aim the Molchan approach (MOLCHAN 1990(MOLCHAN , 1991(MOLCHAN , 1997MOLCHAN and KAGAN 1992) has been preferred (even if it usually applies in the absence of observational gaps that, because of clouds, we cannot avoid) to likelihood tests (e.g., ZECHAR and JORDAN 2008) which cannot be applied to models with a non-probabilistic alarm function (SHEBALIN et al. 2014). In our case the Molchan error diagram was implemented (like in SHEBALIN et al. 2006) by plotting the fraction m of missed earthquakes (i.e., apparently non preceded/followed by SSTAs) against the fraction of alerted space-time volume s:  For the computation of m(L) all earthquakes (with magnitude CM) occurred within the investigated area augmented for an area extending up to 1°from its borders, were counted into the denominator of (2). The numerator of (2) reports instead the number of earthquakes (with magnitude CM) not preceded (followed) by an STA in the previous 30 (following 15) days (missed) independently on the actual observation conditions (clouds, data gaps) possibly affecting continuity of observations (see par. 3.3). For this reason computed values of m(L) have to be considered as an upper limit due to the used observational technology which is not able to monitor with continuity the occurrence of all possible STAs.
The numerator of (3) has been computed, for each class of earthquake magnitude, as the sum of all

Figure 7
Example of SSTA identified on 2 and 3 October 2013. As in Fig. 3, TIR anomalies [ DT (r,t) C 4] are depicted in different colors according to RETIRA index values and the colored contours correspond to different space/magnitude correlation windows. Significant TIR anomalies appeared in the Peloponnesus area 10 days before a large seismic event (M = 6.2) happened on 12 October 2013 and before/after some earthquakes with lesser magnitude. As in Fig. 4 the temporal gap of earthquakes from the first appearance of corresponding SSTAs is indicated by a number in parentheses 298 space 9 time volumes alerted by the 62 identified SSTAs: where: t f is the day of last STA appearance plus 30 days, t i is the day of the first STA appearance minus 15 days, -Alerted Area (for each SSTA) is the convolution of all the potentially affected areas (by past/future earthquakes with magnitude CM) corresponding to all STAs belonging to the considered SSTA.
For instance, in the case of the SSTA over Crete island shown in paragraph 3.3 (line 20 in Fig. 5) the potentially affected (alerted) space 9 time volume has been computed by multiplying the total possible affected areas by future/past earthquakes obtained by the convolution of the contours drawn in Fig. 4 for all STAs belonging to that SSTA-for a time period corresponding to the temporal range since 15 days before the first STA appearance (i.e., 27 June 2007) up to 30 days after the last STA appearance (i.e., 29 June 2007) so, in the considered case, the temporal period was 48 days.
To evaluate the alerted space 9 time volume only before the possible occurrence of an earthquake (prevision mode) the numerator of (3) has been computed as follows: where t f is the day of last STA appearance plus 30 days and t 0 is the day of first STA appearance.
In both cases, the denominator of (3) is the total space 9 time volume computed by multiplying the whole area covered by TAMs (contoured in red in Fig. 3) plus 1°in neighboring areas for all the analyzed days (i.e., 3531).

Figure 8
Molchan error diagram analysis computed for different class of magnitude and SSTAs on the whole study period (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). Full circles refer to earthquakes occurred only after the appearances of SSTAs (pre-seismic anomalies). Empty circles refer to earthquakes occurred before or after the appearances of SSTAs (earthquake-related anomalies). Dashed lines indicate points with equal probability gains G (reported as labels) in comparison with random guess results which are represented by the black continuous line Vol. 173, (2016) Long Term RST Analysis of TIR Anomalies in Greece 299 In Fig. 8 are reported the results achieved for different earthquake magnitudes (and spatial correlation window), separately considering the cases of earthquake preceded by SSTAs (full circles) from earthquakes preceded or followed by SSTAs (empty circles). In addition, the probability gain (computed as G = (1 -m)/s, e.g., AKI 1989;MOLCHAN 1991;MCGUIRE et al. 2005) has been also computed and plotted.
Notwithstanding previous considerations on the systematic underestimation of m values (which reflects in a corresponding underestimation of measured gains) from Fig. 8 it is easy to recognize that: • A non-casual correlation actually exists among observed SSTAs and earthquake (M C 4) occurrence with a probability gain (compared with a random guess) in between 1, 8 and 3, 2; • The prognostic value of SSTAs is also non-casual with a probability gain which is in between 1, 5 (for M C 6 earthquakes) and 3, 7 (for M C 5 earthquakes) if only pre-seismic SSTAs are considered; • In both previous cases the maximum gain is achieved when earthquakes with M C 5 are considered.

Conclusions
In this paper the Earth's thermally emitted radiation measured from TIR satellite sensors has been evaluated, on a 10-year long testing period, for its possible relations with earthquake occurrence. 10 years (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) of daily TIR images acquired by the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board the Meteosat Second Generation (MSG) satellite, were used to identify TIR fluctuations in a possible space-time relation with M C 4 earthquakes occurred in Greece in the same period. A refined RST data analysis approach and RETIRA index were used to preliminarily define and then to identify SSTAs. On the base of specific validation rules (in line with the ones used by the Collaboratory for the Study of Earthquake Predictability-CSEP-Project), based on physical models and previous results obtained by applying RST approach to several earthquakes all around the world, the correlation analysis showed that more than 93 % of all identified SSTAs occur in the prefixed space-time window around time and location of occurrence of earthquakes (with M C 4), showing a prevalence for SSTAs (i.e., C66 %) to appear only before the occurrences of earthquakes with M C 5. The overall false positive rate is \7 %. Even if the presence of meteorological clouds in the TIR scenes do not allow to give continuity to the observations (producing, in this way, a possible overestimation of missed events), Molchan error diagram analysis gave a clear indication of non-casualty of such a correlation. A probability gain (compared with a random guess) from 1.5 up to 3.7 was achieved as far as only SSTAs preceding earthquakes with M C 5 are considered. Such a non-casual correlation together with an extraordinary low rate of false positives, both registered on a so long testing period, seems already sufficient to qualify SSTAs among the parameters to be seriously considered (for now at least for Greece) within a multi-parametric system for time-Dependent Assessment of Seismic Hazard (t-DASH, TRAMUTOLI et al. 2014) able to improve short-term (from days to weeks) seismic hazard forecasting.