InSAR Time-Series Analysis for Monitoring Ground Displacement Trends in the Western Hellenic Arc: The Kythira Island, Greece

Abstract

Kythira Island is situated at the western Hellenic Arc, which is a region of very high seismicity and tectonic activity. On 8 January 2006, a large seismic event of M_w = 6.7 occurred close to Kythira, in association with the Hellenic subduction zone. We present an extensive multi-temporal interferometry study of ground displacements in Kythira Island exploiting SAR data of a decade, from 2003 to 2009, and from 2015 to 2019. The line-of-sight displacement field for the 2003–2009 time period presents small-scale displacements that do not exceed −3 mm/y, identified to the East of Kythira’s airport, and 3 mm/y in the central part of the island. We exploit then Sentinel-1 data from 2015 to 2019, of both descending and ascending Sentinel-1 SAR imaging modes to decompose the line-of-sight measured deformation to vertical and east–west motion components. Higher vertical displacement rates characterize the central part of the island, with a maximum value of 5 mm/y. This short-term uplift rate exceeds by an order of magnitude the long-term geologic uplift rate of ~0.13 mm/y found in the island during the Quaternary. We discuss possible regional geophysical explanations for this discrepancy. Strike-slip components are detected in the N–E coast of the Potamos village, where horizontal displacements occur, with an east-ward trend and a maximum value of −3 mm/y. These insights are valuable input for the systematic monitoring of this high seismic risk island and the dynamic assessment of its hazard potential.

1. Introduction

The Hellenic Subduction Zone (HSZ) is the most active geotectonic structure in the entire Mediterranean region. This is due to the active subduction of the Nubian lithosphere beneath the southern margin of the Eurasian lithospheric plate at the Aegean Sea area from about SW to NE e.g., [1] (Figure 1). In the region of the HSZ, the main morphological elements of tectonic origin from SW to NE are the Hellenic trench, consisting of a series of trenches, and the Hellenic arc, consisting of an outer sedimentary arc and an inner volcanic arc [2]. The HSZ region is characterized by very high seismicity, which geographically follows the arcuate shape of the structure. The active deformation along the HSZ has been extensively studied from different perspectives, including seismicity (e.g., [3,4,5]), focal mechanisms of earthquakes (e.g., [6,7,8,9]), tectonic observations (e.g., [10,11,12]), GPS measurements (e.g., [13]), and lithospheric modeling [14,15].

However, Interferometric Synthetic Aperture Radar (InSAR) techniques for the study of active deformation in the HSZ region have found only limited application, so far (e.g., [16]). One of the less studied areas in terms of its active deformation is the area of Kythira Island in the southwestern segment of the HSZ (Figure 1), which is historically characterized by very high seismicity [17]. The last strong earthquake (M_w = 6.7) that occurred in the island is that of 8 January 2006 [18].

In the last decades, new methods that exploit multiple SAR acquisitions, designed for long time-series analysis have been established as a powerful tool to efficiently map the evolution of the Earth’s surface deformation, applicable in a variety of physical processes [16,23,24], but also in human-induced deformation [25,26].

In this study, we present a SAR interferometry time-series analysis for Kythira Island covering a time interval of more than a decade. The objective of the study is to provide insights on the crustal deformation and the active tectonics of the area and better understand their association with the seismicity of the area. All available data-sets from 2003 to 2019 were analyzed and ground deformation velocity maps were produced. In addition to the displacements observed in the line-of-sight direction (LOS), the estimation of vertical and east–west motion components for the 2015–2019 time period provides an important asset for the interpretation of ground deformation in Kythira Island. A comparison with previous tectonic studies on the island is provided, aiming to better interpret active deformation and to highlight new displacements examined by the exploitation of Sentinel-1 imagery and their association with ongoing tectonic processes on the island.

2. Tectonic Setting

