Using RST approach and EOS-MODIS radiances for monitoring seismically active regions : a study on the 6 April 2009 Abruzzo earthquake

In the last few years, Robust Satellite data analysis Techniques (RST) have been proposed and successfully applied for monitoring major natural and environmental risks. Among the various fields of application, RST analysis has been used as a suitable tool for satellite TIR surveys in seismically active regions, devoted to detect and monitor thermal anomalies possibly related to earthquake occurrence. In this work, RST has been applied, for the first time, to thermal infrared observations collected by MODIS (Moderate Resolution Imaging Spectroradiometer) – the sensor onboard EOS (Earth Observing System) satellites in the case of Abruzzo (Italy) earthquake occurred on 6 April 2009 ( ML∼5.8). First achievements, shown in this work, seem to confirm the sensitivity of the proposed approach in detecting perturbations of the Earths emission thermal field few days before the event. The reliability of such results, based on the analysis of 10 years of MODIS observations, seems to be supported by the results achieved analyzing the same area in similar observation conditions but in seismically unperturbed periods (no earthquakes withML≥5) that will be also presented.


Introduction
Since the eighties, a growing number of studies (see for example Gorny et al., 1988;Qiang and Dian, 1992;Tronin, 1996;Qiang et al., 1997;Tronin et al., 2002;Ouzounov and Freund, 2004) have reported the appearance of space-time anomalies in TIR (Thermal Infra-Red) satellite imagery before (from weeks to days) severe earthquakes.
In order to explain the appearance of anomalously high TIR records near the place and the time of earthquake oc-Correspondence to: N. Pergola (pergola@imaa.cnr.it)currence, different authors attributed their appearance to the increase of green-house gas (such as CO 2 , CH 4 , etc.) emission rates, to the modification of ground water regime (see Hamza, 2001) and/or to the increase of convective heat flux (Qiang et al., 1991;Tronin, 2000).Other, more complex models have been proposed, that include the increase of near surface temperature, among the other expected pre-seismic phenomena (e.g.Pulinets et al. 2002Pulinets et al. , 2006Pulinets et al. , 2007;;Ouzounov and Freund, 2004).
However, such a claimed correlation has been considered, up to now, with some cautions by scientific community.The main problems of the mentioned studies were the lack of a rigorous definition of anomalous TIR signal fluctuations, the absence of a convincing testing step based on a validation/confutation approach and the scarce attention paid to the possibility that other causes (e.g.meteorological) different from seismic activity could be responsible for the observed TIR variations.
Starting from these considerations, a different approach called Robust Satellite Techniques (RST; Tramutoli 2005Tramutoli , 2007) ) -initially named RAT (Robust AVHRR Techniques, Tramutoli, 1998) -was proposed to investigate possible relations between earthquake occurrence and space-time fluctuations of Earth's emitted TIR radiation observed from satellite.
The RST analysis is based on a statistical definition of "TIR anomalies" and an original method for their identification even in very different local (e.g.related to atmosphere and/or surface) and observational (e.g.related to time/season, but also to solar and satellite zenithal angles) conditions; it has been carried out always by using a validation/confutation approach, in order to verify the presence/absence of anomalous space-time TIR transients in the presence/absence of seismic activity.
Published by Copernicus Publications on behalf of the European Geosciences Union.
In this paper, for the first time RST approach is applied to data acquired by MODIS, the optical sensor aboard polar EOS satellites (Terra and Aqua), over Italy.10 years of MODIS data have been processed and RST products analyzed, in particular, in correspondence of the Abruzzo earthquake (6 April 2009, M L ∼5.8).Results obtained by RST at the time of this earthquake will be compared with those obtained by an identical analysis (confutation) performed for the same area in a different, quite sesimically unperturbed (i.e.characterized by the absence of earthquakes of similar magnitude over the same area) year, in order to verify the presence/absence of anomalous space-time TIR transients in presence/absence of significant earthquakes in similar observation conditions.
RST is based on a preliminary multi-temporal analysis of several years (from 4 to 10 depending on the availability of historical data-sets) of homogeneous satellite TIR records, co-located in the space-time domain, devoted to characterize the TIR signal (in terms of its expected value and natural variability range) for each pixel of the satellite image to be processed.On this basis, anomalous TIR patterns are identified as a deviation from those "normal" conditions, using a specific index, RETIRA (Robust Estimator of TIR Anoma-lies, Filizzola et al., 2004;Tramutoli, 2005), to be computed on the image at hand as in Eq. ( 1): where: r ≡ (x,y) represents location coordinates of the pixel centre on a satellite image; t is the time of image acquisition with t∈τ , where τ defines the homogeneous domain of satellite imagery collected in the same time-slot of the day and period (month) of the year; T ≡ (r,t) is the value of the difference between the punctual value of brightness temperature T (r,t) at the location r≡(x,y) and at the acquisition time t and its spatial average T (t) (i.e.T = T (r,t) − 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 r is located on the sea and only land pixels if it is located on the land).Note that the choice of such a differential variable T (r,t) instead of T (r,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; µ T (r) time average value of T (r,t) at the location r ≡ (x,y) computed on cloud free records belonging the selected data set (t∈τ ); σ T (r) standard deviation of T (r,t) at the location r ≡ (x,y) computed on cloud free records belonging the selected data set (t∈τ ).
By this way ⊗ T (r,t) gives the local excess of the current T (r,t) signal compared with its historical mean value and weighted by its historical variability at the considered location.Both, µ T (r) and σ T (r), are computed, once and for all, for each location r, processing several years of historical satellite records acquired in similar observational conditions.They are two reference images describing the normal behaviour of the signal and of its variability at each location r in observational conditions as similar as possible to the ones of the image at hand.The difference T (r,t)−µ T (r) 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 σ T (r) 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.
Being independent form a specific sensor/satellite system, RST can easily be exported on different satellite data, and applied to analyze different events in the world.However, a wider and more detailed description of the RST approach and Nat.Hazards Earth Syst.Sci., 10, 239-249, 2010 www.nat-hazards-earth-syst-sci.net/10/239/2010/As said before, in this work we apply the RST to the thermal data acquired by MODIS, the multispectral (36 bands in the optical range of the electromagnetic spectrum) radiometer aboard EOS-Terra since November 1999 and EOS-Aqua since June 2002.More in detail: we have used the data acquired by MODIS TIR channel 31 (10.780-11.280µm) at 1 km of spatial resolution, so that the T in Eq. ( 1) is the brightness temperature measured by MODIS at these wavelength.A full data-set of 10 years of MODIS observations collected in between 24:00 and 02:00 GMT has been processed, combining records acquired by both (EOS-AQUA and EOS TERRA) satellites orbiting with a delay of about three hours each other.

The Abruzzo earthquake
The Abruzzo earthquake (M L ∼5.8) occurred on 6 April 2009 at 01:32:39 GMT, having its epicentre at 42. 334N and 13.334E (INGV, 2009) in central Italy.The mainshock occurred as a result of normal faulting on a NW-SE oriented structure in the central Apennines.Since January 2009 the zone has been object of frequent seismic events (Fig. 1) with characteristics of seismic swarm with hundreds of modest entity shake.The seismicity was confined in the upper crust and interested an area about 30 km long and strikes in the NW-SE direction, parallel to the Apennine mountain axis and to the main fault structures known in the area.
In order to apply the RST processing scheme, ten years of TIR MODIS imagery, all acquired from 2000 to 2009 at around midnight GMT (from 24:00 to 02:00 GMT) during the months of March and April, were used for computing µ T (r) and σ T (r) reference fields shown in Fig. 2.
On this basis, RETIRA index has been computed for a set of EOS-MODIS imagery in order to perform the validation/confutation analyses.For validation purposes, the months of March and April 2009 have been considered, while, in the confutation analysis, the months of March and April 2008 were selected, being 2008 a seismically "unperturbed" (i.e.no earthquakes with M≥5, in the same region and in the same months but in a different year) one in the whole used data set.

Validation analysis
RETIRA index (⊗ T (r,t)) has been computed, for validation purposes, for all available days from 15 March to 15 April 2009.All scenes at hand have been analyzed and scenes affected by a wide cloudy coverage (>80% of the total scene, e.g. 29 March 2009) are not shown as well as images that could produce TIR anomalies as a consequence of computation of RETIRA index only on a very few number (not statistically representative) of clear pixels (Fig. 3, more details on the so called cold spatial average effect problem can be found in Aliano et al., 2008a).The results are shown in Fig. 4, pixels with ⊗ T (r,t)≥3 (i.e.T (r,t) − µ T (r) excess greater than 3σ  ) are depicted in red and hereafter, for sake of simplicity, we will refer to them as "TIR anomalies".Clouds have been detected by using the OCA approach (Cuomo et al., 2004) and are represented as gray pixels.μ ΔT (r) ) are depicted in red and hereafter, for sake of simplicity, we will refer to them as "TIR anomalies".Clouds have been detected by using the OCA approach (Cuomo et al., 2004) and are represented as gray pixels.) are depicted in red and hereafter, for sake of simplicity, we will refer to them as "TIR anomalies".Clouds have been detected by using the OCA approach (Cuomo et al., 2004) and are represented as gray pixels.for sake of simplicity, we will refer to them simply as "TIR anomalies").As clouds completely mask Earth's emitted radiation in the TIR spectral region, in our analysis clouds affected pixels, identified by using OCA (One-channel Cloudyradiance-detection Approach, Cuomo et al., 2004), have been excluded, as missing data, from whatever further processing and analysis.
Looking at the sequence of pictures of Fig. 4 it is possible to note that, in this case, pixels with ⊗ T (r,t) ≥ 3 appear in the northern part of Italian peninsula on 24 March (in Piedmont region) and on 3 April (in the Padania plain), while affected the Central Italy on 30 March and 1 April, the Calabria's West Coast on 1 April.On 13 April, TIR anomalies are clearly visible in the north part of the Adriatic Sea and Tyrrhenian Sea.TIR anomalies appear in the Balkan region on 30 and 31 March and 1, 3, 4, and 6 April.Other isolated anomalies (1/2 pixels) are present in some of the remaining analyzed days.In order to make identifiable also these cases the scenes presenting even just 1 pixel with RETIRA>3 have been enclosed in a red box in ) signal anomalies are depicted in different colours.Cloudy locations are depicted in gray and the main Italian's tectonics lineament overlayed to scenes are depicted in green.

