From high temporal resolution to synthetically enhanced radiometric resolution: insights from Night Thermal Gradient results

Electromagnetic emissions in thermal infrared bands are an important research topic on pre-earthquake studies. Satellite thermal data have been investigated by many independent research groups looking for their anomalous behaviour before main earthquake's occurrences. Among them, geosynchronous satellite data are reported as less prone to artefacts during data processing. In this work, the Night Thermal Gradient (NTG) algorithm is presented, which has been specifically proposed for geostationary thermal infrared data processing. NTG method relies on the exploitation of high temporal resolution data to find coherent low frequency components of a hypothetical precursory signal of seismic activity. In this paper, the method is presented by giving details about the applied procedures, steps, theoretical assumptions and results obtained during the studies of L'Aquila 2009 earthquake and the seismic activity of Central Italy and Sardinia.

Abstract. Electromagnetic emissions in thermal infrared bands are an important research topic on pre-earthquake studies. Satellite thermal data have been investigated by many independent research groups looking for their anomalous behaviour before main earthquake's occurrences. Among them, geosynchronous satellite data are reported as less prone to artefacts during data processing. In this work, the Night Thermal Gradient (NTG) algorithm is presented, which has been specifically proposed for geostationary thermal infrared data processing. NTG method relies on the exploitation of high temporal resolution data to find coherent low frequency components of a hypothetical precursory signal of seismic activity. In this paper, the method is presented by giving details about the applied procedures, steps, theoretical assumptions and results obtained during the studies of L'Aquila 2009 earthquake and the seismic activity of Central Italy and Sardinia.
Data collected by geostationary satellites have proved to be often the most effective [20] and allowed researchers to propose and develop specific algorithms for data processing by temporal interpolation over specific amounts of samples [9,15,17]. Nevin Bryant [9,15] proposed a first interpolation algorithm, showing the possibility to use more than one sample per night to highlight night anomalous behaviours of the observed thermal images pixels.
Piroddi [17,21] introduced an algorithm named Night Thermal Gradient (NTG) that allows the detection of anomalous behaviours by means of the interpolation of pre-processed geostationary data generated from measurements collected during several consecutive nights. The NTG method was applied to the L'Aquila's earthquake a e-mail: lucapiroddi@yahoo.it

112
The European Physical Journal Special Topics in Italy (M L 5.9 M W 6. 3). This test demonstrated the good time-accordance with previously published thermal infrared (TIR) studies [12] and allowed subsequent geospatial interpretation by different research groups [22][23][24].
Differently from former studies, in this paper, the NTG method is presented in detail, and its potential strengths and weaknesses are discussed in depth. In addition, cases of application of the method on different geographic areas are presented in detail: L'Aquila, Central Italy, and Sardinia.
The paper is organized as follows: -Section 2: the NTG method is summarized by, first, looking at the processing steps and, second, highlighting theoretical assumptions and consequences on final results of the processing choices; -Section 3: the application context is introduced presenting the case study pecu- 2 Night Thermal Gradient method 2