The island of Kythira, along with the Peloponnese and Crete Island (Figure 1), is part of the outer sedimentary arc in the southwestern segment of the Hellenic arc (e.g., [10]). The tectonics of Kythira is characterized by Late Cenozoic faults striking mainly NNW-SSE (e.g., [10,22]). However, no active faults have been detected on the island. In the historical and instrumental periods of seismicity, the broader area of Kythira Island experienced large magnitude earthquakes (Figure 1). The island was damaged by the earthquakes of 1750 and 1866, which likely have been shallow events. The area of Kythira strait was strongly affected by the (M_w = 6.6) intermediate-depth earthquake of 11 August 1903. Heavy damage was caused in Kythira Island, mainly in the Mitata village, while less damage was reported from the western Crete Island [17,27]. On 8 January 2006, another strong (M_w = 6.7) intermediate-depth earthquake [18] (Figure 1) caused damage and various ground failures (e.g., landslides) in Mitata and in other villages of Kythira (e.g., [22]).

3. Materials and Methods

Envisat and Sentinel-1 missions with their SAR instruments provided a dense and uniform data archive with an exceptional temporal coverage, suitable for deformation monitoring in Kythira Island. InSAR time-series were available with the processing of SLC/level-1 data, provided by the European Space Agency for the time interval from 2003 to 2019. A summary of data-sets information is provided in

Table 1. Data-set information.

	SENSORS
	Envisat ASAR	SENTINEL-1
No. of Acquisitions	18	37	36
Relative orbit	415	7	102
Time interval	2003–2009	2015–2019	2015–2019
Orbit pass	Ascending	Descending	Ascending

, while Figure 2 presents the satellite footprints of all available tracks covering the area of interest and employed for time-series analysis in the present study. We analyzed 21 Envisat ASAR SLC images from the ascending satellite pass no. 415, from 5 August 2003 to 24 February 2009. Recent data sets from Sentinel-1 sensor, enabled the combination from both ascending and descending satellite passes in Kythira Island. Sentinel-1 SLC images of Interferometric Wide (IW) swath mode were selected to map deformation for a time span of five years, i.e., between 2015 and 2019. Specifically, 37 images from descending track no. 7 and 36 images from the ascending track no. 102 were used in the implementation of the Persistent Scatterer Interferometry (PSI) technique.

Delft object-oriented interferometric software (DORIS version 4.06-beta2) [28], developed by the Delft Institute of Earth Observation and Space Systems, was used for interferogram generation with Envisat data. Our processing chain is named Parallelized Persistent Scatterer Interferometry (P-PSI) [29], a distributed PSI processor for the assessment of line-of-sight ground velocities with minimal operator interventions, tailored to Sentinel-1 data. Our P-PSI implementation invokes two main software packages that have been modified to parallelize the execution of the different processing tasks. We use the InSAR Scientific Computing Environment (ISCE) with its topsStack processor [30] for creating a stack of coregistered Single Look Complex (SLC) dataset. The SLC stack is then used for PSI analysis with a parallelized implementation of the StaMPS/MTI (Stanford Method for Persistent Scatterers—Multi-Temporal InSAR version 4.1-beta) [31,32], using the combination of Persistent Scatterers (PS) [33] and Small Baseline Subset (SBAS) [34] techniques. An SRTM DEM 1-arc seconds was used for the subtraction of the topography contribution [35]. Finally, we applied an atmospheric phase-based linear correction with the open-source Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) [36], in all interferograms for both time-spans investigated. Line of sight (LOS) displacements were estimated, where positive rates correspond to a motion towards the satellite and negative rates to a motion away from the satellite.

Data acquired by Envisat sensor were quite limited and a sufficient number of historical images to produce time-series were available only for one satellite pass. While all available tracks provided a good spatial coverage of the island, the only limiting factor was the number of scenes per track, which combined with the physical features of the island, restricted the final number of selected data for multi-temporal InSAR analysis. Although a typical number of 15–20 images from C-band is sufficient for PSI processing, bad image quality, for some of the scenes, resulted in stacks of less than 15 images per track. This restriction combined with the geomorphology and vegetation of the island could not ensure the desirable accuracy in the estimation of LOS velocities. Therefore, with the use of past data deformation, analysis on Kythira Island is limited to one satellite pass from Envisat sensor. The same limitations excluded the processing of SAR images, like data from ERS sensor, acquired before 2002.

