Satellite-Observed Thermal Anomalies and Deformation Patterns Associated to the 2021, Central Crete Seismic Sequence

Abstract

Nowadays, there has been a growing interest in understanding earthquake forerunners, i.e., anomalous variations that are possibly associated with the complex process of earthquake evolution. In this context, the Robust Satellite Technique was coupled with 10 years (2012–2021) of daily night-time MODIS-Land Surface Temperature remote sensing data to detect thermal anomalies likely related to the 27 September 2021, strong onshore earthquake of magnitude Mw6.0 occurring near the Arkalochori village in Central Crete, Greece. Eight intense (signal-to-noise ratio > 3) and infrequent, quite extensive, and temporally persistent thermal signal transients were detected and characterized as pre-seismic anomalies, while one thermal signal transient was identified as a co-seismic effect on the day of the main tectonic event. The thermal anomalies dataset was combined with tectonic parameters of Central Crete, such as active faults and fault density, seismogenic zones and ground displacement maps produced using Sentinel-1 satellite imagery and the Interferometric Synthetic Aperture Radar technique. Regarding the thermal anomaly of 27 September, its greatest portion was observed over the footwall part of the fault where a significant subsidence up to 20 cm exists. We suggest that the thermal anomalies are possibly connected with gas release which happens due to stress changes and is controlled by the existence of tectonic lines and the density of the faults, even if alternative explanations could not be excluded.

1. Introduction

Recently, a growing interest in the physics of precursors that could be used to develop a system for short-term (days to weeks) forecasting of strong earthquake events and support the understanding of the evolution of seismicity exists [1]. To create a multiparametric prediction system, it would first be necessary to identify the parameters whose anomalous variants might be potentially related to the complex and dynamic tectonic processes [2,3].

Among the various parameters, variations in the Earth’s emitted thermal radiation measured by space-born sensors operating in the thermal infrared (TIR) spectrum and the variations on the Earth’s surface deformation recorded by satellite radar interferometry have long been suggested as potential precursors [3,4]. Nevertheless, although much research has been undertaken on the investigation of several precursors simultaneously (e.g., [5]), the analysis of patterns of the different parameters has attracted less attention.

Several satellite sensors, which have channels in the infrared portion of the electromagnetic spectrum, can monitor Earth’s thermal field and generate a Land Surface Temperature (LST) product with variable spatiotemporal resolution depending on the satellite/sensor system chosen [6,7,8,9]. Therefore, LST can be used as a critical guide to the understanding of land surface processes at various scales [10,11,12,13,14,15,16,17,18,19,20] and environments, especially in areas with a variability in geophysical parameters such as vegetation, topography, geology and geomorphometry.

On 27 September 2021, a rare and strong onshore earthquake of magnitude Mw6.0 occurred in Central Crete, near the village of Arkalochori at 06:17:21 UTC [21]. The epicenter of the seismic event was located at latitude 35.15 N and longitude 25.27 E, while the focal depth was about 10 km and associated with the Kastelli Fault Zone [22] whose mean recurrence interval has been estimated at approximately 800 years [23]. Since the beginning of June, almost 4 months earlier, a significant number of foreshocks had beenrecorded in the broader area, while a rich aftershock sequence was observed in the months after the main shock occurrence [21,22,24,25]. On 24 July, a strong preseismic event appeared with magnitude Mw4.8. After that, more than 250 earthquakes with magnitudes up to 3.8 were recorded in an approximately N-S-oriented, 15-km-long area west of Arkalochori until the occurrence of the main event [24,26]. Eight aftershocks of magnitudes equal to or greater than M4.2 occurred during the day of the main event, while the strongest aftershock occurred on 28 September with Mw5.3. The earthquake had devastating consequences in the area near the Arkalochori village as one person lost his life, while many buildings collapsed affecting thousands of people.

The objective of the present work is to identify thermal anomalies associated with the rare major tectonic deformation event of Central Crete and to correlate them with co-seismic regional crustal deformation. Thermal and crustal deformation anomalies will be calculated using data mining techniques over a period of three months before and one month after the occurrence of the earthquake; similarities or discrepancies between the behavioral changes of these anomalies will then be analyzed. Based on a preliminary analysis, we note that in the period May 2021 to July 2021, where the first foreshock appeared, no thermal anomalies were observed in the region.