.1 The algorithm
With the assumption that a thermal signal can be a precursor of seismic activity, the Night Thermal Gradient (NTG) method consists in the application of a workflow to enhance the Signal-to-Noise Ratio of a thermal signal and reduce the effects of weather conditions (cloud movement, but, also, cloudy/clear day or daily temperature dynamics) [17,21]. Since clouds or anomalous weather conditions are not stable over time and homogeneous in space, their effects on single pixels signal can be reduced through time-lapse measurements.
As shown in Figure 1, the NTG approach for the weather effect attenuation can be divided into the following three steps: -1st step: calculation of the mean temperatures for specific instant (sample) of the day, by means of a moving window lasting several days; -2nd step: extraction of the overnight averaged data; -3rd step: extraction of the temperature derivative and intercept from the selected data obtained in the 2nd step. Figures 2, 3 and 5 illustrate how the processing workflow is implemented on the temperature timeseries for each pixel of the thermal images. The individual plots in these three Figures show the mean daily temperature-versus-sample for a specific pixel (the mean temperature is calculated over N thermal images, i.e. N days). So, in short, each curve represents the temperature variation on a specific pixel along a day of acquisition.
The main processing steps are applied as follows: -1st step: on raw data, and similarly to many geophysical methods, a stacking procedure is applied averaging the contributions to the same pixel, at the same time, over a time-window of N days (Fig. 2). A family of raw data is selected from the analysed day backwards for a stack-window of N days (blue curves) and the At the 1st step of the processing, stack-temperature is obtained (T N ) averaging the T values over the considered N -day window (N is a constant value; in this case, N = 9 was selected); this new average temperature is the 4D function of the newest date involved in stack processing (DAY stack ), position and sample; at this step, T N data per (stack-)day are still 96 as in the initial dataset (clearly, the specific T N values depend, in any case, on the choice of the window duration N ). At the 2nd step, only stack-temperature T N data occurring during the night-time are selected (41 samples in the specific cases in this paper). At 3rd step, after linear regression of T N data, two indices are obtained from each stack-night (41 samples): slope (dT /dt or NTG index) and intercept temperature; they are 3D functions of position and dates involved in stacking.
average of the temperature values is calculated for the specific phase condition. The result is the stack-day (red curves). This procedure is repeated on the entire dataset, with a constant stack-window (in a sliding window average fashion). For the analyses presented in this paper, the stack window has been chosen 9 days long. This stack-window width is the result of several attempts aiming at optimizing the robustness of the retrieved temperature without (over-)smoothing important features.
-2nd step: analysing the curves obtained from different measurements. Day-time data were more affected by random environmental noise than night-time ones. As a matter of fact, during the solar energizing phase, the effects of the variable cloud cover are temporally amplified with respect to the effects that the same cloud cover has during the night. In addition, real temperature curves (raw data) are more often not continuous and interrupted by lack of data due to bad meteorological conditions. Consequently, night-time samples were preferred, being a more stable signal with easier anomaly recognition. Thus, two consecutive stack-days data are joined into a consistent stack-night (Fig. 3).
-3rd step: the thermal signal is processed with a least squares regression in order to extract few representative indices of the whole night. In this step, a preliminary analysis of the data-fitting algorithms was performed to choose the polynomial 114 The European Physical Journal Special Topics Fig. 2. NTG processing flow: the stacking processing. Curves, here and in the figures above, report the stacked data of some pixels of analysed dataset: they are used for a symbolic representation of daily curves for both raw and processed data with colour as their processing condition classification. One thicker line is used to evidence the processing operations. function better fitting the sample dataset. To do this, different tests have been performed, looking at the Root Mean Square Error (RMSE) as an index of the fitting representativeness (Fig. 4). Since the results were substantially equivalent for all the fitting functions, the linear regression was preferred because of the minor number of indices resulting after its application. In addition, a robustness constraint [25] was tested, leading to clustering of constrained and non-constrained functions  into two families. Linear robust interpolation was the final choice because it was less affected by outliers. By choosing the linear fitting, two indices are extracted for each stack-night, the slope (dT /dt) and the intercept temperature (T 0 ). The slope (i.e. the Night Thermal Gradient index -NTG index) was considered the most indicative of the actual system dynamics and the least affected by external factors like the varying diurnal energisation (Fig. 5).
A pixel behaviour is considered anomalous if NTG shows a positive slope (an increasing temperature is strange if occurs at night). After processing of the experimental dataset, a unique colour scale range was selected for all maps (±15 • C/minutes * 1/1500), and a pixel was considered anomalous when assumes values near to saturation (about 60% of NTG positive value saturation, 80% of the full values-scale) or upper [21].