Confutation Analysis
In order to verify the absence of TIR anomalies in a relatively seismically unperturbed period, the confutation phases has been performed by considering the same period (15 March-15 April) but in a different year (2008).The selection of 2008 for confutation purposes has been done consulting the INGV (2009) seismic catalogue: within the range 2000-2009 no seismic event with magnitude greater than (or equal to) 5 is reported over the investigated area during the months of March and appreciate the time evolution of the TIR anomalies before observed.Low intensity anomalies (Fig. 5) generally follow the stronger ones depicted in Fig. 4, noticeably enlarging the anomaly area and filling gaps both in the space (among isolated anomalous pixels) and time domains.
As already discussed in previous works (Filizzola et al., 2004;Tramutoli et al., 2005), the RETIRA index is intrinsically not protected from the abrupt occurrence of signal outliers related to particular natural (e.g.local warming due to night-time cloud passages) or observational (e.g.errors in image navigation/co-location process) conditions.For this reason, TIR anomalies sequences have been subjected to a space-time persistence analysis in order to discriminate actually significant anomalous spacetime transients from outliers.This is, for instance, the case of the anomalous pixels which appear in Southern Italy on 3 April.A local warming effect originated by the night time passage of a cloudy system and/or by particular industrial emissions coming from the Italsider steel plant (the most important in Italy located in Taranto) have been considered among the possible causes of this spatially extended but not time persistent anomaly.Spatial extension and persistence in time are in fact the further requirements to be satisfied (together with relative intensity) in order to preliminarily identify significant TIR anomalies.
Taking these requirements in mind, four sequences can be discriminated, as also reported in the panels of Fig. 5: -since 30 March -i.e. 7 days before the main shock of Abruzzo earthquake and a few hours before its strongest foreshock (M L ∼4.1) occurred on 30 March at 13:38 UTC -TIR anomalies are visible near seismic epicenters; -in the following days, these anomalous pixels slowly weaken until they disappear on 4 April (areas bordered by black dashed lines in Fig. 5); -anomalous pixels with ⊗ T (r,t) ≥ 2 affect the Padania plain in correspondence the tectonic lineament, for some days beginning from 21 March until 3 April (areas bordered by purple dashed lines in Fig. 5).A seismic event with magnitude M L ∼4.6 occurred in this area (Forlì earthquake) on 5 April 2009 (INGV, 2009); -in the Piedmont region, beginning from 16 March up to 28 March, anomalous values at different levels of intensity are present in correspondence of a tectonic lineament with variable persistence in the space and time (areas bordered by blue dashed lines in Fig. 5).A seismic event with magnitude M L ∼3.9 occurred in this area (Bra earthquake) on 19 April 2009(INGV, 2009); -also in the Balkan region from 22 March to 6 April lower intensity anomalies are present with a variable spatial distribution (areas bordered by brown dashed lines in Fig. 5).The regions has been affected by seismic events with magnitude M L ∼4.2 on 31 March 2009 (INGV, 2009).

