Toward the development of a multi parametric system for a short-term assessment of the seismic hazard in Italy

With the aim to develop a multi-parametric system devote to improving our present capability to assess seismic hazard particularly in the short-medium term, the preliminary step is to identify those parameters that have demonstrated their non-casual relation with earthquake occurrences and whose anomalous variations can be associated to the complex seismic process. Since 80, fluctuations of Earth thermally emitted radiation, measured by satellite sensors in the Thermal InfraRed (TIR) spectral region, have been reported by several authors in some connection with the occurrence of earthquakes; they and can be considered as one of the possible parameters to include within a multi-parametric system. In this paper, the Robust Satellite Techniques (RST) approach has been exploited to highlight Significant Sequences of TIR Anomalies (SSTAs) possibly related to seismic events happened in Italy in the period June 2004 - December 2014. In particular, we evaluated the level of correlation between occurrence of earthquakes with M≥4 and RST-based TIR anomalies using two different spatial resolutions of MSG-SEVIRI (Meteosat Second Generation -Spinning Enhanced Visible and Infrared Imager) images. Results of the two RST analyses show that the 82% of the identified SSTAs does not disappear when downscaled spatial resolution SEVIRI TIR records were used. For both analyses, more than 60% of SSTAs is apparently in connection with the occurrence of M≥4 earthquakes; more than the 80% of them has tendency to anticipate the occurrence of earthquakes. Although about the 40% of SSTAs is not apparently related to documented seismic activity (false positives), results of random tests (i.e. Molchan’s error diagrams) indicate a non-casual correlation between SSTAs and earthquake occurrences. These results confirm that the parameter “RST-based satellite TIR anomalies” is one of the possible candidates to be included in a multi-parametric system for time-Dependent Assessment of Seismic Hazard (t-DASH).