The Robust Satellite Techniques (RST) [27,28,29] was used to detect the spatial evolution of thermal anomalies probably connected with the evolution of the seismic sequence. During the last two decades, the RST has been widely used to distinguish thermal signals possibly related with earthquake evolution from variations in the Earth’s thermal emission due to other phenomena [27,28,29]. For the RST, any anomaly in the space–time domain is considered as a deviation from a “normal” behavior. The latter can be calculated by processing a multiyear dataset of cloud-free satellite images of the same month, the same spectral channel (s) and the same acquisition times.

In addition, the regional co-seismic deformation was mapped using the Interferometric Synthetic Aperture Radar (InSAR) technique. The term SAR stands for Synthetic Aperture Radar [30]. SAR systems that record microwave radiation are called active because of their ability to emit pulses of microwave radiation lasting a few minutes, thus being autonomous sources of energy. The recording of the pulse return time by the satellite sensors determines the position of the earth features on the image. Compared to conventional geodetic methods, the ability of radar to measure distances through time and at an angle, over long distances from either two different systems or with repeated shots of the same system was a new technique, Interferometric SAR (InSAR). Differential Interferometry (Differential InSAR technique-DInSAR) is an advanced technique [31] aimed at detecting surface movements due to geophysical phenomena or human interventions. Since the 1990s, the DInSAR technique has proven to be an interesting tool for measuring and observing ground deformation suitable for analyzing geodynamic processes (e.g., [32,33,34,35]). In this context, we used SAR images of the Sentinel-1 Copernicus satellite in both ascending and descending orbits to create the differential interferograms and to produce deformation patterns through phase unwrapping, using open ESA’s SNAP software.

The contribution of the present work is to study thermal anomalies observed before a very rare and strong earthquake event. Its innovation is that it associates the observed thermal anomalies with regional tectonic characteristics. We note that the thermal anomalies mapped in the area are coupled with co-seismic deformation patterns observed by the InSAR method.

2. Geotectonic Setting

Crete is situated at the southernmost part of the Hellenic Island arc and has high seismic activity due to the subduction of the African plate underneath the Eurasian plate [36,37,38]. Both extensional and compressional stresses that exist in this area result in the creation of an extremely complex tectonic environment characterized by significant horizontal and vertical movements [39]. The Heraklion tectonic graben to the north and the Messara tectonic graben in the southwest dominate the central Crete region on shore (Figure 1).

The Heraklion graben is bounded on the west by the Ida Mountains and on the east by the Dicti Mountains along the eastern Psiloritis and Kastelli fault zones, respectively, [39,40,41]. Both fault zones exhibit an almost NE–SW direction, with the eastern Psiloritis fault zone being more prominent due to the topography of the Ida Mountains. The Messara graben is bounded to the north by the Ida Mountains and to the south by the Asteroussia Mountains, along two EW oriented fault zones (Figure 1). According to previous studies [38,41,42] the faults on Central Crete can be classified in four groups. The first group consists of E–W trending faults that cut mainly bedrock or bound bedrock and Miocene sediments. The second group consists of large and medium N–S striking faults that cut the aforementioned group. The third and fourth groups are trending towards NE–SW and NW–SE, respectively. The geology of the area of interest is composed mainly of Miocene to Pliocene sediments which cover the central and western part of the study area overlaid by Quaternary deposits forming depositional plains. The Neogene deposits were formed in terrestrial to deep marine environments, giving evidence for significant changes in depositional depths attributed to vertical movements and climatic changes [43]. Carbonates of the Tripolis nappe are exposed in the northeastern part of the area. Dissected hills of phyllites and quartzites coupled with Plattenkalk limestones are mainly outcroppings to the eastern part of the area of interest [43].

3. Materials and Methods

3.1. Robust Satellite Techniques

