Multi-Sensor SAR Geodetic Imaging and Modelling of Santorini Volcano Post-Unrest Response

Volcanic history of Santorini over recent years records a seismo-volcanic unrest in 2011–12 with a non-eruptive behavior. The volcano deformation state following the unrest was investigated through multi-sensor Synthetic Aperture Radar Interferometry (InSAR) time series. We focused on the analysis of Copernicus Sentinel-1, Radarsat-2 and TerraSAR-X Multi-temporal SAR Interferometric (MT-InSAR) results, for the post-unrest period 2012–17. Data from multiple Sentinel-1 tracks and acquisition geometries were used to constrain the E-W and vertical components of the deformation field along with their evolution in time. The interpretation of the InSAR observations and modelling provided insights on the post-unrest deformation pattern of the volcano, allowing the further re-evaluation of the unrest event. The increase of subsidence rates on Nea Kameni, in accordance with the observed change of the spatial deformation pattern, compared to the pre-unrest period, suggests the superimposition of various deformation sources. Best-fitting inversion results indicate two deflation sources located at southwestern Nea Kameni at 1 km depth, and in the northern intra-caldera area at 2 km depth. A northern sill-like source interprets the post-unrest deflation attributed to the passive degassing of the magma intruded at 4 km during the unrest, while an isotropic source at Nea Kameni simulates a prevailing subsidence occurring since the pre-unrest period (1992–2010).


Introduction
Santorini volcano is considered to belong in caldera-forming systems undergoing long-term periods of quiescence, of approximately 20,000 years ( Figure 1). Although its last big eruption dates back 3600 years, its volcanic activity up to the most recent eruption in 1950 [1] was intertwined with the building of the intra-caldera islets of Palea and Nea Kameni. The latest volcano's reactivation was followed by the restless period of 2011, which did not culminate in an eruption. Increased microseismic activity [2] and significant ground uplift [3] reaching 14 cm from March 2011 to March 2012 at Cape Skaros (Figure 1), and 9 cm at Nea Kameni islet [4], underlined the seismo-volcanic unrest. Most studies based on geodetic data interpret the episode with a single inflation source due to magma intrusion, at the northern part of the caldera, estimated within 3-4 km depth [4][5][6][7][8]. An alternative model based on GPS data [9], proposes two inflationary magmatic sources that relocate in depth and geographic location throughout the unrest. The evolution of the shallow magma body beneath Recent advances in geodetic imaging techniques and the systematic availability of Synthetic Aperture Radar (SAR) data from spaceborne missions provided the necessary tools and data, in order to monitor the deformation field of Santorini volcano over the post-unrest period, and to re-evaluate the deformation mechanism during the unrest. MT-InSAR techniques were employed to generate and analyze the SAR deformation time series. A comprehensive analysis of multi-sensor SAR acquisitions of Copernicus Sentinel-1, Radarsat-2 and TerraSAR-X satellite missions was performed to investigate the evolution of ground deformation from 2012 to 2017. The Copernicus Sentinel-1 and Radarsat-2 measurements were further analyzed in an inversion framework to estimate the source parameters that can best resolve the observed displacements, while TerraSAR-X data provided an additional source of information to verify the deformation pattern. Finally, insights on the post-unrest response of the volcano and the interrelationship between the post-unrest interval and the unrest event were investigated, while ERS-1 and -2 and ENVISAT data from 1992 to 2010 [3,4] were additionally used to interpret the similarity between the pre-and post-unrest volcano's deformation.

Santorini Volcanic Setting
Santorini volcano is part of the Hellenic Volcanic Arc in the Southern Aegean Sea (Figure 1), and it is partly situated on a SW-NE trending tectonic horst, the Amorgos Ridge, whereas the northwestern Remote Sens. 2019, 11, 259 3 of 18 half of the volcanic field lies within the Anydros Basin [15]. Extensional tectonics seems to have a profound effect on Santorini volcano. Regional fault systems, namely the NE-SW striking Kameni and Columbos Fault Zones (KFZ, CFZ) [16][17][18] mark the alignment of several eruptive vents [16,19] and have been interpreted as major normal faults moving in response to a NW-SE extension [20].
Volcanism in Santorini evolved with time from an early formation of volcanic cones, two successive explosive cycles of pyroclastic eruptions, and the later collapse of the caldera with the development of shield volcanoes within the caldera [16]. Since the caldera-forming Minoan eruption in the late 1600s BC, a peculiar volcanic configuration was finally set, with two active volcanoes, the Nea Kameni volcano located at the center of the caldera and the Columbos submarine volcano located almost 7 km NE of Thera [21][22][23] (Figure 1). The largest volcanic eruptions date to 197 BC, 1866, 1925 and 1949-50 and are all associated with Nea Kameni volcano. However, an eruption in 1650 AD took place offshore at Columbos volcano [24,25]. Aseismic and small-scale intrusions are inferred to have been emplaced on the northern caldera between 1994 and 1999 [26]. Marine geophysical surveys in the same area provide evidence of shallow intrusions of the post-Minoan activity [27]. The most recent activity of the volcano was recorded at the beginning of 2011 lasting up to the first half of 2012 [3][4][5][6][7][8][9]. This activity was marked by an increase of low magnitude (ML < 3.2) earthquakes, mainly aligned along the Kameni fault [28], deformation and changes in the geochemical parameters of the gas emissions. These unrest signals were concentrated within the caldera sector surrounding Nea Kameni and Thera [2,29]. Since the ending of the unrest, seismicity associated with caldera volcanic activity has decreased ( Figure 1).