Introduction
The occurrence of even destructive earthquakes in areas that until then were not considered at high seismic hazard, on the basis of probabilistic maps, highlighted the weakness of approaches such as the Probabilistic Seismic Hazard Analysis [PSHA; Cornell, 1968], widely used for almost 50 years but never validated by objective testing [Mulargia et al., 2017]. Events of Sumatra ( [Geller, 2011;Kossobokov and Nekrasova, 2012;Stein et al., 2012;Kagan and Jackson, 2013] which caused heavy human and economic losses. The increasing awareness that seismic hazard cannot be assessed exclusively through such traditional approaches has progressively led to exploring new ways. Integration of multi-parametric observations seems to be the direction that scientific community is recently taking in the attempt to improve our present capability to assess seismic hazard particularly in the short-medium term. Building up a system based on continuously updated observations (i.e. timedependent), provided that they are adequately selected and analysed, may represent a hopeful research line to actually exploit accumulated knowledge and improve quality of seismic hazard. Isolated attempts of multi-parametric approaches were done during specific conferences (e.g. the 2018 Workshop on Electro-Magnetic Studies of Earthquakes and Volcanoes; EMSEV 2018 website, 2018) or within the framework of international research projects such as EU-FP7
Increased reliability and precision can be expected from multi-parametric t-DASH (time-Dependent Assessment of Seismic Hazard [Tramutoli et al., 2014]) systems. In fact, considering just the intersection of space-time volumes individually alerted by each parameter, a significant reduction of both, the (combined) probability of false alarms and alerted space-time volumes can be achieved. Seismic hazard levels associated to specific space-time volumes can be then systematically updated on the basis of the specific alerted space-time windows and predictive capabilities of the considered parameters. However, the selection of suitable (chemical, physical, biological …) parameters to be integrated within such a t-DASH system is not a trivial task.
In principle, a candidate parameter suitable for integration in a t-DASH system should: a) be selected on the basis of experimental observations and/or convincing physical models (e.g. Scholz et al., [1973]; Tronin, [1996];Freud, [2007a]; Pulinets and Ouzounov, [2011]; Huang, [2011b]; Tramutoli et al., [2013]) that could justify their appearance within the complex process of an earthquake preparation phase; b) be measurable with sufficient space-time continuity; in this context, satellite-based observations may have a great potential because they guarantee information over large areas (up to a global scale) in a quite continuous way. c) exhibit space-time transients (potentially related to EQs) identifiable through clear, scientifically founded and repeatable, data analysis techniques. In fact, even if a huge amount of scientific papers reports the appearance of variations of geophysical parameters (e.g. Tronin, [2006] and reference therein; Cicerone et al., [2009]) in some relation with the occurrence of seismic events, used methodologies and the definition itself of "anomaly" concept are not always clearly stated. d) Show a non-occasional relation between space-time transients and earthquake occurrence, which should be preliminarily established by specific long-term (plurennial) correlation analyses. An actual correlation with earthquake characteristics (e.g. time, position, magnitude) can be appropriately evaluated only if long-term correlation analyses are performed. Only from such analyses it will be possible to quantitatively characterize the predictive capabilities (in terms of successful forecasts, false positive rates and/or missed events) of each parameter with reference to specific regions, alerted space-time volumes, expected earthquake magnitudes, etc.
Good examples of such long-term correlation analyses are represented by studies concerning variations of the plasma frequency at the F2 peak foF2 of ionosphere during 1994-1999 , Ultra Low Frequency (ULF) geomagnetic signals over a 10-year period [Han et al., 2014], ionospheric ion density recorded by DEMETER satellite during more than 6 years [Li and Parrot, 2013] and, more recently, Earth's thermally emitted radiation [Eleftheriou et al., 2016]. e) Provide an appreciable improvement of the EQ forecast capability of the whole t-DASH system.
With all that in mind, in this paper we consider the Earth thermally emitted radiation, measured by satellite sensors in the Thermal InfraRed (TIR) spectral region, as a suitable candidate for a t-DASH system. Since the '80s, several Nicola Genzano et al.
The application of the RST technique to tens of earthquakes, occurred in four different continents and geo-tectonic contexts, with magnitudes ranging from 4.0 to 7.9, using both polar (NOAA-AVHRR, EOS-MODIS) and geostationary (MFG-MVIRI, MSG-SEVIRI, GOES-IMAGER, MTSAT-IMAGER) satellite data, permitted us to find correlations between TIR anomalies and seismic event occurrence, mainly selecting a space window based on the Dobrovolsky et al. [1979] distance and a time window ranging from 30 days before to 15 days after a seismic event. At the same time, such longterm correlation analyses allowed us to highlight the appearance of spurious TIR anomalies, for example, in case of quite rare events which are able to increase (for more than 1 day) measured TIR signal due to an increase of surface temperature (e.g. extended forest fires) or emissivity (e.g. very large floods). However, it is expected that large-scale phenomena (e.g. Earth degassing) continue to produce TIR anomalies also in large size pixels, while small-scale events (e.g. forest fires) can give a decreasing contribution to the TIR signal as the pixel size increases.
In this paper, we evaluate the level of correlation between occurrence of earthquakes with M≥4 and RST-based TIR anomalies using two different spatial resolutions of SEVIRI images. In particular, the same TIR SEVIRI time series (from June 2004 to December 2014) is analysed at full SEVIRI pixel resolution as well as at a downscaled spatial resolution (3 × 3 pixel box). Results will be compared and random tests (Molchan diagram) will be performed in both cases to assess the lack of randomness in the appearance of RST-based TIR anomalies.