The Robust Satellite Technique (RST) [27,28,29] is applied herein to distinguish anomalous thermal variations probably related to earthquake activity from normal variations in Earth’s thermally emitted radiation due to other phenomena. The RST is a well-known multitemporal remotely sensed data analysis technique having as its main advantage the reduction of “natural noise” [28]. The cause of this “natural noise” can be the atmospheric weather or several variations in vegetation (natural or planted) or man-made constructions, topography and climate over space and time. As shown in previous works (e.g., [2,19,20,27,28,44,45]), the advantage of the RST methodology compared to other approaches is its ability to isolate residual thermal variations which are probably connected to tectonic activity. Over the past two decades, the RST method has been applied worldwide in several geotectonic regimes and to earthquakes with magnitudes greater than 4.0. Both polar and geostationary satellite images have been used, and connections between thermal anomalies and tectonic activity in a temporal window from 30 days before to 15 days after earthquakes occurrence have been found ([2] and references therein). Among others, TIR anomalies associated with the 7 September 1999, M5.9 earthquake in Athens [44], to the M7.8 Izmit earthquake on 17 August 1999 [28] and to the M5.8 earthquake in Abruzzo (Italy) on 6 April 2009 [46] have been reported. Moreover, the 10 years analysis of continuous data over Greece by [45] showed the general correlation between earthquakes and the thermal anomalies computed by the RST. A similar investigation of RST capabilities recently took place in Turkey and confirmed a nonrandom correlation between RST-based thermal anomalies and earthquake occurrence [20].

In the present work, the RETIRA-index computed as [28]:

$${\otimes }_{\mathsf{\Delta }\mathsf{{\rm T}}}\left(r,\mathrm{t}\right)=\frac{\mathsf{\Delta }\mathsf{{\rm T}}\left(r,\mathrm{t}\right)-{ \mathsf{\mu }}_{\mathsf{\Delta }\mathsf{{\rm T}}}\left(r\right)}{\begin{array}{c}{\mathsf{\sigma }}_{\mathsf{\Delta }\mathsf{{\rm T}}}\left(r\right)\\ \end{array}}$$

(1)

where r = (x,y) represents the accurate location of each pixel (x,y) on the satellite image, t is the acquisition time of the satellite image, ΔΤ(r,t) refers to the difference (T(r,t) − Τ(t)) of the observed TIR signal value T(x,y,t) with the spatial mean value T(t) of all the pixels of the satellite image. T(x,y,t) is measured for each pixel (r), while T(t) is calculated in place on the satellite image, excluding the cloudy pixels$;{ \mathsf{\mu }}_{\mathsf{\Delta }\mathsf{{\rm T}}}\left(r\right)$ defines the time average and ${\mathsf{\sigma }}_{\mathsf{\Delta }\mathsf{{\rm T}}}\left(r\right)$ the standard deviation of ΔΤ(r,t) measured at site r, computed only on cloud-free pixels of the homogenous satellite images.

The RETIRA index calculates the local spatiotemporal excess of the current ΔT(r,t) signal over its historical mean value, weighted by its variability σ_ΔΤ(x,y) at the given location. Using ∆T(r,t) instead of T(r,t) reduces the possible contributions due to meteorological variations (e.g., periods of abnormally hot weather) and/or seasonal time changes. Note that the signal (S) is evaluated by comparing it to the standard deviation σ_ΔT(r); therefore, the signal-to-noise ratio (S/N) can be used to characterize the intensity of anomalous thermal transients.