SAR Interferometric Processing
Differential Interferometric SAR (DInSAR) techniques are widely used to detect and monitor subtle ground displacements associated with volcanoes, such as dome growth and subsidence [30,31]. By combining SAR observations from different sensors and viewing geometries, over common acquisition periods, a proper characterization of the spatial and temporal deformation patterns can be robustly achieved. Advanced Multi-Temporal InSAR (MT-InSAR) time series analysis techniques offer the capability to monitor the temporal evolution of ground deformation phenomena. These techniques lie on the capability to utilize large series of SAR imagery to measure small deformation signals, to millimeter level of accuracy, on individual Persistent Scatterers (PS) point targets (human infrastructures, natural scatterers). Several MT-InSAR techniques have been proposed (e.g., References [32][33][34][35][36][37][38][39][40][41][42][43][44][45]), each following different approaches for the identification of point scatterers and resolving displacement histories. A detailed review of proposed MT-InSAR techniques in the literature is presented in Reference [46].
In this study, the MT-InSAR approach proposed by Reference [47] was adopted, by using the GAMMA s/w packages. The implemented Terrain Observation by Progressive Scans (TOPS) acquisition mode [48] on Sentinel-1 mission required specific/particular interferometric handling to ensure proper coregistration of burst-type data compared to standard stripmap acquisitions [49]. Sentinel-1 captures 250 km in three sub-swaths in range, each divided into nine bursts in azimuth. Burst synchronization ensures the interferometric capability of the mission. Since accuracy at 0.005 pixels is required [50], coregistration is performed on a pixel basis using initial orbital information and DEM-assisted cross-correlation methods, while an iterative refinement procedure is followed, by using the Enhanced Spectral Diversity (ESD) technique [49,51]. The coregistration scheme was further adapted to avoid temporal decorrelation effect on ESD estimates, especially since the burst overlap areas for Santorini are spatially limited over land. The estimation of ESD offsets for each scene were based on the temporally closest already coregistered slave, in such a way that coregistration proceeds successively, starting from the scenes closest to the master to the most temporally distant ones.
Interferometric processing with multi-looking factors of 6x2 in range and azimuth, respectively, was considered in order to mitigate signal variability and to obtain comparable pixel spacing to previous ERS-ENVISAT results. Topographic phase was simulated and subtracted based on a 20 m resolution Digital Elevation Model (DEM) generated from 1/50.000 scale topographic map of the Hellenic Military Geographical Service (HGMS).
Since no variability of deformation within a 6-day interval was expected, only acquisitions from Sentinel-1A satellite were considered (12-days repeat cycle), offering a sufficiently dense temporal stack, to minimize unwrapping artifacts due to decorrelation effects, and characterize the spatio-temporal behavior of the deformation signal. More specifically, three Sentinel-1 orbit geometries were used (Table 1), covering the period from October 2014 up to December 2017. Data used for the ascending relative orbit A029 amounts to 91 acquisitions, while for the descending relative orbits D109 and D036 correspond to 93 and 92 acquisitions, respectively ( Figure 2). Since no variability of deformation within a 6-day interval was expected, only acquisitions from Sentinel-1A satellite were considered (12-days repeat cycle), offering a sufficiently dense temporal stack, to minimize unwrapping artifacts due to decorrelation effects, and characterize the spatiotemporal behavior of the deformation signal. More specifically, three Sentinel-1 orbit geometries were used (Table 1), covering the period from October 2014 up to December 2017. Data used for the ascending relative orbit A029 amounts to 91 acquisitions, while for the descending relative orbits D109 and D036 correspond to 93 and 92 acquisitions, respectively ( Figure 2).   A key point of the processing, given the size of the data stack, was to derive a redundancy of all possible interferometric pairs. The final selection of the interferograms was done by limiting the temporal baselines, as well as by using an upper threshold on the maximum normal baseline value. Although the orbital tube of Sentinel-1 is well-tuned for interferometric applications (within 120 m radius), the presence of steep slopes along the caldera walls, of height differences up to approx. 300 m, led to the further control of the baselines. Interferometric stacks of 160, 160 and 242 pairs were generated for the A029, D109 and D036 orbits accordingly, for the final MT-InSAR configuration ( Table 2). The Singular Value Decomposition (SVD) approach was used to obtain a solution for the phase time series, based on a set of multi-reference point differential interferograms. Deformation phase time series were estimated using a weighted least-squares algorithm that minimized the sum of squared weighted residual phases [36]. The phase model included a height related term proportional to the derivative of the interferometric phase with respect to height, allowing the update of PS heights during the interferometric processing. To mitigate the atmospheric contribution, filtering based on assumptions of atmospheric statistics [33] was applied in both temporal and spatial domains, in a way to minimize the effect of spatially correlated and temporally uncorrelated signal.
Redundancy of interferometric pairs reduces uncertainties in the time series mainly due to phase errors introduced by decorrelation and residual topography. With respect to the selected pairs ( Figure S1), few scenes per orbit were excluded from the analysis; still temporal sampling was sufficiently dense to examine temporal variability in PS displacement histories. For each solution the linear displacement phase rate, the standard deviation of the phase time series relative to the linear fit, as well as the unwrapped phases for each PS target were stored.