Observations of the GPS station KITH, located near Kythira’s airport, were provided by Dionysos Satellite Observatory of National Technical University of Athens (NTUA). An area around this station served as a reference in the processing of all three data-sets. The location of KITH station is included in all deformation maps in the sections to follow. KITH station belongs to the NTUA/Oxford GPS/GNSS network and provided continuous observations from 2004 to 2010. The intermittent operation of KITH station since 2010 justifies its use as a reference point and not as a validation method for the deformation estimated during the 2003–2009 period.

Estimation of Vertical and Horizontal Motion Components

The decomposition of the LOS velocity vectors of the 2015–2019 time period to vertical and horizontal displacements contributes to identifying the origin of the observed motion with respect to the sensor and adds important information on ground displacements in Kythira Island. We use the methodology described previously [37] for the generation of Synthetic Permanent Scatterers that assumes that clusters of neighboring pixels belong to the same deformation regime, and the decomposition geometry.

Based on a previous study [37], a grid covering Kythira Island was formed, extending from 22.80° to 23.10° in longitude and 36.10° to 36.40° in latitude. Patches (50 m × 50 m) of the grid located in the sea, with zero PS scatterers, were excluded from the decomposition analysis. From a total of 33,007 scatterers from Sentinel-1 ascending track no. 102 and 53,661 scatterers from descending track no. 7, the decomposition resulted in 11,232 synthetic scatterers for both vertical and horizontal displacements. For decomposing the LOS displacements to vertical and horizontal components, the system of Equations (1) and (2) for the descending azimuth look angle direction [38], was solved for every synthetic scatterer. A uniform downsampling method was suitable for the type of deformation that Kythira Island experiences. The uniform type of grid enabled the assessment of deformation, not only in regions where pixel sampling is dense, but in the wider area of Kythira Island and resulted in vertical and horizontal velocity maps, which give an overview of the deformation mechanisms on the island. Although the selected grid size and the methodology followed lead to an undersampling of the data, this cell size was maintained, as the deformation in Kythira Island is smooth and subtle and the loss of information was not critical for the determination of displacements on the island.

For the heading angle rates, an average heading value for all pixels, which is around 349.31° for the ascending track and 190.32° for the descending track, were used.

$$\left[\begin{array}{c}{d}_{LOS}^{asc}\\ d_{LOS}^{desc}\end{array}\right]=A\left[\begin{array}{c}{d}_{Up}\\ d_{ALD}\end{array}\right]$$

(1)

$$A=\left[\begin{array}{cc}\cos \left({\theta }^{asc}\right)& \frac{\sin \left({\theta }^{asc}\right)}{\cos \left(\Delta \alpha \right)}\\ \cos \left({\theta }^{desc}\right)& \sin \left({\theta }^{desc}\right)\end{array}\right]$$

(2)

where d_LOS denotes the LOS deformation, d_Up is the vertical velocity, d_ALD is the projection of horizontal deformation in the descending azimuth look direction [37], θ is the incidence angle, and Δα is the satellite heading difference between ascending and descending pass. Moreover, the exponent asc and desc in d_LOS and θ, represent the ascending and descending satellite passes respectively.

The assessment of vertical and horizontal motion components is followed by an approximation of their accuracy. Considering the PS method’s theoretical accuracy limits of 1mm/y of the PS LOS velocities for both tracks as reported previously [39], an approximate accuracy of the vertical and ALD components, was estimated using the Equations (3) and (4). In the propagation of uncertainty law, the covariance is assumed to be zero.

For the error analysis, a mean value of the incidence angles was used, for the ascending and descending tracks, respectively, as well as the same mean heading values, already used in the decomposition analysis. It is also assumed that horizontal velocities are mainly due to east–west motion, owing to InSAR low sensitivity to the north component. For the error analysis, we used a mean value of 27.75° and 31.05° for the incidence angles of ascending and descending tracks, respectively, as well as the mean heading values, already used in the decomposition analysis. The estimated precision of the vertical and horizontal motion components is equal to 1.01 mm/y and 1.05 mm/y, respectively, and does not deviate a lot from the theoretical precision of the PS technique. The mean standard deviation thresholds of the vertical and horizontal displacements are also low and equal to 0.88 mm/y and 1.55 mm/y, respectively.