RST methodology and RETIRA index
The method proposed by Tramutoli et al. [2001], [2005] , called RST (Robust Satellite Techniques) is mostly based on the general approach named RAT (Robust AVHRR Technique; Tramutoli, [1998]). Being the RAT methodology exclusively based on satellite data at hand (no ancillary data are requested), it is intrinsically exportable on different 4 satellite packages, reason why the original name RAT was changed in the more general RST [Tramutoli, 2005;2007].
The RST approach is devoted to isolating possible anomalies from those signal variations, that are related to known, but also unknown, natural and/or observational factors [e.g. surface spectral emissivity, temporal variations of satellite view angle, etc.; see Tramutoli et al., 2005], which can be responsible for "false alarm" proliferation. It is based on a multi-temporal analysis of a historical dataset of satellite observations, which were acquired in similar observational conditions (e.g. same month of the year, same hour of the day, same sensor, etc.), devoted to characterizing 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, space-time anomalies can be identified by comparing the measured signal, in the geographical position r at time t', with a preliminarily computed signal behavior.
In order to identify abnormal fluctuations of the measured signal (i.e. anomalies), the general index ALICE (Absolutely Llocal Index of Change of the Environment) 1 has been proposed for the first time by Tramutoli [1998] and it is reported below: (1) In order to better characterize the analyzed signal a historical, homogeneous (same platform, period of the day, period of the year, etc.) and long-term dataset of satellite images should be constructed. In fact, selecting satellite images acquired in the same period of the year and in the same hour of the day (e.g. midnight), it is expected a reduction of the signal variability linked to the daily (diurnal variation of the temperature) and annual (seasonal variation of the temperature and emissivity, which is mainly related to the different vegetations coverage) solar cycle. Moreover, in the case of studies on Earth's thermally emitted radiation, using only nighttime images, which are less influenced by effects related to soil-air temperature differences, normally higher during other hours of the day, and less sensitive to local solar exposition, further elements of variability of the signal (independent of the seismicity) can been reduced. In order to identify possible TIR anomalies associated to seismic activity, the RETIRA (Robust Estimator of TIR Anomalies) index has been used. The RETIRA index belongs to the most general class of the ALICE indexes, it has been introduced for the first time by Filizzola et al. [2004] and subsequently has been always applied in all RST applications [Tramutoli et al. , 2015b[Tramutoli et al. , 2018bCorrado et al., 2005;Aliano et al., 2007Aliano et al., , 2008aAliano et al., , 2008bGenzano et al., 2007Genzano et al., , 2009aGenzano et al., , 2009bGenzano et al., , 2015Lisi et al., 2010;Pergola et al., 2010;Eleftheriou et al., 2016] devoted to the thermal monitoring of seismic areas. It can be computed as follows: (2) where: • ≡(x, y) represents location coordinates on a satellite image; • ' is the time of acquisition of the satellite image at hand, with ∊ where defines the homogeneous domain of satellite imagery collected in the same time-slot (hour) of the day and period (month) of the year; • ( , ')= ( , ')-< ( ')> is the difference between the current ( = ') TIR signal ( , ') measured at location , 1 Since Tramutoli [1998] the double l is used to make reference not only to a specific place but also to a specific time '.
and its spatial average < ( ')>, computed in place on the image at hand, discarding cloudy pixels and considering only sea pixels, if is located on the sea, only land pixels, if is located over the land; • ( ) and ( ) are the time average and standard deviation values of ( , ), at location , computed on cloud-free satellite records belonging to the selected homogeneous dataset ( '∊ ).
The choice of such a differential variable ( , ') instead of ( , ') 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.
As proposed by Eleftheriou et al., [2016], a further reduction of signal variability due to particular meteorological conditions, such as the extensive cloudiness and the presence of an asymmetrical distribution of clouds within the analyzed scene, could be achieved selecting only the satellite images that are not affected by these conditions.
In this study, the same criteria adopted in the study of Eleftheriou et al., [2016] have been taken in account. Briefly, land (or sea) portions of the TIR images that are excessively covered by cloud coverage (i.e. ≥ 80% of the total land/sea pixels) or affected by the cold spatial average effect [Aliano et al., 2008a;Genzano et al., 2009a] have been removed from the datasets construction, and consequently excluded from the computation of the RST reference fields (i.e. time average and standard deviation).
Land (or sea) TIR portions possibly affected by the cold spatial average effect are identified on the basis of the following expression ( ') ≤ -2 (where ( ') is the TIR spatial average computed only on cloud-free pixels of the scene acquired at time t' and belonging to the same land/sea class, and are respectively the temporal average and standard deviation of ( '), which are computed using all TIR images of the homogeneous dataset.
In this work, with the aim to reduce the proliferation of TIR anomalies due to sporadic and localized events (for example due to fires or to industrial accidents), and concurrently trying to preserve and emphasize TIR anomalies with an appreciable spatial extension, such as TIR anomalies possibly related to seismic activities, a new index called RETIRA box has been introduced.
For each image, which represents the differential variable ( , '), a spatial filter has been applied. Taking in account the schematic representation of Figure 1, every three pixels starting from the pixel at position row R2 and line L2, the spatial mean value has been computed only if 5 out of 9 pixels within a region of 3 × 3 pixels centered on the considered pixel were not affected by clouds. As a result, a downgrade of the image spatial resolution has been obtained 2 . On this basis, the RETIRA box index can be computed as follows: ( 3) where: • box ( , ') is the spatial average of ( , ') within a region of 3 × 3 pixels centered at location . It is computed only if 5 out of 9 pixels are clear or not close to clouds; • box ( ) is the time average value of box ( , ') at the location computed on cloud-free records belonging to the selected data set ( '∊ ); • box ( ) is the standard deviation value of box ( , ') at the location computed on cloud-free records belonging to the selected data set ( '∊ ).
Note that only cloud-free radiances have been taken in account in the RST elaborations. So, cloudy pixels, which have been identified by OCA (One channel Cloudy radiance detection Approach; Cuomo et al., [2004], and the first neighboring pixels around them, which can be partially affected by clouds (edge clouds) which are not identified by the cloud detection analysis, have been excluded. 2 The pixel size of SEVIRI image over Italian region is averagely 3.5 (longitude) x 4.5 (latitude) km 2 . After the application of the spatial filter the size of new images is ~ 10.5 x 13.5 km 2 .

Significant Sequence of TIR Anomalies (SSTA) and their possible relations with earthquake occurrence
In this study, nighttime (i.e. time slot 00:00-00:15 GMT) images acquired by the sensor SEVIRI (Spinning In addition to Thermal Anomalies (TAs) possibly related to seismic activities, other natural and/or observational conditions could generate a proliferation of TAs in the TAMs, such as night-time cloud passages and errors in image navigation/co-location processes [Aliano et al., 2008a]. In both cases, these TAs have a tendency to vanish with the time and/or to be restricted in some geographical regions (i.e. along the coastlines). This particular character allows us to distinguish them from TAs possibly related to an impending earthquake, which are spatially and temporally persistent . Like the previous RST applications to satellite thermal monitoring of seismic areas, a space-time persistence analysis has been carried out with the aim to qualify TAs highlighted by the RST methodology as Significant Thermal Anomalies [STAs; Eleftheriou et al., 2016;Genzano et al., 2015;Tramutoli et al., 2015b], which could be part of a Significant Sequence of Thermal Anomalies (SSTA).
So, each single location ( ) has been considered (at the time of observation ') to be part of an SSTA, if it Nicola Genzano et al.
6 Figure 1. Schematic representation of the applied spatial filter. Dashed lines represent the original full resolution SEVIRI pixels; instead, continuous lines the downscaled resolution SEVIRI pixels. Every three pixels the spatial mean value has been computed only if 5 out of 9 pixels within a region of 3 × 3 pixels centered on the considered pixel were not affected by clouds.
satisfies the following requirements: • relative intensity: the TA should have a value of the RST-based index 3 ≥ K (in this study K=3.5); • spurious effects are discarded: portions (i.e. land or sea) of the TAMs affected by a wide cloudy coverage (in this study we considered useful only land/sea portions of the scene having a fraction of cloudy pixels less than 80%), navigation errors  and known spurious effects [e.g. cold spatial average effect, in this study identified by this condition ( ) ≤ -2 ; see Aliano et al., 2008a;Genzano et al., 2009aGenzano et al., , 2015 Eleftheriou et al., 2016 for more details] are discarded from the subsequent analyses; • persistence in space-time domain: the existence of a group of STAs which cover at least 150 km 2 within an area of 1×1° around (spatial persistence) and reappear at least one more time in the 7 days preceding/following t (temporal persistence).
The • it belongs to a previously identified SSTA; • an earthquake of M≥4 occurs 30 days after its appearance or within 15 days before (temporal window); • an earthquake with M≥4 occurs within a distance D, from the considered STA, where D=R D being R D =10 0.43M the Dobrovolsky et al., [1979] distance (spatial window). On the basis of the previous discussion, it is assumed D=150km for M≤5. SSTAs have been identified by using the RETIRA index and 28 SSTAs using the RETIRA box index (    North Africa ( 4 (0), 4.5 (+6), 4 (+6), 6.3 (+7), 5 (+7), 4.8 (+7), 4.7 (+7), n.  As announced in the previous paragraph, the possible relation between the highlighted SSTAs and occurred M≥4 earthquakes has been verified on the basis of pre-established rules. The long-term correlation analysis highlights that in both RST analyses more than 60% of SSTAs is in a space-time relation with the occurred earthquakes, with a remain ~40% of SSTAs apparently not related to documented seismic activity. Moreover, looking at Figure 3, it is possible to note that more than the 80% of SSTAs associable to earthquakes has tendency to anticipate the occurrence of earthquakes (more than the 50% of them appears only before earthquakes and not before and after).

Results
With the aim to reveal possible periodicities and characteristics between seismic events and the SSTAs identified by the RST analysis and the RETIRA box index in Italy, Superposed Epoch Analysis (SEA) has been performed. Here, the SEA approach has been implemented by comparing the time of the occurrences of earthquakes associable to SSTAs (i.e. that satisfy the used empirical rules) with the time of the occurrences of First STA (FSTA) appearances.
Here, clusters of earthquakes have been considered as daily single events and the temporal interval between earthquake occurrence and FSTA has been aggregated on the basis of two day counts. Moreover, in order to verify particular tendency of different types of earthquake (i.e. events due to normal faulting or due to reverse faulting), results of the SEA analysis have been differentiated on the basis of earthquake focal mechanisms.  [Dziewonski et al., 1981;Ekström et al., 2012] and the European-Mediterranean RCMT catalogue [Pondrelli et al., 2002] for Balkan and North Africa regions. For the seismic events where moment tensor solution is not computed, the type of faulting has been established in accord with the local tectonic regime of the specific geographical region where earthquake has happened, on the basis of seismic source zone model ZS9 proposed by Meletti et al., [2008]. Results of the SEA analysis shown in Figure 4 highlight that more than 45% of considered earthquakes occurs within a temporal window of two weeks after the first appearances of the significant thermal anomalies and that the 50% of them is due to normal faulting, which corresponds to the 66% of all events due to normal faulting. In order to verify the non-causality of achieved results, an adapted Molchan error diagram approach [Molchan 1990[Molchan , 1991[Molchan , 1997Molchan and Kagan, 1992]  Taking in mind that a proliferation of missed events could occur, due to large and persistent meteorological clouds, as well as the data missing, which prevent from fully appreciating and monitoring with continuity the occurrence of all possible STAs, looking at Figure 5 it is possible to note that a non-casual correlation between SSTAs and earthquakes exists. In fact, by comparison with random guess line (i.e. diagonal of Molchan diagram), both in the case of earthquakes 5 Confidence limit is not computed for M≥6 because only 1 event does not allow us to attribute any statistical significance.

11
To develop a multi parametric system for short-term SHA in Italy that follow the SSTAs (full circles/triangles) and in the case of earthquakes that follow or precede the SSTAs (empty circles/triangles), all considered strategies (i.e. various classes of earthquake magnitude and spatial correlation windows) are in the optimal strategy zone (i.e. below the random guess line). In both cases, achieved statistics do not reject the null-hypothesis of random guessing at the significant level =5% (i.e. the accepted significance level for this kind of study). In fact, in this study it is =14,8% in the case of analysis at full resolution and =17,2% in the case of analysis at downscaled resolution. Finally, although the measurement of the probability gain (G) is not trivial [see Molchan, 1991], an attempt to estimate G has been done following the Aki's paper [Aki, 1989]. Results highlight that a probability gain, which was computed as G=(1-)/ ) up to 4,9 (for M≥5 earthquakes), can be achieved if only pre-seismic SSTAs identified by using the RETIRA index are considered, and up to 6.2 6 for the same classes of earthquake magnitude but considering SSTAs identified by using the RETIRA box index.
In comparison with results achieved by using the traditional RETIRA index, when the RETIRA box is used to identify STAs: • more than 81% of SSTAs survives; • more than 60% of SSTAs is apparently in connection with the occurrence of M≥4 earthquakes; • an improvement in terms of probability gain up to 6.2 as compared to random guess is achieved.

Conclusion
With the aim to define a multi-parametric system able to improve our capabilities to provide short-term (from days to weeks) seismic hazard forecasting, in this paper we have studied the potentiality and the capabilities of the 6 It should be noted that for magnitudes greater than 6 a higher probability gain is achievable, but the little data sample does not allow us any statistical significance to be attributed to it. Nicola Genzano et al.
12 while on the right side the corresponding Molchan error diagram computed to SSTA identified by applying the RST approach on SEVIRI TIR images at downscaled resolution. Full circles and full triangles refer to earthquakes occurred only after the appearances of SSTAs (pre-seismic anomalies). Empty circles and empty triangles refer to earthquakes occurred before or after the appearances of SSTAs (earthquake-related anomalies). Red line and blue line indicate the 95% confidence limit curves surrounding the diagonal of random guess for earthquakes with magnitudes greater than 4 and 5, respectively.
parameter "Earth's thermally emitted radiation measured by satellite sensors in the Thermal InfraRed (TIR) spectral region" to provide useful information about the seismic process.
To this purpose, on the basis of 11 years of MSG/SEVIRI TIR records collected over Italian peninsula, we have evaluated the level of correlation between the occurrence of earthquakes with magnitude greater than 4 and satellite TIR anomalies highlighted by Robust Satellite Techniques (RST). Moreover, for the first time, with the aim to reduce the proliferation of spurious TIR anomalies, due to other causes independent from the seismic activity (e.g. extended forest fires or very large floods) an ad-hoc version of the RETIRA index (i.e. RETIRA box ) was presented.
To evaluate its accuracy, the same TIR SEVIRI time series (from June 2004 to December 2014) has been analysed computing both the RETIRA and the RETIRA box index. Results of the two RST analyses have shown that the 82% (i.e. 27 of 33) of the identified SSTAs does not disappear when a downscaled spatial resolution SEVIRI TIR records was used. These circumstances confirm that the low spatial resolution is not a constraint to monitor possible preseismic thermal anomalies; in fact, large-scale phenomena (e.g. Earth degassing) continue to produce TIR anomalies also in larger size pixels.
Results of correlation analyses performed applying pre-established rules highlight that in both the RST analyses: • more than 60% of SSTAs is in a space-time relation with the occurred earthquakes, with a remaining ~40% of SSTAs apparently not related to documented seismic activity (false positives); • more than the 80% of SSTAs associable to earthquakes has tendency to anticipate the occurrence of earthquakes, with a 50% of them appearing only before earthquakes and not before and after the seismic events; • superposed epoch analysis indicates that more than 45% of earthquakes associable to SSTAs occurs within a temporal window of two weeks after the first appearance of the significant thermal anomalies.
Taking in mind the intrinsic technological limitation represented by the presence of meteorological clouds, which prevent from appreciating space-time continuity of observed thermal anomalies and consequently which can produce an overestimation of missed events, the results of Molchan's error diagrams show a precursory correlation in both performed analyses. Moreover, moving from the analysis of TIR SEVIRI time series at full pixel resolution (i.e. SSTAs identified by using the RETIRA index) to the analysis at a downscaled spatial resolution (i.e. SSTAs identified by using the RETIRA box index), an increase of probability gain has been obtained.
Finally, results achieved in this work confirm that satellite TIR anomalies highlighted by the RST approach and the RETIRA-based indexes provide useful information on seismic processes and oncoming earthquakes. Although a long-term correlation analysis shows an intermediate number of false positives, on the basis of random test analysis (i.e. Molchan diagram) we can conclude that the parameter "RST-based satellite TIR anomalies" can be included in the framework of multi-parameter systems devoted to short-term earthquake forecast [Tramutoli et al., 2014], which should be based on well-founded observational methodologies and data analysis techniques. So, thanks to the systematic integration of the short-term observations with the traditional approaches (e.g. based on seismological catalogues) devoted to defining a probabilistic assessment of the seismic hazard in the long-term, it could be expected an increase in terms of reliability and precision of seismic hazard assessment, which can be systematically updated on the basis of the specific alerted space-time windows and predictive capabilities of the considered parameters.