NTG theoretical assumptions: index stabilization and expected resolution increase
From a theoretical point of view, the application of the stacking (step 1) can be more effective if the searched hypothetical thermal signal (anomaly) has a wavelength larger than the duration of the filter span. In fact, in this case, we can assume that random noise is strongly reduced by the averaging procedure, while the coherent component of the signal is less attenuated and made recognizable (Fig. 6a). After the application of the stack operator, the resulting maxima of the thermal functions are, in addition, forward-shifted (Fig. 6a) reducing the time for possible alert/precursory. As shown in Figure 6b, the derivative nature of NTG operator moves back the alert time because the positive NTG values anticipate the thermal curve maximum (which corresponds to zero value of the NTG index), increasing the time for possible alerts.
The stack-window for L'Aquila case was dimensioned on a few months dataset with the attempt to optimize signal stability in terms of space and time persistency, reducing high frequency anomalies, but keeping visible the feature preceding main shock. The results obtained with this procedure are consistent with the results obtained from other studies [12].
The extraction of the information from raw data, obtained with the application NTG algorithm to exploit the high temporal resolution of geostationary satellites, is expected to act as a synthetical enhancement of radiometric resolution. In fact, the discussed approach, consisting in up to 41 samples repeated every night for 9 days (in hypothetical clear sky conditions), allowed obtaining two indices representative of the studied thermal signal, the NTG index and the intercept temperature. In this way, together with the change of datasets formats from 16 bits integer to 64 bits floating numbers, which reduces the losses due to numerical approximations, experimental noise is minimized, and coherent information, in space and time, emphasised. It is known that in remote sensing three types of resolutions (spatial, spectral, radiometric) are closely related [26]. The fourth, temporal resolution, is often considered transversal to the other three. Classical definition of radiometric resolution refers to the ability of the sensor to differentiate among subtle variations in brightness that is the energy measured by the instrument [27]. With the application of the NTG algorithm, we synthetically obtain on final data the same effect of an increase of radiometric resolution on the sensor, provided that the models of data merging (stack and linear regression) are representative of the studied phenomena. On the other side, the application of NTG algorithm can produce temporal-frequency changes of the original signal components contained in the raw data. In addition, the two merging steps above described (stacking and linear regression) can produce artefacts on the final maps due to the time distribution of the weather anomalies.
NTG index (dT /dt) has been designed to potentially be used together with other existing and new thermal parameters for the study of land surface temperature dynamics related (but not limited) to earthquakes. In fact, enhancing radiometric resolution leads to apparently improved spatial resolution and stability of thermal patterns for low frequency thermal signals [17,21,28]. The NTG results should also be compared with raw data to reduce artefacts risk. As reported in [21], the anomaly recognition (dT /dt > 0) could also be done applying less physical (i.e. losing physical units in final indices) and more statistical strategies such as frequency analyses [29][30][31][32], standard deviation normalisation [9,15,33] or other robust techniques [20,34].
3 The application context (experimental data and seismic activity) 3.1 First experimental setup and data NTG index has been applied for L'Aquila earthquake monitoring. For the calculation of this index geostationary satellite data have been used starting from Land Surface Temperature (LST) data estimated and provided by the agency Satellite Application Facility on Land Surface Analysis (LSA SAF) of the EUMETSAT consortium. LST data were based on the acquisitions in the TIR channels of the SEVIRI sensor on a) b) Fig. 6. Theoretical assumptions of NTG operator. Effect of the stacking step on (high frequency) random noise and a supposed (low frequency) coherent component of the thermal signal (a): the highest frequency components are more attenuated than the lowest ones, the latter being readable at the stacked signal. Setting up the stack filter to work from a day to N days before, in the stacked signal the effect of the hypothetical coherent anomaly lasts N days more than in the original dataset and its maximum is forward-shifted of N/2 days. Decomposition of the expected stacked thermal signal, land surface temperature (LST) data of the observed scenario, in two regular components as they can be assumed after the work of an optimal stacking processing (b): a variable regional field given by the repetition of the diurnal cycle, and a transient local field given the supposed thermal signal generated by pre-earthquake dynamics. As it is possible to observe from the central plot, the fact that the slope angle of the line tangent to the signal is zero at the maximum of the curve, means that the definition of positive slopes as anomalous implies the recognition of the anomaly before its maximum intensity. Meteosat9. LSA SAF algorithm used for the estimation of LST is based on the split window method applied to the bands of 10.8 and 12.0 µm [35]. LST raw data have been used at their full sampling rate (one data every 15 min). The spatial resolution is variable over pixels of the scene observed (a mean value could be around 6 km). In order to compare the index behaviour in presence and in absence of relevant seismic activity, two regions (close each other) have been selected: Sardinia and Central Italy (Fig. 7a). The first is characterised by a stable absence of seismic activity and is sufficiently wide to capture expected extensions of thermal phenomena (55 × 65 pixels). The second is characterised by the presence of seismic activity and includes the maximum coverage of Apennines mountain chain (150 × 150 pixels). Furthermore, Central Italy has been divided into two regions based on the distance from L'Aquila epicentre as the major seismic activity was clustered there: looking at seismic catalogues, the most active region has been chosen with latitude between 42 • and 43 • N , while external regions were in intervals 39 • -42 • and 43 • -45 • N (Fig. 7b). The proximity of the investigated areas is crucial for the validation of the methodology due to similar meteorological, climatic and observational (i.e. satellite angle of view and atmospheric layers thickness) conditions. These regions were observed, and related data were processed for two distinct periods selected based on the seismic activity. For both regions, the first period, characterised by only background seismic activity, is from October 1st to 31st, 2008; while the second period, characterised by major seismic activity, goes from January 27th to May 8th, 2009 (see Tab. 3). Data collected on Central Italy in the first period (October 2008) show no strong variations in the seismic activity; consequently, the two areas shown in Figure 7b have been merged (see Tab. 3) for October 2008 analysis.