$${\sigma }_V^{Up}=\sqrt{{(\frac{{\theta }_V^{Up}}{{\theta }_V^{asc}})}^2}{\left({\sigma }_V^{asc}\right)}^2+{(\frac{{\theta }_V^{Up}}{{\theta }_V^{desc}})}^2{\left({\sigma }_V^{desc}\right)}^2$$

(3)

$${\sigma }_V^{ALD}=\sqrt{{(\frac{{\theta }_V^{ALD}}{{\theta }_V^{asc}})}^2}{\left({\sigma }_V^{asc}\right)}^2+{(\frac{{\theta }_V^{ALD}}{{\theta }_V^{desc}})}^2{\left({\sigma }_V^{desc}\right)}^2$$

(4)

where ${\sigma }_V^{Up }$is the precision of vertical velocities, ${\sigma }_V^{ALD}$ is the precision of ALD velocities, ${\sigma }_V^{asc}$ is the precision of LOS velocities of the ascending satellite pass, and ${\sigma }_V^{desc}$is the precision of LOS velocities of the descending satellite pass.

4. Results

The combination of PS and SBAS, implemented on data-sets from Envisat sensor, provided a significant increment of spatial sampling leading to an overall of 38,525 permanent scatterers, despite the small data sample. The PSI processing of Sentinel-1 data resulted in 33,007 scatterers for the ascending track and 48,425 for the descending track. The line-of-sight deformation rates estimated for Kythira Island, spanning from 2003 to 2019, are low.

4.1. The 2003–2009 Period

The mean LOS velocity map for Envisat data, covering the time span between 2003 and 2009, as obtained by the combination of PS and SBAS methods, is presented in Figure 3. The implementation of the SBAS method, significantly increased the spatial pixel sampling, considering the constraints in PSI processing, associated with limited data and the geomorphology of the island. A general observation is that almost 90% of the scatterers on the island show LOS deformation lower than 2 mm/y. The four islets situated to the east of Kythira island share the same deformation pattern with that at the eastern coast of Kythira and are characterized by low mean deformation rates.

Small-scale displacements dominate in the majority of Kythira Island, not exceeding −2 mm/y and 2 mm/y, while the mean standard deviation threshold is 1.04 mm/y. Yet, an interesting observation is the small-scale, localized displacements towards the satellite occurring in the center of the island and the eastern coast. Positive deformation rates are observed in the region around the Mitata village and northwest of the Livadi village, with a maximum velocity rate of 3 mm/y. Displacements with a maximum rate of 3 mm/y are observed in the coast to the east of KITH GPS station. Finally, an area which is characterized by a maximum deformation rate of 2 mm/y towards the satellite is the northeast coast of the Potamos village. This site also appears to be an uplifting area, during the 2015–2019 investigation time span, but with lower deformation rates (Figure 4).

4.2. The 2015–2019 Period

The LOS displacement maps (Figure 4) as obtained by the implementation of the StaMPS [31], for both ascending and descending imaging geometries, also denote small scale displacements in Kythira Island for the 2015–2019 period. The displacements assessed from Sentinel-1 interferograms (Figure 4) are consistent with the deformation pattern estimated from Envisat data processing (Figure 3). In all stacks from 2003 to 2019, Kythira Island is experiencing low rate deformation. Mean LOS velocity rates also extend from −2 mm/y to 2 mm/y, in most parts of the island and the four islets to the east of Kythira, with a mean standard deviation threshold of 0.8 mm/y for the ascending track and 0.5 mm/y for the descending track.