For the RETIRA index calculation, the MOD11A1 Version 6.1 product, available from 24 February 2000 was used. It provides daily per-pixel “Land Surface Temperature and Emissivity” with 1 km pixel size in a 1200 by 1200 km grid. The pixel temperature value is derived from the MOD11_L2 (https://0-doi-org.brum.beds.ac.uk/10.5067/MODIS/MOD11_L2.006, accessed on 1 November 2021) swath product [47]. The retrieved MODIS Land Surface Temperature uses the MODIS cloud mask product and as a result contains only cloud-free pixels. Images from LST were preferably taken at night because they are affected to a lesser extent by temperature differences between the ground and air than images taken at any other time of day. In addition, they are also less sensitive to local solar radiation variations and shade, an important source of land surface temperature variability not related to the seismic activity. Specifically, 10 years of remotely sensed data (from 2012 to 2021) were included in the calculation of the RETIRA index since the use of multi-year time series of satellite data is crucial to establish the required regional thermal background. Initially, the “ Night-time and Surface Temperature” layer dataset was isolated for the 10-year period and then spatially subset to the Central Crete area. Monthly reference fields μ_∆T(x,y) and σ_∆T(x,y) were then generated for the four months period, i.e., from July to October. Nevertheless, it has been shown that in cases where the cloudy fraction of the scene is greater than 80%, the remaining values lead to a non-reliable signal [45]. For this reason, such images were excluded from the computation of the final reference fields.

To identify any possible variations, a RETIRA index was calculated for each day of the period of 3 months before and a month after the main earthquake. We note that RETIRA is a Gaussian standardized variable, and the selection of its relative threshold value, quantitatively reflects the rareness and the significance of the identified anomalies. In this work, to define a thermal anomaly, we used a RETIRA index greater than 3, following a very strict rule as a distribution value greater than 3 or 4 suggests an anomaly with a probability 99.7% or 99.99%, respectively. In some cases, the computation of thermal anomalies is affected by relatively rare physical phenomena such as heatwaves, wildfires, cloud coverage and/or observational changes as in the case of inaccurate image navigation/co-location. The RETIRA index is based on time-averaged quantities and as shown in previous research [28,44,45] is sensitive to the abrupt appearance of extreme signal due to the abovementioned natural [48] or anthropogenic [44,49] phenomena. However, these signal variations are characterized by specific spatiotemporal patterns such as limited duration, similar and expected spatial distribution and known date of occurrence, allowing their exclusion from further analysis [28]. As a result, the RETIRA index anomalies related to seismicity were isolated from those due to wildfire and those induced by cloud-coverage or/and imprecise image navigation/co-location.

3.2. Interferometric Synthetic Aperture Radar (InSAR)

The Sentinel-1 mission of the European Commission (EC) and the European Space Agency (ESA), comprises a constellation of two polar-orbiting satellites, Sentinel-1A and Sentinel-1B, operating 24 h per day with a short revisit cycle of 6 days. They are performing C-band synthetic aperture radar imaging under all weather conditions [50]. To map the co-seismic deformation due to the 27 September 2021 earthquake, we used one ascending and one descending SAR image pair (with acquisition dates of 23 September 2021–29 September 2021 and 25 September 2021–1 October 2021, respectively) from ESA’s Sentinel-1A and Sentinel-1B satellites. In both cases, the reference image was the one before the earthquake occurrence, while the repeat image was the one that refers to the date after the event.

Each reference–repeat pair was processed using the ESA’s open SNAP software, and two individual interferograms were generated. The topographic phase was subtracted using the SRTM 1 arc-second Digital Terrain Model, a 30 m resolution Shuttle Radar Topography Mission Digital Elevation Model (USGS 1 ARC-second SRTM DEM, https://doi.org/10.5066/F7DF6PQS), while the signal-to-noise ratio was enhanced by applying the adaptive power spectrum filter of [51] with a coherence threshold of 0.4.

4. Results

We proceed now to present the results achieved after applying the RETIRA index and to compare them with the co-seismic crustal deformation pattern.

For a period of 10 years (i.e., 2012–2021) and for each of the 4 months spanning the period before and after the main event (i.e., July to October), 1230 images were obtained. Each cloud-free image contains 3525 pixels. The cloud-covered pixels in the images are recognized as “no data”, and as a result, pixels reduce as the cloud coverage increases. Finally, 1086 images were used in the calculation of the monthly reference fields μ_∆T(x,y) and σ_∆T(x,y). Then, the RETIRA index was computed for 105 images (i.e., for 105 of the 123 days in the period July to October 2022). Among the images used, 31 were in July, 31 in August, 28 in September and 15 in October. The remaining 20 images during September and October 2021 presented cloud cover greater than 80% and were therefore excluded from further analysis. The daily analysis led to the detection of 18 (out of the 105) images containing pixels with RETIRA values greater than three. Applying the previously described methodology to all MODIS-LST images from July 2021 to October 2021, we finally identified nine thermal anomalies, three for each of the months of July, August and September 2021 (Figure 2, Figure 3 and Figure 4). During October 2021, there was no observed RETIRA index following the post-seismic thermal anomaly criteria.

In detail, eight quite intense (Signal/Noise > 3.0) and rare, spatially extensive TIR signal transients were identified before the main earthquake. The first was recorded on 2 July, a few hours after the occurrence of a M4.1 earthquake and three days before the 5 July M4.2 event. The second thermal anomaly was recorded on 30 July, 6 days after the 24 July M4.8 earthquake. A new thermal anomaly was mapped one day later (Figure 2). During 5, 10 and 29 August, three thermal anomalies were identified (Figure 3). We note that during August, the earthquake activity appeared with lower magnitude events, with the strongest on 8 August with a magnitude M3.8. The last two pre-seismic thermal anomalies were identified on 18 and 20 September, nine and seven days before the main event on 27 September, respectively, and the last thermal anomaly was detected on the day of the main earthquake (co-seismic anomaly) (Figure 4). As can be seen in Figure 2, Figure 3 and Figure 4, in some cases the RETIRA index exceeded the value five (see dark brown pixels), and this could be used as an evaluation of the intensity of the thermal anomalous events. As for the spatial distribution of the observed thermal activity in Central Crete, the thermal anomaly of 2 July occurred to the south and east of the village of Arkalochori and occupied the southern part of the Kastelli seismogenic zone. In the northern part, the thermal anomaly spatially coincided with Mount Yuchtas, where a conjunction of two active tectonic faults exists. The 31 July thermal anomaly was also located near Mount Yuchtas and south of Arkalochori, while the thermal anomaly observed on 30 July was located to the south and was bounded by the two large faults, oriented east–west (Figure 2). On 5 August, one more thermal anomaly appeared near Mount Yuchtas northwest of Arkalochori. An extended thermal anomaly was identified on 10 August and occupied the previous area located further on north, south and southwest of Arkalochori. On 29 August, a much weaker, spatially limited thermal anomaly was located near Houdetsi village (Figure 3).

On 18 September, a thermal anomaly occurred north of the earthquake epicenter. Two days later, on 20 September, a RETIRA thermal anomaly with values up to 6.0 was observed over Mount Yuchtas, with the largest values occurring east of the mountain and extending to the south (Figure 4).

On the day of the earthquake occurrence, the thermal anomaly was located near an east–west oriented zone south of the epicenter, occupying an extensive area bounded by the villages of Metaxochori to the north and Ligortynos to the south, reaching as far as the vicinity of Arkalochori (Figure 4).

To visualize the frequency of thermal anomalies and their spatial distribution, Figure 5a was constructed. With the aim to define the areas where the thermal anomalies were spatiotemporally persistent, a map showing the frequency of the occurrence of thermal anomalies (i.e., the number of times a pixel participated in one of the nine identified thermal anomalies during the period from July to September 2021) was created. The frequency map of thermal anomalies was then compared with the normalized map of fault density obtained by calculating the density of linear features (faults) in the neighborhood of each output raster cell according to [52]. More specifically, a normalized lineament density map was created by counting the faults per unit area (number/km²), dividing the output by the maximum calculated density and then plotted in the respective grid centers.

The frequency ranges from one to four indicating that especially in the vicinity of Mount Yuchtas, the thermal anomalies were particularly persistent, as they occurred on four of the nine days in total. The frequency map of thermal anomalies was then compared with the normalized fault density map (Figure 5b). It is worth noting that the region north of the earthquake epicenter exhibits a medium-high value of fault density and the highest persistence of thermal anomalies, in contrast to the area south of Arkalochori, which shows anomalies with medium persistence and low fault density (Figure 5).

We proceed now to present the two co-seismic interferograms that emerged for the Arkalochori earthquake of 27 September 2021. The two wrapped (ascending and descending) interferograms are of good quality due to the low temporal geometric baselines and the high temporal resolution of Sentinel-1. As a result, the interferograms exhibit high coherence, with 80% and 68% at ≥0.4 and 0.6, respectively. They contain the phase difference between master and slave images produced by the main seismic event and its aftershocks until 29 September (ascending orbit) and 1 October (descending orbit). A clear pattern of six fringes, i.e., six concentric “circles” quite similar in shape, forming a lobe (Figure 6a,b), which is misaligned mainly due to the different geometry of the images in the two different orbits, is evident, indicating subsidence. Each interference fringe is a phase change which corresponds to motion of 2.8 cm in the satellite line of sight. The two wrapped interferograms were then used in the second stage; to calculate the terrain displacement, an unwrapping process was performed, and the phase unit was transformed into distance units in the line of sight (LoS) (differential interferograms) for every interferometric pair.

The InSAR results show that the area surrounding Arkalochori moved roughly downwards and are in a good agreement with previous works [21,22,53]. Note that dip-slip earthquakes have a very asymmetrical displacement pattern: i.e., the subsidence will be larger than the uplift for normal-faulting earthquakes. The terrain displacement products are similarly shifted. In the ascending orbit, we observe the epicenter of the main earthquake in the center of the ground deformation which has a maximum value of 18 cm (Figure 6c) which is the result of the superimposed effect of all the seismic events, mainly between 27 September and 29 September. In the descending orbit, the ground deformation has a maximum value of about 20 cm which results from the seismic events mainly between 27 September to 1 October and is located mainly east of the epicenter of the main earthquake (Figure 6d).

Finally, to map the vertical (up–down) and horizontal (east–west) displacement, decomposition of ascending and descending LOS displacement products was performed. After the decomposition, the ground deformation in the vertical (up–down) direction showed subsidence up to 20 cm, while no uplift displacement was detected (Figure 7a). The horizontal (east–west) displacement map indicated a maximum displacement of about 8 cm towards the east of the footwall block and about 6 cm towards the west of the hanging wall (Figure 7b) in agreement with that observed in [21,22] and modeled in [53], probably attributed to the strike–slip component of the earthquake focal mechanism as suggested in [21].

5. Conclusions

In the present work the RST (Robust Satellite Technique) was used to detect and map thermal anomalies possibly associated to the 27 September 2021, strong, rare, and unexpected earthquake of magnitude Mw6.0 occurring onshore of Central Crete (Greece), in the vicinity of Arkalochori village. A total of 10 years (2012–2021) of daily night-time Land Surface Temperature (LST) remote sensing data from Moderate Resolution Imaging Spectroradiometer (MODIS) were used to create the background pattern of the area. The seismic activity began several months before the main event, and nine intense (signal-to-noise ratio > 3), infrequent, spatially extensive, and temporally persistent TIR signal transients were identified from July to September 2021. Eight of them could be characterized as pre-seismic thermal anomalies, while the last one as co-seismic, since it was observed the day the main event occurred. Concerning their spatial distribution, five of the pre-seismic thermal anomalies coincided with the Mount Yuchtas area, while two (observed on 2 July and 10 August, respectively) were located south of Arckalochori village in the proximity of the Kastelli seismogenic source. Similar spatial distribution existed in the thermal anomaly observed during the day the main event occurred. The latter possibly implies that the two thermal anomalies observed south of the earthquake epicenter before the earthquake occurrence could be related to the future co-seismic displacement of the region due to the activation of the Kastelli zone.

Moreover, using InSAR patterns, we confirmed that the Mw = 6.0 earthquake and its strongest aftershock of 28 September 2021 (Mw5.3) produced mainly subsidence up to 20 cm combined with a less important horizontal eastward movement up to 6 cm and a westward movement up to 8 cm to the west and to the east of Arckalochori town respectively.

The superimposition of the frequency thermal anomaly map to the ground deformation maps (produced using Sentinel-1 satellite images) (Figure 8) reveals the existence of a common area in the south of Arkalochori village. As far as we know, this is the first time that the spatiotemporal correlation between the two different phenomena, i.e., thermally emitted radiation anomalies and the Earth’s surface deformation, has been revealed.

Concerning the day the main event occurred, the two datasets partially overlap with a great portion of thermally anomalous pixels to be contained by the area that significantly deformed. This is more evident in the case of horizontal displacement where a significant portion of the thermal anomaly coincided with the westward-moving terranes and the southern part of the Kastelli fault zone.

What is worth noting is the fact that part of the deformed area was characterized as thermally anomalous twice during the four months (Figure 8). This observation, although needing further investigation, could be of great importance for the development of a short-term forecasting system of strong earthquake events.

Herein, we attempt to propose a preliminary physical mechanism for the observed thermal anomalies in agreement with the deformation pattern associated to the Mw6.0 main event. Thermal anomalies seem to exhibit a selective spatial distribution, as they appear to be preferentially concentrated in areas with high fault densities, such as areas where two faults are joined, as in the case of Mount Yuchtas. This significant increase of thermal radiation can be attributed to greenhouse gas (mainly CO₂ and CH₄) reservoirs whose spatial relationship with seismogenic fault systems has been reported (e.g., [54,55]). As stated by [54], an increase of stress can cause the release of gases, such as CO₂, CH₄ and N₂, that are trapped in the pores of the rocks. The gases are released from the Earth’s surface with the activation of cracks since micro-fracturing can lead to permeability changes and gas diffusion, which is supported by their high mobility [55,56]. In some cases, and especially in extensional tectonic environments, CO₂ is produced from the melting of subducted carbonates and accumulates in crustal reservoirs triggering earthquakes at crustal depth [57,58].

An alternative explanation which could be active in parallel with that proposed above is based on the activation of electrical charges (electrons and positive holes) in rocks by the activation of peroxy defects during the increasing of tectonic stresses prior to any major seismic activity [59]. The positive holes are able to flow out of the stressed rock volume through the surrounding unstressed rock traveling fast, while electrons remain trapped locally [60]. The presence of hole charge carriers that form a surface charge layer could be an alternative explanation since when holes accumulate in the surface or near surface structure, they can recombine and return to the self-trapped positive hole pair from which they exited before their stress activation at depth. Furthermore, laboratory experiments till failure on a constant rated load anorhosite sample (a nominally monomineralic feldspar rock) suggests an IR emission that appears and disappears as the stress build-up progresses, interpreted as a hole pair recombination process ([59,60] and references therein).

The hypothesis of thermal anomalies due to the release of gases is suggested for the Central Crete earthquake sequence, and it is reinforced by the recent report of the presence of biogenic gases in water wells in the Messara basin [61] and the broader Arkalochori area as well [62] (Figure 9). The extended thick Neogene sediments are suggested as the host of the gas whose migration follows the general dip direction of the Miocene strata and the faults charging the sediments located several kilometers away to the NNW of the basin depocenter. This could explain the intense and persistent thermal anomalies on Mount Yuchtas. Given the fact that the Neogene sediments are extremely extended in the study area, with the western and the eastern parts of the Messara basin exhibiting the greatest depths, while the central part is significantly more elevated, one could similarly explain the existence of thermal anomalies southeast of Arkalochori and the absence of the phenomenon to the north and northeast of Arkalochori (Figure 9). The continuously increasing stress during the preseismic phase resulted to an increase in degassing activity. As the stress evolution creates a network of micro- and meso-scaled cracks that superimposed to the already existing faulting system, a pathway is created that supports the discharge of trapped gas buffers. As the earthquake approaches, due to the creation of a significant fracture in the rupture zone, a new rise in degassing activity can occur before normality is gradually restored [20]. The absence of thermal anomalies during October could indicate that after the earthquake occurrence, the stress field relaxed ([63] and references therein) and thus reduced the above-described processes. A crucial note is that the appearance of thermal anomalies was in the time period where the system presented critical dynamics [25].

Further analysis of the proposed mechanism could help in understanding the earthquake preparation process and may indicate a hopeful pathway of research. The RST technique in combination with InSAR images and geological data can play a very important role in this direction. In conclusion, the present work is an attempt that will increase the scientific effort to define earthquake forerunners, and its future fusion with other state-of-the-art methods (e.g., [5]) will support our knowledge for earthquake preparatory processes.