Hosted Processing on Geohazards Exploitation Platform
The Geohazards Lab initiative, originated by the European Space Agency (ESA) with the support of several other space agencies of the Committee on Earth Observation Satellites (CEOS), is based on a group of interoperable platforms with federated resources providing Earth Observation (EO) data access, hosted processing and e-collaboration capabilities to animate and support the geohazards user community [52]. One of its precursors is the Geohazards Exploitation Platform (GEP) (https: //geohazards-tep.eu), developed in the framework of the ESA Thematic Exploitation Platforms (TEP) initiative and has been available since 2016. Today, GEP has primary focus on mapping hazard-prone land surfaces and monitoring terrain deformation. The platform has been expanded to include a broad range of products and services, currently available or under development on cloud processing resources like the GEP, to support experts and users to better understand geohazards and relevant risks.
In the framework of the ESA funded project "Disaster Risk Reduction using innovative data exploitation methods and space assets", ground deformation maps were produced on GEP for the Santorini volcano. Processing was based on a customized implementation of the StaMPS Persistent Scatterer Interferometry (PSI) software [40]. For the period 2012-16 average Line-Of-Sight (LOS) velocities were based on 20 descending Radarsat-2 images, while TerraSAR-X data (25 scenes in descending orbit) were used for the 2012-13 period (Table 1). Datasets are available via Zenodo repository [53].
GEP PSI results are provided as raw velocities generated from automatic processors and have not undergone any post-processing. The adjustment of the reference point was addressed to ensure compatibility with Sentinel-1 results obtained herein, and visual inspections were performed to exclude regions potentially affected by isolated unwrapping errors.

Interferometric Results
A multi-temporal interferometric analysis approach was implemented to explore the dense temporal series of Sentinel-1 data, and to provide the displacement time series for each identified PS target. The MT-InSAR results are shown in Figure 3, presenting the measured LOS displacement rates for each of the Sentinel-1 acquisition geometries. Common deformation patterns are retrieved, indicating LOS motion rates up to −7 and −9 mm/yr over Nea Kameni islet. Variations of maximum observed motion could be partly attributed to the different incident angles between the acquisition geometries. Despite the large amount of data, the relatively short observation period (3.2 years) does not permit substantial decrease of the measurement uncertainties, though they were maintained at a level of 1.5-2.0 mm/yr ( Figure S2). Additionally, the relatively high coherence levels, due to the short span of the interferometric pairs used, especially over Nea Kameni, underline the robustness of the obtained results. It is worth mentioning that deformation does not seem to follow exactly the same pattern among the different orbits ( Figure 3). Therefore, the southwestern part of Nea Kameni seems to be most affected in the ascending orbit, while in the descending ones, maximum deformation is shifted to the central and northern parts of the islet. This spatial variability implies the presence of horizontal motion, also mentioned in previous study [4].
The displacement time series for a PS target over Nea Kameni islet is shown in Figure 4. It is possible to observe the presence of a periodic (annual) component of the displacement, identical both in terms of period and amplitude for all Sentinel-1 datasets ( Figure 5). The amplitudes of the oscillations (up to 5 mm) are not negligible with respect to the expected linear deformations rates. Such cyclic motions can therefore significantly affect, by under-or over-estimation, the deformation rates calculated based on linear regression of the displacement time series over a few years' timespans. To compensate for this source of error, a linear interpolation to the time series was initially applied to produce a regular 6-day sampling, in order to fill gaps in the temporal series. To quantify quasi-annual motion, we then calculated a Fourier transform of the time series by only selecting periods between 200 and 500 days. Other frequency intervals were set to zero. The maximum amplitude of the resulting spectrum was then identified. The corresponding cyclic function (d) of time (t) was retrieved on the basis of the amplitude (A), period (T) and phase ( ) of the Fourier transform for the frequency of the spectrum's maximum (Equation (1)): It is worth mentioning that deformation does not seem to follow exactly the same pattern among the different orbits ( Figure 3). Therefore, the southwestern part of Nea Kameni seems to be most affected in the ascending orbit, while in the descending ones, maximum deformation is shifted to the central and northern parts of the islet. This spatial variability implies the presence of horizontal motion, also mentioned in previous study [4].
The displacement time series for a PS target over Nea Kameni islet is shown in Figure 4. It is possible to observe the presence of a periodic (annual) component of the displacement, identical both in terms of period and amplitude for all Sentinel-1 datasets ( Figure 5). The amplitudes of the oscillations (up to 5 mm) are not negligible with respect to the expected linear deformations rates. Such cyclic motions can therefore significantly affect, by under-or over-estimation, the deformation rates calculated based on linear regression of the displacement time series over a few years' time-spans. It is worth mentioning that deformation does not seem to follow exactly the same pattern among the different orbits ( Figure 3). Therefore, the southwestern part of Nea Kameni seems to be most affected in the ascending orbit, while in the descending ones, maximum deformation is shifted to the central and northern parts of the islet. This spatial variability implies the presence of horizontal motion, also mentioned in previous study [4].
The displacement time series for a PS target over Nea Kameni islet is shown in Figure 4. It is possible to observe the presence of a periodic (annual) component of the displacement, identical both in terms of period and amplitude for all Sentinel-1 datasets ( Figure 5). The amplitudes of the oscillations (up to 5 mm) are not negligible with respect to the expected linear deformations rates. Such cyclic motions can therefore significantly affect, by under-or over-estimation, the deformation rates calculated based on linear regression of the displacement time series over a few years' timespans. To compensate for this source of error, a linear interpolation to the time series was initially applied to produce a regular 6-day sampling, in order to fill gaps in the temporal series. To quantify quasi-annual motion, we then calculated a Fourier transform of the time series by only selecting periods between 200 and 500 days. Other frequency intervals were set to zero. The maximum amplitude of the resulting spectrum was then identified. The corresponding cyclic function (d) of time (t) was retrieved on the basis of the amplitude (A), period (T) and phase ( ) of the Fourier transform for the frequency of the spectrum's maximum (Equation (1)): The obtained cyclic component was finally removed from the initial time series and the LOS displacement rates for each PS were re-computed as the slopes of the corrected displacement time series ( Figure 5). Such periodic signal can be attributed to the topography-dependent atmospheric component of the SAR, also reported for several other volcanoes [54,55]. By compensating for the seasonal signal in the time series, a noteworthy decrease in the displacement rates is apparent for the major part of the volcano ( Figure 6). The new estimation of the displacement rates provided the actual deformation that was used to model the source dynamics. It was not in fact possible for such correction to be derived from Radarsat-2 PSI results, as data provided via GEP platform [53] concerned only LOS velocities and not actual time series.  To compensate for this source of error, a linear interpolation to the time series was initially applied to produce a regular 6-day sampling, in order to fill gaps in the temporal series. To quantify quasi-annual motion, we then calculated a Fourier transform of the time series by only selecting periods between 200 and 500 days. Other frequency intervals were set to zero. The maximum amplitude of the resulting spectrum was then identified. The corresponding cyclic function (d) of time (t) was retrieved on the basis of the amplitude (A), period (T) and phase (ϕ) of the Fourier transform for the frequency of the spectrum's maximum (Equation (1)): The obtained cyclic component was finally removed from the initial time series and the LOS displacement rates for each PS were re-computed as the slopes of the corrected displacement time series ( Figure 5). Such periodic signal can be attributed to the topography-dependent atmospheric component of the SAR, also reported for several other volcanoes [54,55].
By compensating for the seasonal signal in the time series, a noteworthy decrease in the displacement rates is apparent for the major part of the volcano ( Figure 6). The new estimation of the displacement rates provided the actual deformation that was used to model the source dynamics. It was not in fact possible for such correction to be derived from Radarsat-2 PSI results, as data provided via GEP platform [53] concerned only LOS velocities and not actual time series.
Furthermore, LOS measurements were decomposed into vertical and E-W motion components to facilitate the better interpretation of motion patterns. The E-W and vertical components were calculated using all three Sentinel-1 geometries in a comprehensive decomposition scheme, considering as well the variability of the incidence angles within the scenes (Figure 7). The procedure involved rasterization of MT-InSAR point vectors and calculation of the two components on common pixels. Then, spline interpolation was applied to fill gaps and obtain a more visually appealing result. It is worth noting that, for the comparison against modelling results and the calculation of misfits, only pixels presenting real and not interpolated data were considered.
By compensating for the seasonal signal in the time series, a noteworthy decrease in the displacement rates is apparent for the major part of the volcano ( Figure 6). The new estimation of the displacement rates provided the actual deformation that was used to model the source dynamics. It was not in fact possible for such correction to be derived from Radarsat-2 PSI results, as data provided via GEP platform [53] concerned only LOS velocities and not actual time series.  Furthermore, LOS measurements were decomposed into vertical and E-W motion components to facilitate the better interpretation of motion patterns. The E-W and vertical components were calculated using all three Sentinel-1 geometries in a comprehensive decomposition scheme, considering as well the variability of the incidence angles within the scenes (Figure 7). The procedure involved rasterization of MT-InSAR point vectors and calculation of the two components on common pixels. Then, spline interpolation was applied to fill gaps and obtain a more visually appealing result. It is worth noting that, for the comparison against modelling results and the calculation of misfits, only pixels presenting real and not interpolated data were considered.
Prevalent subsidence for the entire Nea Kameni is indicated by the vertical motion map of Figure  7, reaching -7 mm/yr. The observed concentric deformation pattern presents an opening towards northern parts of the islet, where high deformation rates are also evident. Motion for the rest of the area has also been reduced after compensating for the seasonal signal, with overall rates (a couple of mm/yr), mostly within the uncertainties of the measurements. Concerning horizontal E-W motion, as expected, the N-S extending zero deformation zone on Nea Kameni is well aligned to the local concentric subsidence pattern, showing opposite motions on both sides and towards the area of maximum deformation. For the rest of the volcano, of significant interest is the relatively high westward motion of Cape Skaros, an area that accumulated the highest inflation rates during the 2011-12 unrest [4]. On the basis of the SAR Interferometry results, Nea Kameni remains the area Prevalent subsidence for the entire Nea Kameni is indicated by the vertical motion map of Figure 7, reaching −7 mm/yr. The observed concentric deformation pattern presents an opening towards northern parts of the islet, where high deformation rates are also evident. Motion for the rest of the area has also been reduced after compensating for the seasonal signal, with overall rates (a couple of mm/yr), mostly within the uncertainties of the measurements. Concerning horizontal E-W motion, as expected, the N-S extending zero deformation zone on Nea Kameni is well aligned to the local concentric subsidence pattern, showing opposite motions on both sides and towards the area of maximum deformation. For the rest of the volcano, of significant interest is the relatively high westward motion of Cape Skaros, an area that accumulated the highest inflation rates during the 2011-12 unrest [4]. On the basis of the SAR Interferometry results, Nea Kameni remains the area exhibiting the highest ground deformation during the post-unrest period, with higher rates compared to the pre-unrest period (average subsidence rate of −5 ± 0.6 mm/yr for the 1992-2010 period) [3,4,56]. In addition to this, the ongoing deformation at Cape Skaros on the main island indicates that the volcano has not fully recovered yet.
Finally, regarding GEP PSI results and specifically those of Radarsat-2 (2012-16), having major overlap with the Sentinel-1 observations (2014-17), they are showing comparable motion in terms of both spatial pattern and magnitude (Figure 8a). Providing a complementary LOS measurement angle and sufficiently large time span for obtaining a robust solution, they were considered in the source modelling. On the contrary, for TerraSAR-X PSI (Figure 8b), and mainly due to the difference in the temporal coverage (2012-13), as well as the level of noise in the data, a decision was made not to include them in the inversion modelling. Nonetheless, they provided valuable information to underline the higher deflation rates observed during the early post-unrest phase. LOS displacement rates up to -12 mm/yr were observed in the center of the Nea Kameni islet with a visible concentric deformation pattern, whereas a similar localized deformation is denoted at the northern part of the islet. Apart from the higher deformation of TerraSAR-X data, due to the temporal proximity to the unrest, spatial deformation patterns are consistent with Sentinel-1 data, indicating the persistence of the deformation maxima areas since 2012.

Source Modelling
InSAR data show subsidence of the inner islets and, partially in the main island of Thera, across the inner walls ( Figure 6). This pattern is quite new with respect to the pre-unrest, since 1992-2010 InSAR results show a subsidence pattern limited to the Kameni islets, and negligible deformation for the rest of the caldera [3,4,56,57]. Furthermore, the subsidence on the Kameni islets after the unrest extends to a wider area with respect to the pre-unrest period. During the unrest, higher uplift was concentrated on the main island at Cape Skaros, whereas Nea Kameni demonstrated comparable uplift rates of the northern part, with an observed tilting towards the south.
SAR modelling of a single deflating source attributed to the caldera-wide deformation yielded residuals in Nea Kameni justifying the use of a second source for the inversion ( Figure S3). In order to take into account the localized subsidence at Nea Kameni and the caldera-wide deformation, two separated sources of deformation were considered. A joint inversion of mean velocity data was performed from ascending and descending Sentinel-1, together with descending Radarsat-2. The datasets were pre-processed by masking the areas undergoing layover issues, affecting in particular the inner cliffs in the ascending orbit of Sentinel-1 data. Afterwards, the four datasets were subsampled with a step of 60 m in the inner islets and 200-300 m in the remaining islands of Santorini.

Source Modelling
InSAR data show subsidence of the inner islets and, partially in the main island of Thera, across the inner walls ( Figure 6). This pattern is quite new with respect to the pre-unrest, since 1992-2010 InSAR results show a subsidence pattern limited to the Kameni islets, and negligible deformation for the rest of the caldera [3,4,56,57]. Furthermore, the subsidence on the Kameni islets after the unrest extends to a wider area with respect to the pre-unrest period. During the unrest, higher uplift was concentrated on the main island at Cape Skaros, whereas Nea Kameni demonstrated comparable uplift rates of the northern part, with an observed tilting towards the south.
SAR modelling of a single deflating source attributed to the caldera-wide deformation yielded residuals in Nea Kameni justifying the use of a second source for the inversion ( Figure S3). In order to take into account the localized subsidence at Nea Kameni and the caldera-wide deformation, two separated sources of deformation were considered. A joint inversion of mean velocity data was performed from ascending and descending Sentinel-1, together with descending Radarsat-2. The datasets were pre-processed by masking the areas undergoing layover issues, affecting in particular the inner cliffs in the ascending orbit of Sentinel-1 data. Afterwards, the four datasets were subsampled with a step of 60 m in the inner islets and 200-300 m in the remaining islands of Santorini.
The non-linear geodetic inversion was performed considering a 2-steps procedure based on the Neighbourhood Algorithm [59,60]. The VSM (Volcano and Seismic source Model) is a FORTRAN code retrieving not only the best-fit model, but also performing a second step with a Bayesian inference in the generated model ensemble. For the reasons described above, two sources were adopted, one for the caldera-wide deformation, and another one for the localized subsidence at Kameni islets. After several attempts aiming to minimize the residuals (a chi-square function) between data and models, we obtained a chi-square equal to 1.2, and the preferred source models are to be a sill-like source [58] for the caldera-wide deformation, and an isotropic source for the localized deformation in Nea Kameni [61,62].
Results are shown in Figures 9 and 10 by combining Sentinel-1 data for the East-West and vertical components. The modelled data is simply the East-West and vertical component of the mean model, without further interpolation. Detailed comparisons for all the datapoints jointly inverted are reported in Figure S4, while the posterior probability density functions are shown in Figures S5 and S6. The caldera-wide source is located below the sea, north of Kameni, at a depth of 2 km. It has a radius of a few hundred meters and underwent a volume variation rate of −12·10 4 m 3 yr −1 during 2012-17 ( Table 3). The volume variation is not inverted, but it is computed from the sill parameters, considering the medium as Poissonian. The second source is located below the area of maximum subsidence at Nea Kameni, at a depth of about 1 km, and characterized by a volume variation rate of −1.4·10 4 m 3 yr −1 . The surface projection of the source center coincides with the change of direction of the East-West component, as expected for isotropic sources. The residuals show a trend of underestimation in the northern islands of Santorini for the horizontal component. However, the horizontal displacement is affected by larger errors than the vertical, amounting to few mm/yr. The vertical component shows some residuals across the northern cliffs, but as that area was affected by layover issues in the ascending orbit, it is therefore characterized by a larger uncertainty than the rest of the dataset.  Table 3. Data comparisons concern East-West and vertical components. In the first column are the East-West The non-linear geodetic inversion was performed considering a 2-steps procedure based on the Neighbourhood Algorithm [59,60]. The VSM (Volcano and Seismic source Model) is a FORTRAN code retrieving not only the best-fit model, but also performing a second step with a Bayesian inference in the generated model ensemble. For the reasons described above, two sources were  Table 3. Data comparisons concern East-West and vertical components. In the first column are the East-West

Discussion
The temporal evolution of ground deformation at Santorini, after the unrest episode, seems to be controlled by a set of two shallow sources, favoring a more complex scenario of the volcano dynamics ( Figure 11). The model proposed here, probed with four different C-band SAR observations from Sentinel-1 and Radarsat-2 missions, led to better constrain the deformation sources acting during the post-unrest period. The two-source model considers two deflating sources with a 4 km distance in between. The first source is a deflating, Mogi-like source located at about 1 km depth and centered below Kameni island. The second source is a sill-like body at 2 km depth located just above the deeper (about 4 km) source of the 2011-12 unrest episode [4][5][6][7][8][9]. According to the available geochemical data, an increase of CO2 concentration has been observed after the unrest period [63]. This increase has been interpreted as an increase in the permeability of the volcanic pile along the Kameni fault resulting from the seismic activity during the unrest. According to these data, and taking into account the results of recent models of degassing from magma chambers [64], we propose that dynamics of our sill-like deflating source is related to the passive degassing (second boiling) of the magma intruded at 4-5 km depth during the unrest. This interpretation is consistent with the shallower depth (about 2 km) and geometry (sill-like) of the top of the magma reservoir, in which a volatile saturated zone is formed at the top by accumulation of gases associated to initial states of fractional crystallization and cooling [64][65][66]. This interpretation is also supported by the available geochemical data, which show a significant increase of CO2 concentration values from the unrest, uplift period (May 2010-February 2012, CO2 = 400 mmol/mol) to the post-unrest (deflation) period (March-July 2012; CO2 = 800 mmol/mol) [62]. Such increase excludes the occurrence of drain-back of magma because, in that case, a decrease in the CO2 concentration values should be observed.

Discussion
The temporal evolution of ground deformation at Santorini, after the unrest episode, seems to be controlled by a set of two shallow sources, favoring a more complex scenario of the volcano dynamics ( Figure 11). The model proposed here, probed with four different C-band SAR observations from Sentinel-1 and Radarsat-2 missions, led to better constrain the deformation sources acting during the post-unrest period. The two-source model considers two deflating sources with a 4 km distance in between. The first source is a deflating, Mogi-like source located at about 1 km depth and centered below Kameni island. The second source is a sill-like body at 2 km depth located just above the deeper (about 4 km) source of the 2011-12 unrest episode [4][5][6][7][8][9]. According to the available geochemical data, an increase of CO 2 concentration has been observed after the unrest period [63]. This increase has been interpreted as an increase in the permeability of the volcanic pile along the Kameni fault resulting from the seismic activity during the unrest. According to these data, and taking into account the results of recent models of degassing from magma chambers [64], we propose that dynamics of our sill-like deflating source is related to the passive degassing (second boiling) of the magma intruded at 4-5 km depth during the unrest. This interpretation is consistent with the shallower depth (about 2 km) and geometry (sill-like) of the top of the magma reservoir, in which a volatile saturated zone is formed at the top by accumulation of gases associated to initial states of fractional crystallization and cooling [64][65][66]. This interpretation is also supported by the available geochemical data, which show a significant increase of CO 2 concentration values from the unrest, uplift period (May 2010-February 2012, CO 2 = 400 mmol/mol) to the post-unrest (deflation) period (March-July 2012; CO 2 = 800 mmol/mol) [62]. Such increase excludes the occurrence of drain-back of magma because, in that case, a decrease in the CO 2 concentration values should be observed.
Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 18 Figure 11. Schematic diagram summarizing SAR geodetic results presented herein and in Reference [4]. A shallow caldera-wide source (sill) at around 2 km depth reflects the post-unrest deflation of the deeper magma chamber responsible for the 2011-12 unrest episode. The volumetric spherical/spheroidal unrest source currently seems to undergo a depressurization due to degassing from its upper parts. A shallower deflating source (Mogi), at around 1 km depth, is also shown below Nea Kameni, exhibiting similar behavior for both the 1992-2010 pre-unrest and 2012-17 post-unrest periods.
Therefore, the degassing from this shallower source, which easily develops along the preferred pathway of the Kameni fault zone (e.g., References [28,63]) is responsible for the observed deflation of the 2 km deep sill-like source. The volume variation involved, as computed from the elastic models is in the order of 10 5 m 3 /yr, which is 1/100 of the volume variation involved during the 2011-12 unrest. Regarding the nature of the degassing magma, Reference [67] suggest that a 3 He-and volatile-rich mafic magma intruded during the unrest. According to the conceptual model proposed above, such volatiles migrated in the uppermost portion of the reservoir and were passively released to the surface along the Kameni line. A component of deformation related to magma cooling can also be invoked for the observed subsidence, however, according to Reference [64], it plays a minor role.
Prior studies based on GPS and SAR observations for the unrest period [4][5][6][7][8][9], yielded a single inflating magma source in the northern intra-caldera area. While nearly all solutions follow almost N-S dispersion, they still indicate the same unrest center, if we take into consideration their location uncertainties ( Figure 10). On the contrary, Reference [9] using GPS time series determined two source locations, one offshore, about 1 km northern of the N-S cluster solutions, with comparable depth, and a deeper one (7-8 km) on Nea Kameni islet. Our post-unrest deflating sill at the northern intra-caldera area is in good agreement with the unrest Mogi source of Reference [9], for their better constrained solution (period P1, Jan-Jun 2011), in terms of number of GPS stations and larger deformation gradients. On the whole, the inefficacy to well constrain magma sources using single SAR LOS observations, or GPS measurements that suffer from poor spatial sampling should be underlined. The utilization of multiple geodetic assets is thus beneficial and should be certainly taken into account when designing volcanic monitoring systems.
The InSAR results and modelling for the Nea Kameni deformation source of this study demonstrate a comparable pattern to previous ERS and ENVISAT observations covering the 1992-2010 pre-unrest period [3-4,56] (Figure 12a). The source's identical location and volume variation rate between the two intervals points to the hypothesis of a steady deformation signal. By assuming the pre-unrest deformation rate (ERS/ENVISAT) as background motion, and subtracting it from the Figure 11. Schematic diagram summarizing SAR geodetic results presented herein and in Reference [4]. A shallow caldera-wide source (sill) at around 2 km depth reflects the post-unrest deflation of the deeper magma chamber responsible for the 2011-12 unrest episode. The volumetric spherical/spheroidal unrest source currently seems to undergo a depressurization due to degassing from its upper parts. A shallower deflating source (Mogi), at around 1 km depth, is also shown below Nea Kameni, exhibiting similar behavior for both the 1992-2010 pre-unrest and 2012-17 post-unrest periods. Therefore, the degassing from this shallower source, which easily develops along the preferred pathway of the Kameni fault zone (e.g., References [28,63]) is responsible for the observed deflation of the 2 km deep sill-like source. The volume variation involved, as computed from the elastic models is in the order of 10 5 m 3 /yr, which is 1/100 of the volume variation involved during the 2011-12 unrest. Regarding the nature of the degassing magma, Reference [67] suggest that a 3 He-and volatile-rich mafic magma intruded during the unrest. According to the conceptual model proposed above, such volatiles migrated in the uppermost portion of the reservoir and were passively released to the surface along the Kameni line. A component of deformation related to magma cooling can also be invoked for the observed subsidence, however, according to Reference [64], it plays a minor role.
Prior studies based on GPS and SAR observations for the unrest period [4][5][6][7][8][9], yielded a single inflating magma source in the northern intra-caldera area. While nearly all solutions follow almost N-S dispersion, they still indicate the same unrest center, if we take into consideration their location uncertainties ( Figure 10). On the contrary, Reference [9] using GPS time series determined two source locations, one offshore, about 1 km northern of the N-S cluster solutions, with comparable depth, and a deeper one (7-8 km) on Nea Kameni islet. Our post-unrest deflating sill at the northern intra-caldera area is in good agreement with the unrest Mogi source of Reference [9], for their better constrained solution (period P1, Jan-Jun 2011), in terms of number of GPS stations and larger deformation gradients. On the whole, the inefficacy to well constrain magma sources using single SAR LOS observations, or GPS measurements that suffer from poor spatial sampling should be underlined. The utilization of multiple geodetic assets is thus beneficial and should be certainly taken into account when designing volcanic monitoring systems.
The InSAR results and modelling for the Nea Kameni deformation source of this study demonstrate a comparable pattern to previous ERS and ENVISAT observations covering the 1992-2010 pre-unrest period [3,4,56] (Figure 12a). The source's identical location and volume variation rate between the two intervals points to the hypothesis of a steady deformation signal. By assuming the pre-unrest deformation rate (ERS/ENVISAT) as background motion, and subtracting it from the Sentinel-1 post-unrest deformation rates, we were able to reveal the post-unrest response signal ( Figure 12). Thereafter, to simplify the comparison between different LOS incident angles, we only considered the vertical component of the displacement. The outcome specifies a ramp tilting towards the north with rate of 5.5 ± 2.7 mm/yr. Modelling results adopting only the sill (the caldera-wide source) show an identical tilt at Nea Kameni (3.4 mm/yr), due to the deflation of the sill source at the northern caldera ( Figure 12d). The computed tilt rate is within the uncertainty of SAR estimates. This demonstrates that local subsidence at Nea Kameni can be considered of the same amount during 1992-2010 and 2014-17, while the residual tilt of 2014-17 at Nea Kameni is attributed to the caldera-wide deflating source of the post-unrest period.
The existence of the local post-unrest shallow deformation source at Nea Kameni, having the same behavior as in the pre-unrest period in terms of both deformation rate and spatial pattern, suggests its continuous action during the unrest. Modelling misfits over Nea Kameni for the unrest period indicate a systematic over-estimation of both SAR and GPS data, as demonstrated in other studies only considering a single offshore magma source [4][5][6]9]. Therefore, the same volcanic activity is supposed to occur at Nea Kameni during the unrest, consistent with the pre-/post-unrest steady subsidence of -5 mm/yr [3,4]. However, such deformation cannot be isolated in the uplift data, since a total amount of 9 cm was measured at Nea Kameni during the unrest. In order to quantify the effect of the local source to the main source of the unrest, we ran a model with two sources. The first one was the unrest source as in Reference [4], and the other one was the Nea Kameni source as constrained in this work. The percentual change of the volume variation rate in the case of a single active source, as opposed to the case where both sources are active, amounts to 0.2%. So, in case the deflation at Nea Kameni is taken into account during the unrest, the inflating magma source is still located at the same position and depth (3.7 km), undergoing an increase of only 0.2% of the volume variation rate. This is not surprising, since the unrest generated large displacements that were actually masking out the deformation at Nea Kameni. In our opinion this prevents an independent constraint of the Nea Kameni source from data, since its action is partially concealed by the unrest. The unrest source has 7 3 4 3 The outcome specifies a ramp tilting towards the north with rate of 5.5 ± 2.7 mm/yr. Modelling results adopting only the sill (the caldera-wide source) show an identical tilt at Nea Kameni (3.4 mm/yr), due to the deflation of the sill source at the northern caldera (Figure 12d). The computed tilt rate is within the uncertainty of SAR estimates. This demonstrates that local subsidence at Nea Kameni can be considered of the same amount during 1992-2010 and 2014-17, while the residual tilt of 2014-17 at Nea Kameni is attributed to the caldera-wide deflating source of the post-unrest period.
The existence of the local post-unrest shallow deformation source at Nea Kameni, having the same behavior as in the pre-unrest period in terms of both deformation rate and spatial pattern, suggests its continuous action during the unrest. Modelling misfits over Nea Kameni for the unrest period indicate a systematic over-estimation of both SAR and GPS data, as demonstrated in other studies only considering a single offshore magma source [4][5][6]9]. Therefore, the same volcanic activity is supposed to occur at Nea Kameni during the unrest, consistent with the pre-/post-unrest steady subsidence of −5 mm/yr [3,4]. However, such deformation cannot be isolated in the uplift data, since a total amount of 9 cm was measured at Nea Kameni during the unrest. In order to quantify the effect of the local source to the main source of the unrest, we ran a model with two sources. The first one was the unrest source as in Reference [4], and the other one was the Nea Kameni source as constrained in this work. The percentual change of the volume variation rate in the case of a single active source, as opposed to the case where both sources are active, amounts to 0.2%. So, in case the deflation at Nea Kameni is taken into account during the unrest, the inflating magma source is still located at the same position and depth (3.7 km), undergoing an increase of only 0.2% of the volume variation rate. This is not surprising, since the unrest generated large displacements that were actually masking out the deformation at Nea Kameni. In our opinion this prevents an independent constraint of the Nea Kameni source from data, since its action is partially concealed by the unrest. The unrest source has a volume variation of the order of 10 7 m 3 /yr, while the Nea Kameni source amounts to 10 4 m 3 /yr, so it fits the effect of the Nea Kameni source on the unrest source being of the order of 10 −3 .
The Nea Kameni deformation source is suggested by Reference [4] as the effect of variations within the shallow hydrothermal system, the existence of which is reported by Reference [68] at 800-1000 m depth. A different mechanism is suggested by Reference [57], whereas slow subsidence could reflect the thermal cooling and the load-induced relaxation of the substrate due to lava flows emitted between 1866 and 1870. In view of our SAR results, and the identification of the same slow subsidence signal before and after the unrest, this latter scenario could not also be ruled out. Based on geodetic data [9] and petrological studies [69,70], deeper volcanic-related processes beneath Nea Kameni are also proposed. If, over and above, we take into account that historic lava flows built up the intra-caldera islets of Palea and Nea Kameni and that the last lava emplacement took place in 1950, the involvement of still-ongoing deep processes beneath Kameni islets could support this hypothesis.
As a final remark, our results show that a change (decrease) in the subsidence rate at Nea Kameni could reflect the onset of a renewed unrest, and as a consequence, the geodetic monitoring of ground deformation represents a key to deciphering the dynamics of the Santorini shallower plumbing system, even for periods not accompanied by seismic activity or increased degassing. The Kameni deformations are therefore of primary importance to correctly interpret future unrest episodes and evaluate the volcanic hazard.