Seismic context
L'Aquila main shock happened on 6th of April 2009 at 3:32 local time (1:32 UTC). It was classified as M W 6.3 and M L 5.9 with a depth of 8.8 km by the Istituto Nazionale Table 1. Seismic sequence dynamics of the area (50 km radius) surrounding L'Aquila epicentre before the main shock, through Gutenberg-Richter law slope and earthquakes daily frequency (after [39] [36,37]. From October 2008 to the mainshock, local seismic activity gradually increased [38,39]. As summarized in Table 1 [39], the system (chosen in a radius of 50 km from the epicentre) remained in background activity condition form January 2006 until the end of October 2008. After that, weak foreshock activity was observed until about ten days before the mainshock. The data reported in Table 1 were obtained using a catalogue with earthquakes with M L > 1.3 and by looking at their rate (events per day) and the b-value (dimensionless), slope of the well-known Gutenberg-Richter or magnitude-frequency law [40]: where N M is the number of recorded earthquakes with a magnitude larger than a threshold M . The slope, b, should be around one for background conditions and drop as soon as foreshock activity is becoming mature, being significantly lower in foreshocks than in aftershocks/background seismicity conditions [39]. For the indicated time windows and magnitudes, seismicity rate gradually grows up while the b-value drops down. System acceleration ten days before L'Aquila earthquake [39] was confirmed by the analyses based on an Accelerated Strain Release approach [41].
The seismic activity in the last ten days before the mainshock had its largest event on the 30th March (M L 4.1) and led to the main shock in the night between 5th and 6th of April. The following seismic activity was further clustered and decayed as usual aftershock activity [36,37]. Table 2 shows that the major aftershock activity (M L > 4) was concentrated in the seven days later the mainshock, and in particular the two biggest aftershocks respectively, occurred one (M W 5.6, M L 5.4) and three days (M W 5.4, M L 5.1) after the mainshock.
The whole seismic activity was mostly concentrated in the areas immediately surrounding the mainshock epicentre. For a map visualization of the main seismic sequence (foreshocks and aftershocks), see for instance [38]. Nevertheless, the day before the L'Aquila mainshock (at 20:20 UTC), a significant event (M L 4.6, depth 28.2 km) occurred outside the highest active area, toward north, near Forlì (latitude 44.236 and longitude 11.999): this event could also be linked to thermal transients as discussed in the followings.

Introduction to NTG results
The results obtained from the parallel analysis of the two areas are summarized in Table 3, where the rows identify the periods and the areas of observation, whereas   columns indicate the number of earthquakes with magnitude larger than 3 and the number of thermal anomalies (NTG [dT /dt] anomalies). In particular, the NTG anomalies were defined as several adjacent pixels characterized by values close to colour-scale saturation. In addition, NTG anomalies were classified: -As "case-1" when occurring up to 15 days before the earthquake and within a 60 km radius from the epicentre; -As "case-2" if occurring without a close (in space and time) earthquake event.
If an anomalous thermal pattern is present for more days, it is considered as a unique NTG anomaly.
In Central Italy, as shown in Table 3, during October 2008 only 1 NTG anomaly is reported, but this seems not to be connected to the unique earthquake larger than the magnitude threshold detected in that period. Paragraph 4.2.1 and Figure 9 elaborate on this anomaly (case-1 condition).
In Spring 2009 (from 27 January to 8 May), for the area in Central Italy, but far from the most active area, 11 NTG anomalies are observed, 6 of which follow the space-time constraints (case-1) and 5 don't (case-2). One of the 6 NTG anomalies of the case-1 type is classified as the potential precursor of Forlì's seismic event happened on April 5th (paragraph 4.2.2).
In the same period for the area close to the mainshock region, a completely different pattern can be seen: characterised by 243 seismic events (M L > 3) but only 2 NTG anomalies, one of which is potentially linked with the 6 of April main shock. These two anomalies are classified as case-1 and are explained in paragraph 4.2.3.
Concerning the Sardinia region, in October 2008, none NTG anomaly was reported, while in Spring 2009 there were 2 NTG anomalies coupled to no seismic activity (case-2).
Regarding the observation of the most active seismic condition, it is important to highlight that, as introduced in Figure 7b and paragraph 3.2, the area of Central Italy had much more seismic activity than other areas for the same period. The related seismic events are often interconnected and dependent on each other: because of the preliminary nature of the statistical analysis presented, in this paper, seismic events were not declustered and are always considered if they are bigger than the threshold M L value of 3. Looking at the two Spring 2009 Central Italy data datasets it was possible to empirically evidence the different sensitivity, of highly active regions and less active ones, to thermal anomalies' appearance (in fact moving from lower to higher active regions, the NTG anomalies, typically but not generally, trend to be spatially and temporally associable to higher magnitude events). Moreover, the classifications of case-1 and case-2 conditions, considering relatively small datasets, have the main objective of summarizing the results presented in next paragraphs and helping in distinguishing different space-time relationships of observed NTG anomalies with raw seismic data. For any statistical analysis of precursory correlations, seismic events should be isolated as independent and seismic catalogue should be declustered, but this was not the purpose of the present article and the procedure was not performed.
In the following, the results are presented in several maps.

Background analyses
Firstly, for the two study areas, the periods in which there was no relevant seismic activity are considered in order to compare the seismic background activity with the measurements of thermal anomalies highlighted by the NTG algorithm applied to the raw LST data provided by EUMETSAT consortium. The background analyses were conducted on Central Italy (October 2008) and Sardinia (October 2008 and Spring 2009). Some exemplifications of NTG dynamics over Sardinia region are presented in Figure 8.
The top-left and bottom-right maps show a common condition of negative NTG values for the stack-nights between 14 and 15 October 2008 (that means raw data collected from 6 October nightfall to the night between 14 and 15 October) and between 21 and 22 of March 2009 (that means raw data collected from 12 October nightfall to the night between 21 and 22 October), respectively. The top-right map (28)(29) October stack-night) shows some very weakly positive NTG index values, but very far from saturation and consequently not classified as anomalies as they remain below the 60% of maximum positive colour scale. Finally, the bottom-left map corresponding to the 8-9 March stack-night shows one of the two patterns recognized as anomalous and classified as case-2.
The NTG monitoring results for Central Italy during October 2008 are quite similar to those reported for Sardinia in Figure 8 up-left and 8 bottom-right where NTG values are generally negative [17]. The NTG dynamics of Central Italy for the period 27 to 30 October are shown in Figure 9, where the time-series highlight an intense anomaly over the central region of the maps. As introduced in paragraph 4.1, this anomaly is classified as case-2 (in agreement to the classification rules adopted for Tab. 3), but the position and the time (end of October) might lead to the conclusion that this anomaly is connected to the change in regional seismic style reported in   Table 1 and based on events with M L > 1.3 [39]. The lower limit of magnitude values of the events used by [39] for computing the seismic activity indices is significantly inferior to the threshold used here for NTG anomalies classification. Therefore, case-2 classification is still compatible with a possible connection to the change of very low magnitude seismic activity.