A general observation from the LOS displacements for the entire period from 2003 to 2019, is the positive velocity rates of PS scatterers in most parts of the island, with a significant increment in the central part of the island to the west of Mitata and Livadi villages (Figure 3 and Figure 4). Higher displacement rates with a maximum value of 5 mm/y are observed in the ascending track no. 102 of Sentinel-1 (Figure 4a), rather than in the descending track no. 7 (Figure 4b) and the ascending track no. 415 (Figure 3) of the Envisat period. To further investigate these observations, time-series plots are produced for a number of scatterers within a 200-m radius around the areas of interest (Figure 5) for all stacks, during both time periods (Figure 6). The dates in the horizontal axis represent the intreferogram dates and each scatterer corresponds to the mean deformation of all scatterers included in the 200-m radius selected.

In all time-series plots, information about the quality of the linear fit is provided. R-squared statistical measure, in Sentinel-1 time-series plots, is equal to 0.7 and 0.8, values that represent a large positive linear association and denote a good linear fit. Moreover, the root mean square error (RMSE) is added, as a measure of the absolute fit. In general, low values of RMSE characterize our time-series plots, indicating that the observed data points are close to the model’s predicted values. Finally, error bars were calculated and added to the plots, as a representation of the variability of our data. In the time-series plot generated for Envisat data (Figure 6e), RMSE is increased and R-squared is decreased, denoting a small positive linear association. We expect this to happen, especially when we compare 2003–2009 with 2015–2019 results. Larger perpendicular baselines in Envisat data, up to 1200 m, are a noise source for SAR interferometry, introducing DEM errors and spatial decorrelation phenomena.

The results verify the observation about the positive LOS displacements in the center of the island, based on the average LOS velocity plots. While, for the center of the island, the deformation is consistent in all stacks, a small shift in the direction of the deformation is observed between the LOS velocities of the ascending (no. 102) and the descending (no. 7) tracks of Sentinel-1 sensor. The coast to the east of the Potamos village (Figure 4 and Figure 7) is characterized by positive LOS displacements in the ascending track of Sentinel-1 data, while, in the descending track, it seems to follow an opposite trend. This alteration in motion very likely underlines the occurrence of horizontal movements along this coast. To further analyze this observation, time-series plots were also generated for scatterers included in a 200-m radius near the Potamos village (Figure 8).

4.3. Horizontal and Vertical Displacements

A common element in all deformation maps produced for Kythira Island during both time periods is the low displacement rate, dominating in most parts of the island. However, the decomposition analysis provided additional details on areas of special interest, initially identified in the LOS velocity maps.

The first outcome of the decomposition analysis is the vertical uplift dominating in the center and to the east of the island (Figure 9a). The second observation is related to the opposite direction of the LOS velocity vectors for the northeast part of Kythira Island, between the ascending satellite pass no. 102 and the descending track no. 7. By comparing the ascending and descending LOS displacements in Kythira island for the 2015–2019 time period (Figure 4), the initial assumption that a shift in the motion direction is associated with horizontal motion contribution is confirmed by the horizontal displacements map (Figure 9b). The coast to the northeast of the Potamos village experiences horizontal deformation with deformation rates with a maximum value of −3 mm/y. Considering the decomposition in the descending look direction, negative deformation rates correspond to a motion towards the East.

5. Discussion

The possible correlation between crustal deformation rates derived from geodesy and geology is a long-lasting issue with results showing usually divergent rates from the two independent approaches. For example, researchers [41] compared geodetic and geologic data from the Wasatch region, Utah, and discussed several alternatives to explain the mismatch between the various data sets at deformation periods from 10 to 10 million years. Comparison of pre-event geodetic and geologic rates associated with three large-magnitude strike-slip earthquakes (Izmit, North Anatolian Fault, 1999, M_W = 7.6; Kokoxili rupture on the Kunlun fault, Tibet, 2001, M_W = 7.8; Denali Fault, Alaska, 2002, M_W = 7.9) revealed a wide range of behaviors [42].

Tectonically, the island of Kythira is dominated by the presence of small-length Late Cenozoic faults striking mainly NNW-SSE [10,22,43]. Detailed mapping of the tectonic structures on the island of Kythira revealed a major extensional detachment fault in the northern part of the island (Figure 9b) named Potamos Detachment Fault [40], which is also shown in earlier tectonic mapping [22]. In addition, small strike-slip faults of E–W direction were found [40] on the east coast of the island, near the Platia Ammos site, as well as to the west of the Livadi village (Figure 9b). However, no active faults have been identified on the island and the geologic displacement rate remains poorly understood. In this sense, a comparison between geodetic and geologic rates in Kythira is not realistic.