Confutation analysis
In order to verify the absence of TIR anomalies in a relatively seismically unperturbed period, the confutation phases has been performed by considering the same period (15 March-15 April) but in a different year (2008).The selection of 2008 for confutation purposes has been done consulting the INGV (2009) seismic catalogue: within the range 2000-2009 no seismic event with magnitude greater than (or equal to) 5 is reported over the investigated area during the months of March and April 2008 (Fig. 6a).As for the validation step, the same cut at 3σ (i.e.⊗ T (r,t) ≥ 3) has been used in order to identify thermal anomalies and scenes affected by a wide (>80%) cloudy coverage are not shown at all (Fig. 6b).
Figure 7 shown the result of confutation analysis.It is possible to note that only isolated (16 March in north Italy, 18 March in northeast of Italy, 20 March in central Italy, 24 March in the Balkan region, 27 March in Sicilia island, 30 March in Tyrrhenian sea, 9 April in southern of Italy, and 12 April in Padania plain) and not time persistent (disappearing just in 1 day) anomalous pixels are detected.Other sporadic anomalies (1-2 pixels) are present in some of the remaining scene marked by red boxes in Fig. 7 to indicate that they contain at least 1 pixel with RETIRA>3.

Conclusions
In this paper, the RST approach has been applied for the first time to thermal image acquired by MODIS sensor on board EOS satellites, in order to verify the possible space-time relationships among TIR anomalies appearance and seismic events occurrence.The Abruzzo earthquake (6 April 2009, M L ∼5.8) has been considered as a test case for validation, while a relatively unperturbed period (no earthquakes with M≥5) was taken for confutation in the same months (March-April) of a different year (2008).Ten years (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) of MODIS TIR observations have been analyzed to characterize the TIR signal behavior in the absence of significant seismic activity.The validation/confutation approach puts in relief that spatially extended and time persistent TIR anomalies (with RETIRA≥3) appear in some space-time correlation with earthquakes of different magnitude occurred in Italy in the considered period (15 March-15 April) and since seven days before the Abruzzo main shock (April 6th 2009, M L ∼5.8).No similar significant (in terms of relative intensity and space-time persistence) TIR anomalies were instead detected during seismically unperturbed periods used for confutation purposes.
The relation between TIR anomalies and earthquake with medium-low magnitude, already found in previous work (e.g.Corrado et al., 2005) appears also confirmed.In fact TIR anomalies persistent in the space-time domain are visible near tectonic feature some days before seismic events of medium-low magnitude (3.9<M L <4.6) occurred in Italy and in the Balkan region in the same period.It should be noted that thermal anomalies of similar relative intensity (S/N≥3) have been never observed by using polar satellites (Aliano et al., 2008a).
Still meteorological clouds seem to limit the potential of such satellite surveys in the TIR spectral regions: -by affecting (reducing usable records) computation of reference fields m(r) and s(r); er to identify anomalous image pixels and scenes affected by a wide cloudy coverage are n (Fig. 6B).Fig. 7 shown the result of confutation analysis.) has been used in order to identify anomalous image pixels and scenes affected by a wide cloudy coverage are not shown (Fig. 6B).Fig. 7 shown the result of confutation analysis.-by generating artefacts (e.g.cold spatial average effect) related to their relative amount and spatial disposition on the scene; -by affecting the space-time persistence analysis reducing the number of identifiable persistences.
The use of existing passive satellite sensors operating in the microwave spectral region could help to overcome such limitations.In fact microwaves, differently from TIR radiation, penetrate (not raining) clouds, allowing us to observe Earth's surface in all weather conditions.Moreover their lower spatial resolution (10-50 km nadir view instead of 1-5 km achievable by TIR sensors) remains largely sufficient to monitor thermal anomalies we always observed at a wider scale around the epicentre zone.