Seismic period -Central Italy -non AQ
In paragraph 3.1 and Figure 7b, it was discussed that, for the observation period during Spring 2009, Central Italy was divided into two regions based on the distance from L'Aquila epicentre: the non-AQ (Central Italy outside the most active region and far from the epicentre) and the AQ (portion of Central Italy closer the epicentre of the earthquake occurred in L'Aquila on 6 April 2009). In this paragraph, the results obtained for the non-AQ area are shown.
The  (Fig. 10). This anomaly can be classified as candidate precursor of the bigger event occurred in non-AQ area during the observation periods (and, in general, in all Central Italy except the mainshock and aftershocks activity), which was about to happen near Forlì on 5 April 2009 (just five hours before L'Aquila main shock) with a magnitude 3) was about to happen on 6 April (Fig. 10).
In Figure 10, another seismic event located in non-AQ, but far from Forlì is shown in the middle of the map, just outside AQ area. In the 15 days before its shake, no relevant NTG anomalies were detected in 10 pixels distance (60 km radius).
The 15 days prior to the earthquake time is in good agreement with previous researches mentioning as precursor event on the 21st, 22nd and 28th of March and 3rd of April [12]. Nevertheless, the most frequent condition in the data analysed for this study is typically significantly shorter. As one can see, in the map it is possible to notice the particular shape of Forlì NTG anomaly, which will be discussed further on (paragraph 4.3).
The other NTG anomalies observed in Central Italy during Spring 2009, they often had linear shape and showed a certain persistence of few days, even if they typically reach a positive saturation level only for a shorter time of 1-2 days (Fig. 11). Figure 11a presents two anomalies classified as case-1, one in the north and the other in the south (which lasted from, at least, 22 to 9 days before the earthquake). In Figure 11b a new NTG anomaly is shown which was classified as a case-1. In Figure 11c we have two anomalies classified as case-1.

Seismic period -Central Italy -AQ
During the observation period in Spring 2009 (27 January -8 May), two relevant anomalies were detected in the areas closer to the epicentre of L'Aquila event. Figure 12 summarises the timeseries of the anomaly connected to the main shock.  Significant temperature variability is observable from the 29th of March to the 6th of April (the night between the 6th and the 7th, i.e. one day after the earthquake). Figure 12a shows the regional trend of the NTG index values close to zero around the 30th of March in the same area of the future epicentre. In Figure 12b, a significant amount of pixels became NTG-positive reaching the maximum in the stack-nights between the 2 and 3 of April (Fig. 12b), and 3 to 4 April (not in the Fig. 12). Figure 12c shows the NTG map of the earthquake night and is characterized by a weak positive anomaly. The days after the earthquake the NTG values come back to zero [17,21].
The days after the 7th of April, the maps come back to general negative NTG values, and maintain this trend up to the end of April, when a new elongated anomaly appears parallel to the line of the Apennines mountains (Fig. 13a). As shown in Table 2, all the most relevant aftershock activity happened before this anomaly, so this pattern cannot be linked to the major activity following the main shock. Despite this, considering the minor seismic activity, it is possible to find an extremely good agreement between the NTG maps (Fig. 9a) and very minor seismic activity (Fig. 13b), [17,21]. The absence of seismic catalogue declustering leads to its classification as case-1.

Some geospatial remarks and discussion
Based on the results obtained for the two anomalies potentially related to the earthquakes of L'Aquila (main shock) and Forlì, it is possible to draw some conclusions.
From the L'Aquila NTG anomaly, three papers, authored by two different research groups, moved to relatively high resolution geospatial conclusions by overlapping the  thermal anomalies of the 3rd of April on layers describing: (i) major tectonic layout (Fig. 14a, [22]); (ii) lithology, faults, digital elevation models (Fig. 14b, [22,23]); (iii) soil coverage (Fig. 14c, [24]). In Figure 14a, it is possible to notice that two linear tectonic structures (the Antrodoco and Olevano Line and the Gran Sasso Thrust) separating major geological units of the area significantly match with the boundaries of the thermal anomalies inferred by the proposed NTG procedure. It is noticeable that, even if the area is generally characterized by extensional normal fault activity (for instance, the Paganica Fault -red circle in Fig. 14a -which breaking led to L'Aquila mainshockviolet circle), the two tectonic lines have a compressive geometry [22]. Moreover, as expected by literature and diffusely explained in [22], Figure 14b [23] the thermal features correspond well with the pronounced topography of the area. In addition, Figure 14c [24] shows the superimposition of NTG anomalies and land cover (left) and the geological and aquifer features (right). Analysing those lines of evidence, the authors suggest that even a degassing activity could have contributed to the NTG anomaly formation. Figure 15 proposes a more detailed visual analysis of topographic effect for the two anomalies of Forlì and L'Aquila, from which a consistent contribution of electronic charges movement up to topographic edges is confirmed as plausible [22]. In