Conclusions
The MT-InSAR analysis of the dense temporal Sentine-1 data series allowed us to measure the post-unrest ground deformation of the Santorini volcano. The results provide additional insights on the recognition of an annual, seasonal periodicity in the displacement time series, affecting the accurate estimation of the deformation gradients. Whether this oscillation pre-existed during the unrest period, or even before, cannot be ascertained due to the non-systematic temporal coverage of former SAR acquisitions. Hence, this new estimation provides an important step for a more reliable assessment of the volcano deformation.
The new volcano state after the unrest period is confirmed by the geodetic analysis of multiple SAR sensor data (Sentinel-1, Radarsat-2 and TerraSAR-X) for the period 2012-17. The post-unrest response to the 2011-12 inflation episode is well explained by a shallow sill-like source at 2 km depth. This source is located just above the~4 km deep inflation source responsible for the 2011-12 uplift. The wide observed deflation extending up to Thera Is. (Cape Skaros and northern caldera walls) provides evidence that the volcano apparently remains still under a recovery state. According to the geochemical and isotopic data on gas emissions, this deflation reflects the passive degassing of the shallow (top) portion of an intrusion of relatively poorly evolved magma emplaced during the unrest. The degassing pathways could be represented by the Kameni fault and fractures, which were re-activated by earthquakes during the unrest period, allowing an increase in local permeability.
The presence of a steady subsidence source at Nea Kameni, in accordance with the pre-unrest period, led to the re-evaluation of the 2011-12 unrest. Our interpretation model suggests the co-existence of the Kameni source during the unrest, having however a lower impact compared to the larger deformations induced by the inflation source.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/3/259/s1, Text 1: Supplementary material of the manuscript. Figure S1: Temporal separation versus normal baseline plots displaying the InSAR pairs (connecting lines) considered in the MT-InSAR processing. Sentinel-1 scenes in blue color were not considered in the solutions. Figure S2: Sentinel-1 MT-InSAR LOS displacement rate uncertainties for (a) Ascending orbit 029; (b) Descending orbit 109; (c) Descending orbit 036. Uncertainties of estimated rates (see Figure 3) are given relative to the reference point, marked by a rectangle. The line-of-sight and azimuth directions of the satellite are displayed by blue and black arrows, respectively. Values in degrees correspond to the incidence angles. Figure S3: Comparison between subsampled data (first column) and modeled data (second column) for the single source model of the post-unrest. The third column reports the residuals (observed data minus modeled data). (a-c) Ascending orbit track 029 (a029); (d-f) Descending orbit track 036 (d036); (g-i) Descending orbit track 109 (d109); (j-l) Descending radarsat-2 data (drs2). The shaded circle is the surface projection of the sill-like source with its actual radius (about 2000 m). The chi-square for this model is 3.5. Figure S4: Comparison between subsampled data (first column) and modeled data (second column). The third column reports the residuals (observed data minus modeled data). (a-c) Ascending orbit track 029 (a029); (d-f) Descending orbit track 036 (d036); (g-i) Descending orbit track 109 (d109); (j-l) Descending radarsat-2 data (drs2). The star is the center of the Kameni source, the shaded circle is the surface projection of the sill-like source with its actual radius. Figure S5: Posterior Probability Density (PPD) functions retrieved by the non-linear joint inversion of the two sources. (a,b) Coordinates of the centre of the sill-like source (S), (c) its depth, (d) radius, (e) potency, i.e., overpressure ∆P vs rigidity µ; (f,g) Coordinates of the spherical source (Mogi source, M), (h) its depth, (i) its volume variation. Easting and Northing are in UTM-WGS84 projection, zone 35. The dashed line is the mean model. Figure S6: Two-dimensional PPD distributions. Each panel represents the 2D distribution of the inverted parameters. Contour lines every 10% confidence. The red cross is the mean model listed in Table 3. Each panel reports the numbering of two parameters for which the 2D-distribution is computed, considering the exact order as reported in Figure S5 and Table 3