On the other hand, the intense relief in the central part of the island is controlled by fault zones [22]. The island is experiencing a long-term constant uplift, which is more intense at the eastern coasts [22]. From 2015 to 2019, the geodetic uplift observed in the central part of Kythira Island appears in the vertical displacements map. Moreover, positive LOS displacements are identified in the 2003–2009 period. The displacements map (Figure 9a) obtained for the 2015–2019 time period, indicates vertical motion up to 5 mm/y in the central part of the island, although, in the area southwest of the Mylopotamos village, the deformation rate is lower. However, a pattern of lower and smooth LOS deformation is observed in the largest part of the island for the 2003–2009 time period (Figure 3). The positive deformation rate identified geodetically in this paper is about one order of magnitude larger than the average geologic long-term uplift rate estimated at ~0.13 mm/y for the Quaternary period [44]. It should be mentioned, however, that long-term geologic uplift is subject to important variability, e.g., in west Crete, the long-term uplift rate of 2.5–2.7 mm/y was found, but this average value includes periods in which the rate varied significantly up to 7.7 mm/y (e.g., [45]).

The main geodetic deformation pattern dominating a large part of the island is uplift, ranging from 2 to 5 mm/y. We suggest that the current deformation pattern identified in Kythira Island from space geodesy is rather explainable by regional deformation processes rather than by deformation in local small-scale, no active faults on the island. Earthquake-cycle modeling is not helpful in this case, since the historical earthquake record of the area is poorly constrained, while the instrumental seismicity is dominated by the two intermediate-depth earthquakes of 1903 and 2006. However, the contribution of post-seismic transients to the deformation pattern and rates is a possibility, e.g., the Early Byzantine tectonic paroxysm in the western Hellenic arc [46]. Aseismic slip is an alternative mechanism for uplift in plate margins. In the Hellenic Subduction Zone, the seismic slip rate represents only a small fraction of plate motion, the largest part up to ~80% being accommodated by aseismic slip (e.g., [3,47]). An additional uplift mechanism that should not be ruled out is sediment underplating. Whatever the causes of the different geologic long-term uplift and geodetic short-term uplift, more research is needed for better understanding the long-term variations of the geological displacement rates in the island of Kythira.

6. Conclusions

The location of Kythira Island in the very active tectonic zone of the western Hellenic arc, combined with the damaging earthquake (M_w = 6.7) of 8 January 2006, stimulated the investigation of ground deformation on the island with the use of SAR interferometry from Envisat and Sentinel-1 satellite data covering the time span from 2003 to 2019. For the 2015–2019 time period, the combination of slant range velocities using both ascending and descending satellite tracks contributed to the estimation of vertical and horizontal displacements.

Horizontal displacements are identified in the coast northeast of the Potamos village with a maximum of −3 mm/y, corresponding to an eastward motion trend. The same type of motion is observed in the region to the west of the Livadi village, with lower deformation rates. Strike-slip components, identified in past tectonic studies [40], were verified and highlighted by horizontal motion trends, derived by the decomposition of SAR displacements estimated to the LOS direction.

The vertical displacement field indicates an upward motion trend in the central part of the island with a maximum rate of 5 mm/y for the 2015–2019 time period. Additionally, lower positive rates are identified SW of the Mylopotamos village. LOS displacements of lower rates, on the order of 2–3 mm/y, characterize the deformation field in the 2003–2009 time period. These uplift rates characterizing the ongoing deformation process in the island is one order of magnitude larger than the average long-term geologic uplift rate of ~0.13 mm/y found for the Quaternary period [44]. The discrepancy between the higher short-term geodetic uplift and the lower long-term geologic uplift could be explained by the variation of the latter, due to several regional deformation mechanisms including post-seismic transients, aseismic slip in the HSZ plate margin, and sediment underplating.