(b) (a) (c)
Fig. 14. Superimposition of major tectonic features (the Antrodoco and Olevano Line and the Gran Sasso Thrust -solid black lines) and the NTG anomalies possibly related to the L'Aquila main shock -violet circle (14a, [22]); 3D visualization (aerial perspective from South) of the digital elevation model with the NTG map superimposed (14b, [23]); Comparison of the NTG anomalies against the land cover (left) and the geological and aquifer features (right) in the area and reported, in literature, as sign of potential degassing origin of L'Aquila NTG anomaly (Fig. 14c, [24]). addition, the topographic effect is also analysed in Figure 16 where transects of maps in Figure 15 are reported. potentially related to L'Aquila main shock. Here the topographic effect is still present (as can also be seen by Fig. 14b) but partially modulated by the effects of Figure 14a. The two diagrams on the right of Figure 16 show two NTG and topographic profiles for the same maps located at rows 52 and 54. Even in this case, positive correlation is often true, as evidenced in Figure 16 by the dashed vertical blue lines, but in few cases NTG local maxima correspond to valleys or mountain sides (dashed magenta vertical lines). A comparison of thermal and higher resolution topographic profiles for the same scene is published also in [22], where it is possible to notice the good agreement between the two curves even for local maxima with negative NTG values close to zero.
Further NTG maps and profiles should be analysed and a general correlation rule between high positive NTG values and mountains edges is for now impossible (only 4/5 of NTG anomalies of Fig. 11 are located on mountains) but, especially for the two main cases discussed in depth (L'Aquila and Forlì), the results are in good agreement with F.T. Freund's theory of positive hole charge carriers, stress-activated deep below, in the hypocentre, which would be able to spread out through the rock column. Upon arriving at the surface of the Earth two phenomena are expected to happen: (1) the positive holes tend to accumulate at topographic highs and (2) the positive holes recombine, returning to the peroxy state, an exothermal reaction, which leads to the emission of IR in the thermal band [42][43][44][45].
Regarding NTG maps with anomalies referred to L'Aquila main quake, it was noted [22] that the two main positive NTG patterns were characterised by different concentration of seismogenic or general fault systems (one had a certain concentration of faults while the other not) with a prevalence of NTG positive values on limestone lithology (76% of all positive values were detected on limestones while only 48% of the lithological map had this classification). Furthermore, most positive NTG values close to saturation were detected in both limestones and sandstones. Anyway, the study of correspondences between lithology and NTG anomalies potentially related to earthquakes should be investigated on more case studies in order to reconstruct a general law of the effect of rock types on NTG anomalies. In the same publication [22] and reported here in Figure 14a, a spatial effect was noted which was brought back to the presence of the two main tectonic discontinuities (Antrodoco-Olevano line and Gran Sasso thrusts) whose geometries were opposite to the mechanisms (normal faults) of the World Stress Map for the area and also to the Paganica Fault rupture of the earthquake of 6th April 2009 on L'Aquila.
As shown in Figure 14c, the NTG anomaly before L'Aquila main shock was also in good agreement with broadleaved deciduous forest distribution for which was observed that forests are a preferential route for degassing activity, especially CO 2 [24]. As for the other geospatial considerations previously discussed, the same analysis should be applied to more case studies to validate the first observations done. Traditionally explained by dilatancy theory [46], earthquake related degassing activity has also been recently brought back to peroxy defects theory by [47].
Furthermore, the similarities of Figure 13a with Figure 13b, in which NTG anomalous pattern for the 30th of April 2009 is compared with very low aftershock activity (M L > 2), are still an open issue. In order to give an explanation of why this phenomenon can be observed here and not in other cases (both for the spatial considerations and for the weakness of seismic activity in comparison with the previous one summarised in Tab. 2) more investigations are needed using larger datasets.
Similarly, even the apparent variation of NTG sensitivity with respect to the seismic activity magnitude as a function of the local background seismic condition and its change appears as an important topic worth of further investigations. In fact, as reported in paragraph 4.2.1, the anomaly of Figure 9 was observed together with a change in the background seismic activity. Moreover, as summarised in Table 3 (paragraph 4.1) and exemplified by the absence of NTG anomaly in the middle of the map in Figure 10 (paragraph 4.2.2), the possible influence of background local seismic activity condition on the magnitude threshold for the appearance of NTG anomalies can be observed. In particular, for the highlighted case in Figure 10, this could be due to the proximity to AQ region and to the observed dependence of NTG sensitivity to the level of local seismic activity. These observations on NTG apparent sensitivity to the background seismic condition are new to TIR anomalies literature and have to be validated on larger datasets in order to evaluate their general validity.
Finally, it is important to underline the complexities of non-seismological monitoring of seismic areas that add variability to experimental results: -Earthquakes are a consequence of stresses on very complex mechanical systems, where not every stress condition leads to the same seismic effects; -Satellite thermal remote sensing is an estimation process of ground thermal parameters based on the compensation of many atmospheric and surface conditions; -Stresses deep below not only influence mechanical effects but do also trigger other physical perturbations which are supposed to be able flowing to the ground surface interacting with the rock volumes gone across and with the ground-atmosphere interface.

Conclusions
In this paper the Night Thermal Gradient method is presented in depth. First, potential strength and weakness points are described. The most important positive aspect is the radiometric resolution enhancement originated from the synthesis of multiple original raw data into only two parameters with respect to the only one-sample-per-night approaches. This expected enhanced radiometric resolution leads to more readable and less noisy maps with thermal patterns that are more stable, spatially and temporarily. Thanks to this, the author and a different research group were independently able to analyse spatial features possibly linked to thermal anomaly patterns with greater detail than it was allowed before. On the opposite, when using NTG techniques, one must be aware that the processing workflow can lead to artefacts due to temporal distribution of raw data or to meteorological conditions (i.e. Foehn winds or only partially clear sky). It is also important to notice that the processing steps change the frequency content of the raw thermal signal and different parameters must be optimised. Furthermore, NTG outputs can be postprocessed with typical statistical techniques of timeseries analyses and their integration with raw data is very important.
The method was applied to study possible earthquake precursors concerning L'Aquila 2009 event. The experimental results include maps related to background and most active seismic condition analyses for different nearby areas: Central Italy in a nonactive period (October 2008), Sardinia region (in both October 2008 and Spring 2009) and Central Italy in an active period (Spring 2009). The related maps report thermal anomalies classified based on the presence of low magnitude events, spatially and temporarily close to the anomalies. In particular, one thermal anomaly is shown for the 21st of March 2009, which is consistent with time and location shown in other researches and that can be linked to an earthquake happened near Forlì the day before L'Aquila main shock. Other thermal anomalies are shortly presented and commented.
Furthermore, some geospatial remarks are discussed, reviewing the literature and adding a new graphical comparison between land elevation and thermal anomalous patterns for the cases of Forlì and L'Aquila. Results mostly show a positive correlation between elevation and thermal anomalies, with higher positive values of NTG corresponding to topographic heights. For the L'Aquila case, probably, the NTG anomalies are not symmetrically distributed around the epicentre partially due to the mountains and forests distribution and to the influence of two thrust structures.
Finally, the first application of the NTG method has proved to be effective and able to offer new potential tools for earthquake precursors applications. The present study can pave the way to future investigations exploiting high temporal resolution data, for instance, reducing environmental noise affecting raw data, better understanding the links between NTG values and mechanical parameters of earthquakes, correlating thermal anomalies and incoming earthquakes.
Open Access funding provided by Universit degli Studi di Cagliari within the CRUI-CARE Agreement. The author is grateful to M. L. Foddis (University of Cagliari) for fruitful discussions and to the two reviewers for the help in improving manuscript quality.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.