Fig. 1 .
Fig. 1.Seismic events with M L >3.5 occurred in March and April 2009.Red star indicates the main shock of Abruzzo earthquakes (INGV, 2009).

Figure
Figure 2. Reference fields (time average

Figure 3 .
Figure 3. Calculation of cloudy coverage on all the scenes processed for the 2009: grey bars represent the percentage of cloudy pixels over the total number of pixels in the image; red bars represent the percentage of anomalous pixels over

Figure
Figure 2. Reference fields (time average

Figure 3 .Fig. 2 .
Figure 3. Calculation of cloudy coverage on all the scenes processed for the 2009: grey bars represent the percentage of cloudy pixels over the total number of pixels in the image; red bars represent the percentage of anomalous pixels over the remaining cloud free pixels in the image.

Figure 3 .Fig. 3 .
Figure 3. Calculation of cloudy coverage on all the scenes processed for the 2009: grey bars represent the percentage of cloudy pixels over the total number of pixels in the image; red bars represent the percentage of anomalous pixels over the remaining cloud free pixels in the image.

Figure 4 .Fig. 4 .
Figure 4. Validation: results of the RETIRA index computation on the investigated area before and after the Abruzzo earthquake (6 April 2009, Ml 5,8).Pixels with ( ) 3 , ≥ ⊗ Δ t T r are depicted in red.Cloudy locations are depicted in gray.Red boxes identify images with detected TIR anomalies (i.e.pixels with ( ) 3 , ≥ ⊗ Δ t T r Fig. 4. Validation: results of the RETIRA index computation on the investigated area before and after the Abruzzo earthquake (6 April 2009, M L ∼5.8).Pixels with ⊗ T (r,t) ≥ 3 are depicted in red.Cloudy locations are depicted in gray.Red boxes identify images with detected TIR anomalies (i.e.pixels with ⊗ T (r,t) ≥ 3).
6. A) Seismic events (Ml≥4) occurred during March-April 2008.B) Calculation of cloudy coverage on a processed for the 2008: grey bars represent the percentage of cloudy pixels over the total number of pixels red bars represent the percentage of anomalous pixels over the remaining cloud free pixels in the image.ing the figure 7 it is possible to note that only isolated (16 March in north Italy , 18 Marc east of Italy, 20 March in central Italy, 24 March in the Balkan region, 27 March in Si , 30 March in Tyrrhenian sea, 9 April in southern of Italy and 12 April in Padana plain) me persistent (disappearing just in 1 day) anomalous pixels are detected.

Figure 6 .
Figure 6.A) Seismic events (Ml≥4) occurred during March-April 2008.B) Calculation of cloudy coverage on all the scenes processed for the 2008: grey bars represent the percentage of cloudy pixels over the total number of pixels in the image; red bars represent the percentage of anomalous pixels over the remaining cloud free pixels in the image.

Fig. 6 .
Fig. 6.(a) Seismic events (M L ≥4) occurred during March-April 2008 (INGV, 2009).(b) Left: calculation of cloudy coverage on all the scenes processed for 2008: grey bars represent the percentage of cloudy pixels over image; red bars represent the percentage of anomalous pixels over the remaining cloud free pixels in the image; Right: example on 8 April 2008 scene of artefacts (spurious TIR anomalies) due to the poorness of cloud-free pixels where RETIRA index can be computed and to the cloud masking most of the warmest part of the scene (cold spatial average effect, see text).