Two-Dimensional InSAR Monitoring of the Co- and Post-Seismic Ground Deformation of the 2021 Mw 5.9 Arkalochori (Greece) Earthquake and Its Impact on the Deformations of the Heraklion City Wall Relic

: Contributing to the United Nations 2030 Sustainable Development Goals (SDGs) within Target 11.4 “Strengthen efforts to protect and safeguard the world’s cultural and natural heritage”, it is critical to monitor the spatial and temporal stabilities of cultural heritages. The study of the interactive relationship between earthquakes and the protection of cultural heritages needs to be strengthened. On 27 September 2021, the destructive Mw 5.9 Arkalochori earthquake occurred ~25 km away from the city of Heraklion (Greece) where the Heraklion City Wall (HCW), a representative cultural heritage of Greece and Europe, was located. This offered a proper case to investigate the shortcomings aforementioned. Here, we intend to set up and answer the following three questions (Whether, Where and What, 3Ws): Whether there were impacts on the HCW caused by the Arkalochori earthquake? Where did the maximum deformation occur? What was the relationship between seismic deformation between the epicenter and the HCW over time? We performed two-dimensional (2D) InSAR measurements for both co-seismic and post-seismic deformations using the ascending and descending Sentinel-1A SAR images. The spatial-temporal characteristics of Up–Down (UD) and East–West (EW) were revealed. The 2D co-seismic deformation ﬁeld showed that the near-ﬁled deformations were dominating compared with the deformations at the HCW, the UD deformation was mainly featured with subsidence with a maximum value of ~21 cm, the EW deformation was ~9 cm westward and ~10 cm eastward. The time-series measurements showed that: (1) temporally, the HCW responded quickly to the Arkalochori earthquake, and the accumulative deformations at the seven different bastions of the HCW showed the same trend as the near-ﬁeld area over time. (2) Spatially, the closer to the Mw 5.9 epicenter, the larger the deformations that occurred. (3) The EW and UD deformation trends of the HCW that were consistent with the Mw 5.9 epicenter were interrupted at the middle time spot (22 January 2022), indicating the inﬂuence of another earthquake sequence consisting of eight earthquakes with magnitudes larger than 3.5 that happened on 16–18 January 2022. Respectively, to summarize and address the aforementioned 3Ws based on the post-seismic analysis accomplished by the MSBAS method, the Arkalochori earthquake did affect the HCW; besides, the inﬂuences of the ~13 km earthquake sequence were also detected; the nearest part to the epicenter suffered the most; the deformation trends of the HCW were approximately the same with the epicenter area of the Arkalochori earthquake both in the UD and EW directions. the epicenter. The results showed that the co-seismic deformation is F.C., W.Z., H.L. J.L.; M.Z.; F.C., H.L. and I.P.


Introduction
Cultural heritages are irreplaceable legacies that witness and record the origin together with the development of human civilization [1]. Threatened by natural disasters, particularly earthquakes, heritage protection is facing severe challenges [2,3]. Due to the degradations through the years and the distortions introduced by earthquakes, cultural heritages are prone to lose tremendous historical wealth, cultural reputation, and artistic significance [4]. Accordingly, it is urgent to strengthen the research on the relationship between the deformations of heritage sites and earthquakes [5,6].
On 27 September 2021, the Mw 5.9 Arkalochori earthquake struck Crete Island causing injuries and facility destruction. According to the relocated reports of [7], the epicenter was located at 35.1413 • N 25.2729 • E, which was~25 km away from the city of Heraklion ( Figure 1). In the downtown area of the city of Heraklion, there was an immortal cultural heritage named the Heraklion City Wall (HCW), which was built from 1462 to the 16th century and was considered to be one of the best-preserved Venetian fortifications in Europe. The main purpose of this study was to investigate the interactions and potential relationships between the earthquake deformations and the induced deformations of cultural heritages, and the Mw 5.9 Arkalochori earthquake and the HCW case offered a proper insight into this problem. InSAR technology has been successfully implemented in heritage health monitoring [4,[8][9][10][11][12][13] and infrastructure health monitoring [14]. Specifically, for detecting the displacements of heritage sites related to earthquake events, refs. [15][16][17] have illus-trated the capability of InSAR with the support of Sentinel-1 SAR images in monitoring archaeological sites and monuments after earthquakes. For the Arkalochori earthquake, in [7,[18][19][20], 2D co-seismic deformations of the near-field caused by the Mw 5.9 Arkalochori earthquake were observed using D-InSAR method, but none of them took the atmospheric phase screen (APS) in the interferograms into consideration. The study by [7] also detected a post-seismic deformation that occurred 2-5 days after the Mw 5.9 event, denoting that the longer the time after the earthquake, the more obvious the post-seismic ground motion was. However, the study by [21] showed that the post-seismic deformation of the following 30 days after the Mw 5.9 event was hard to detect using D-InSAR interferometric pairs. However, studies of the earthquake-induced deformations on the HCW were still lacking. Moreover, the longtime ground deformation extraction of the Arkalochori earthquake by the InSAR time-series analysis had not been reported.
Compared to the D-InSAR technique, the time-series InSAR method is more capable to detect deformations with a small magnitude over a long period [22]. In the meantime, the APS in the interferograms could weaken or strengthen the signal of ground motion, making it even harder to obtain the exact ground behaviors. In addition, the importance of longer-term evolution analysis comes from the recordings of deformation history over time, which is essential for many analyses, e.g., studying the kinematics of a given phenomenon (quiescence, activation, acceleration, etc.), correlation with driving factors, etc. [22]. The existing post-seismic InSAR studies were mainly focused on several days after the Mw 5.9 mainshock, we extend the period after the mainshock to~8 months, not only to monitor the spatial-temporal evolution of the post-seismic deformation but also to investigate if there were potential reactions between the time-series accumulative deformation of the near-field and the HCW which was~25 km away from the epicenter.
Here, we used D-InSAR and multidimensional InSAR time-series methods to measure not only the ground deformation caused by the 2021 Mw 5.9 Arkalochori earthquake but also the deformation of the HCW, aiming to investigate if there were interactive relationships between them. To start with, we decided to find out whether there were effects that earthquakes could introduce to the nearby heritage sites. Then, for large heritage sites like the HCW, where or which part suffered most from the earthquake? Eventually, we intended to reveal how the Arkalochori earthquake influenced the 2D deformations of the HCW and what were the relationships between them.
This study first performed D-InSAR processing with the Generic Atmospheric Correction Online Service for InSAR (GACOS) corrections [23][24][25], implemented to mitigate the tropospheric phase delay in both descending and ascending interferometric pairs. Then, a decomposition of Line-of-Sight (LOS) measurement from both tracks was calculated to obtain the UD and EW deformation in the near-field region. In addition, the Multidimensional Small Baseline Subset (MSBAS) method [26][27][28][29][30][31][32][33][34] was utilized to get the 8 months 2D time-series measurements using 40 Sentinel-1A SAR images from 30 September 2021 to 22 May 2022.

Study Area
Crete Island, the largest island with the most population in Greece, is generally mountainous (Figure 1). It is located at a tectonic active zone where the African Plate subducts beneath the Aegean Sea [21]. Earthquakes occur frequently in this region due to the subduction process along the Hellenic subduction zone [18]. Crete has formed a unique symbol of the cultural heritage in Greece.
After the Mw 5.9 Arkalochori earthquake, an earthquake sequence larger than M 3.5 occurred on 16-18 January 2022, which was~13 km away from the HCW. The sequence consisted of 8 earthquakes with a medium-low magnitude, thus, it was hard to detect in conventional D-InSAR interferograms. The burial depth ranged from 21.2 km to 25.5 km as reported by the local official catalog from the Institute of Geodynamics of the National Observatory of Athens (NOA-IG). The effect of the earthquake sequence was observed in the time-series InSAR observations.
The HCW was also known as the Fortifications of Heraklion. It consists of seven triangle-shaping bastions ( Figure 2). The length of the HCW is about 4.5 km. In the following parts of the study, we mainly focused on the 2D deformation rates and the time-series evolutions over time at the seven bastions.

Materials and Methods
According to the different features of ground deformation, the applicable methods are different. The D-InSAR technique is currently mainly used in co-seismic deformation monitoring, which can extract the deformation field information in a quick, efficient, and economical way. While the InSAR time-series analysis is capable of high-precision measurements of long-term deformation with a relatively small magnitude by analyzing a large number of time-series SAR images to separate different phase components according to the temporal and spatial variation characteristics. In this study, the GACOS-corrected D-InSAR and multi-dimensional InSAR time-series methods were jointly used for 2D deformation measurements of the Mw 5.9 Arkalochori earthquake and their relations with the stabilities of the HCW.

SAR Datasets
The C-band Sentinel-1A SAR images in Terrain Observation with Progressive Scans (TOPS) mode were used in this study, belonging to 3 different tracks. Two tracks were the ascending tracks (Track 102 and Track 029) and another was the descending track (Track 109) (Figure 1). The ascending track was flying in a roughly south-to-north direction while the descending track was flying in a roughly north-to-south direction. Note that the final 2D co-seismic and time-series results were limited to the common spatial coverage of the SAR images used for the decomposition.
Because almost all SAR data are acquired using a right-looking geometry [35], to extract the 2D deformations, at least 2 tracks from one ascending and one descending orbit were needed. For D-InSAR processing, the ascending Track 102 and descending Track 109 were processed. The detailed parameters of these SAR images were shown in Table 1. The 2 frames of the ascending Track 102 were stitched together to get better coverage of the area of interest in this study. The revisit time interval of the Sentinel-1A is 12 days. All the SAR images used in our study were right-looking as marked in Figure 1.
For MSBAS time-series processing, 20 SAR images from the ascending Track 029 and 20 SAR images from the descending Track 109 were processed. The related detail parameters of these images were shown in Table 2.

D-InSAR Processing
The D-InSAR processing was accomplished utilizing the open-source InSAR Scientific Computing Environment version 2 (ISCE2) software [36]. The topographic phase was removed using the 1 arcsecond (~30 m) Shuttle Radar Topography Mission (SRTM). The multi-look factor in range and azimuth was set to be 8:2, respectively. The interferometric phase was filtered using Goldstein filtering [37]. Then, the filtered phase was unwrapped by the integrated correlation and unwrapping (ICU) method, which was an enhanced version of the traditional residue-based unwrapper established in [38]. The ICU unwrapper was a built-in method in the ISCE2 package. By integrating the analysis of phase gradients, amplitude, and coherence information, the ICU unwrapper was guided to determine multiple starting seed pixels [39,40]. We use external GACOS [25] data to conduct the APS correction to the unwrapped phase. Finally, the APS-corrected phase was converted to deformation (unit in cm) and geocoded to the WGS-84 coordinates.

MSBAS Processing
The MSBAS method was developed by [31]. Jointly making use of the InSAR observations from multiple SAR sensors, looking angles, and flight directions, the MSBAS method offers the opportunity to reveal the UD and EW deformation of the ground with enhanced time sampling. Here, we briefly present the key equations of the MSBAS algorithm in (1), (2), and (3). More detailed theories of the MSBAS algorithm could be referred to in [31,33,34].
The matrix form of the MSBAS algorithm can be written as [31]: where∅ consists of InSAR-observed LOS deformations from multiple sources, L is the Tikhonov regularization matrix [41], λ is the regulation parameter [31,33,34], V EW and V UD are the unknown EW deformation and UD deformation that need to be determined through the equation. MatrixÂ is defined based on matrix A, a matrix constructed from the time intervals between consecutive dates of SAR acquisition [31,34]. MatrixÂ can be written as: where {µ EW , µ UD } is formed according to the decomposition using the looking geometries of the LOS measurements for different SAR orbits. Equation (2) can be further written as: where θ is the incidence angle, α is the heading angle of the SAR satellite.
In this study, 157 interferograms derived from 20 scenes of ascending SAR images together with 152 interferograms derived from 20 scenes of descending SAR images were firstly unwrapped using the Statistical-cost, Network-flow Algorithm for Phase Unwrapping (SNAPHU) method [42]. Then, the unwrapped phases were corrected using GAGOS. Furthermore, they were resampled to the same grid and cut to the common area in the WGS-84 geosystem. In the following MSBAS processing, we applied the Tikhonov regularization [41] to solve the ill-posed problem caused by rank deficiency. The regularization was similar to a temporal filter that removed the high-frequency noises [31]; thus, no additional temporal filtering was applied. Then, a 2D low-pass filtering in the spatial domain was implemented to mitigate the residue noise. A trade-off value of 0.25 for the regularization parameter [31,34] was used in this study. The outputs of the MSBAS processing were the 2D deformation rate maps and the accumulative time-series behaviors of certain points. All interferograms were produced by the batch processor of the open-source GMTSAR software [43] based on the standard SBAS network of the interferometric combinations.

Reference Point Selecting
InSAR measurements, including D-InSAR and time-series methods, are relative to a reference point that is stable during the SAR acquisitions [44]. To maintain the homogeneity of the multiple InSAR measurements, in this study, we selected the same reference point to calibrate all the interferograms of D-InSAR and stack calculated from ascending and descending tracks, respectively. Coherence is a common indicator for reference point selection [44]. The coherence threshold of the reference point was commonly suggested to be above 0.85 [45]. The distance between the reference point and the area of interest should not be too large because of potential APS propagations [46].
We obtained 4 coherence maps in total (Figure 3), two of them were from ascending and descending D-InSAR processing. The other two coherence maps are calculated by averaging the 157 ascending interferograms and 152 descending interferograms. Finally, the point that was close to the HCW was selected which was maintained to be highly coherent in all of the 4 coherence maps. The coherence of this reference point was 0.927, 0.939, 0.928, and 0.940, respectively, in ascending D-InSAR (Figure 3a), descending D-InSAR (Figure 3b), and ascending stack (Figure 3c), descending stack ( Figure 3d). The distance between the reference point and the epicenter was~25 km.
The selected reference point was next to the HCW. The SAR coherence had indicated high values in 4 independent processing routines. We also check the location of the reference point in an optical Google Map image. The reference point was located at a low building area as shown in Figure 4.

GACOS Corrections
We applied GACOS corrections to mitigate the APS in D-InSAR interferograms. As shown in Figure 5, the qualities of the interferograms were improved. There were obvious tropospheric phase delays simulated by GACOS (Figure 5b,e). The deformations contributed by the APS reached~6 cm to~8 cm, which was approximately the same magnitude as the co-seismic deformation, which was not ignorable. The simulated APSs from GACOS had the same pattern as the phases in the uncorrected interferograms. After the GACOS corrections, the corresponding phase patterns were mitigated (Figure 5c,f). Note that we also conducted a phase detrending to remove the long wavelength phase delays according to the GACOS manual (www.gacos.net/ (accessed on 16 October 2022)).

The UD and EW Deformations of the Arkalochori Earthquake
One-dimensional LOS deformation was insufficient to reveal the real ground motions [47,48]. We used the GACOS-corrected interferograms for the decomposition calculation to obtain the UD and EW deformations. The deformation of the North-South direction was ignored because current SAR satellites operated in near-polar orbits, the LOS observations measured from both ascending-and descending-looking geometries were strongly sensitive to the deformation along the UD direction, modestly sensitive along the EW direction, and weakly sensitive along the North-South direction [47]. The co-seismic deformation shown in Figure 6 suggested the feature of the ground motions caused by the earthquake was mainly downward (or subsidence). The maximum downward deformation reached~21 cm. The eastward deformation was about 10 cm and the westward deformation was~9 cm. Compared with the near-field deformations, the deformations of the~25 km where the HCW was located were much smaller (Table 3).  The MSBAS processing is an advanced SBAS-based time-series algorithm that required stacks of interferograms from at least two different view geometries. Based on the conventional small baseline pair selection principle in the SBAS method [49], we constructed two stacks of interferograms from the ascending and descending track, respectively. For the ascending stack, the image acquired on 22 May 2022 was selected as the super reference image, as for the descending stack, the image acquired on 22 April 2022 was selected as the super reference. The super reference images were selected since they were located in the middle of the clouds of baselines (Figure 7). All other images were co-registered with the selected super reference images, respectively, for the two stacks. The threshold of the spatial and temporal baseline is 150 days and 150 m. The interferometric combinations were shown in Figure 7. We get 157 and 152 pairs of interferograms, respectively from the ascending and descending tracks. The GACOS plus the spatial-temporal filtering strategy was more capable to get better time-series results [50]. Accordingly, we also applied the GACOS corrections to all the total 309 interferograms to mitigate the APS.

MSBAS Results
The outputs of the MSBAS processing are 2D deformation rate maps and accumulative deformation evolution of a certain point over time. As introduced in detail in [33,51], the MSBAS algorithm made a boundary correction along the timeline of SAR acquisition dates, that was, at both ends of the period, the first and the last interferogram were compensated, making the first and last SAR acquisition date to be assumed the same with the adjacent one. Particularly in our case, there were 40 SAR images in total, while the final outputs of the MSBAS processing consisted of 38 independent time spots throughout the entire period.
Starting with the GACOS corrected interferograms, the MSBAS processing was conducted. By joint decompaction, the UD and EW deformation velocity maps were shown in Figure 8. Note that, the MSBAS method could only extract the common footprint of all the data sets that were used in the processing (Figure 3c,d) [31]. For the epicenter area, the EW deformation rate was less than 1 cm/year slightly toward the west, while the UD deformation rate was~2 cm/year downward, indicating the effect of the mainshock continued with a long-lasting post-seismic ground motion in the following 8 months mainly in the vertical direction. For the HCW, the ground was relatively stable compared with the epicenter area. Both the EW and UD deformation rates were less than 1 cm/year. The database of the faults plotted in Figure 8 was referred to in [7]. Then, we focused on the deformation velocity of the HCW relic. The zoom-in 2D velocity maps were shown in Figure 9. In the EW direction, the westward deformation velocity was smaller than the eastward one. The western part moves faster than the eastern part. The EW deformation velocity of bastions NO. 1, 2, 3, and 4 ranges from~0.25 cm/year to~1 cm/year toward the east, respectively. In the vertical UD direction, the HCW generally showed motion uplift, the southern parts moved faster with a velocity of~0.8 cm/year, while the northern parts moved slower with a velocity of~0.2 cm/year to~0.4 cm/year. The deformation velocity of the Martinengo Bastion (NO. 4) was the fastest compared with the other six bastions, this was because the Martinengo Bastion was mostly close to the Mw 5.9 epicenter. Furthermore, the 2D accumulative defamations were investigated. Note that the near-field time-series deformations were extracted at a point close to the Mw 5.9 epicenter (shown in Figure 8), because at the exact location of the epicenter, the values were not existed due to decorrelations along the more than 300 interferograms. Figure 10 showed the accumulative EW deformations of the seven bastions of the HCW together with the epicenter. It was observed that before the marked time spot, the deformation trend of the HCW was the same as the epicenter, showing a deformation behavior slightly toward the west and then moving toward the east till the marked time spot. The consistent deformation trends between the epicenter and the HCW were interrupted at time spot NO. 19. It was also shown in Figure 8 that there was an earthquake sequence consisting of eight earthquakes with magnitudes larger than 3.5, according to the Institute of Geodynamics of the National Observatory of Athens (NOA-IG), that happened on 16-18 January 2022. The locations of the eight events were~13 km to the east of the HCW according to the records by NOA-IG. Starting from time spot 19, the EW deformation trend of HCW was different from the Mw 5.9 epicenter. Specifically, the EW deformation of the Mw 5.9 epicenter was eastward, with a relatively fast rate (slope of the deformation curve), but for the HCW, five bastions showed obvious westward motions and the other two bastions showed~24 days of stabilities. It was observed in Figure 10 that the bastions NO. 3, 4, 5, 6, and 7 were moving toward the west during the NO. 19 and 23 intervals (24 days) with a rate of~1.56 cm/year while the epicenter was moving toward the east during the same period with a rate of~1.78 cm/year. The location of the earthquake sequence which happened 4-6 days before the time spot numbered 19 was~13 km to the east of the HCW, as a result, the bastions NO. 3, 4, 5, 6, and 7 were influenced more because they were nearer to the epicenters of the earthquake sequence. The previous same trends as the Mw 5.9 earthquake were interrupted at the 19th time spot.  Figure 11 showed the UD accumulative deformations of the seven bastions of the HCW together with the epicenter. Before the marked time spot, the UD motions of the HCW also showed the same trend as the epicenter. This consistent trend was interrupted after the marked time spot, too. After the NO. 19 time spot, the deformation of the Mw 5.9 epicenter area continued to show subsidence, while all of the seven bastions of the HCW were moving uplift during the same time, indicating that the behaviors of the HCW deformations were no more correlated with the Mw 5.9 event. During the whole 8 months observation period after the Mw 5.9 earthquake, the maximum subsidence of the Mw 5.9 epicenter area was~1.06 cm. Maximum uplift of the HCW was detected at the NO.4 Martinengo Bastion, reaching~0.35 cm.

Discussion
Previous research results [52][53][54][55][56][57][58][59][60] have proven that InSAR technology is mature and provides reliable information on co-seismic and post-seismic deformation monitoring. Although we didn't have the Global Navigation Satellite System (GNSS) data for direct comparison and validation, as discussed in [21], InSAR had shown an accuracy comparable to the GNSS measurements for the Mw 5.9 Arkalochori earthquake.
The relationship between the Mw 5.9 earthquake and its induced deformation for the HCW was clarified by the InSAR time-series analysis. Meanwhile, we noticed a sudden change in the time-series deformations (Figures 10 and 11, Spot No.19). What caused that sudden change? Was the correlated same trend interrupted by some other disturbances, or there were no relations? The acquisition dates and number identifications of the 40 SAR images together with the earthquake events were shown in Figure 12.
We checked the earthquake catalog released by NOA-IG (https://bbnet.gein.noa.gr/ HL/databases/database (accessed on 16 October 2022)) and found there was an earthquake sequence that happened~13 km away from the HCW on 16-18 January 2022. This earthquake sequence consisted of eight earthquakes, each of which was larger than M 3.5. As shown in Figure 12a, the eight earthquakes happened between the 18th and the 19th time spot along the timeline. The trend interruptions in both EW and UD accumulative deformation curves were observed in the 19th to 20th interval. Since the earthquake sequence marked in Figure 8 happened 1-3 days after the 18th spot (Figure 12a), and the corresponding deformation trend inflections were not detected before the 19th spot. It could be assumed that there was at least a 4-6 day delay for the HCW to show reactions due to the seismic influence of the earthquake sequence, and these reactions were captured by InSAR within a single interval (6 days) in our observation. According to [61], earthquakes could cause potential deformation to a dam; hence, it was assumed that thẽ 13 km earthquake sequence (Figure 12b) might be the reason that the previous consistent trend between the HCW and the Mw 5.9 earthquake was interrupted. The magnitude of the eight earthquakes in the earthquake sequence was 3.5-3.8 (Figure 12c), and the depth of the epicenter was reported to be 21.2-25.5 km, making it difficult to be detected through D-InSAR measurements. However, it was detectable in the InSAR time-series monitoring thanks to the advantage of the capability and sensitivity in long-lasting micro deformation measuring. The changing signs of the curve trends in both EW and UD components between spots 10 and 13 were assumed to be related to the aftershocks and afterslips. The dominant driven force of the deformations was interpreted to be the earthquake. Once the main shock occurred, the strain and stress were changed, and the elastic crust started to adjust to the new situation. The afterslips and aftershocks continued to influence this progress. From Figures 1 and 8, we could see that in the period between spots 10 and 13, there were continuous aftershocks in the near-field area but there were no aftershocks larger than M 3.5 near the HCW (Figure 12) until the 18th time spot. After the Mw 5.9 mainshock, the dominant driven force started to take over, and the deformations of HCW acted corresponding to it. The trend sign of EW deformation of HCW at the 10th spot, but the near-field area of the Mw 5.9 epicenter started this sign change~2 weeks later, this would be related to the continuous aftershocks around the Mw 5.9 epicenter area, but these aftershocks were relatively weak far away (20 km or more) from the HCW; thus, they had less effects on the deformation of HCW. For the UD deformations, the deformations of HCW did not show clear signs of changes but the trends tended to be smoother. The UD sign change of the "Near-field" area was also assumed to be related to the continuous aftershocks in the 10-13 observation intervals because according to the NOA-IG catalog, there were two aftershocks larger than M 3.5. One happened on 3 December 2021 with a magnitude of 3.7 and was buried at 8.3 km in depth, and the other happened on 7 December 2021 with a magnitude of 3.7 and was buried at 9.4 km in depth. The time of these two aftershocks located between the 12th and the 13th spot in the timeline ( Figure 12) and the signs of changes of the "Near-field" curves in both EW ( Figure 10) and UD ( Figure 11) directions also occurred in the same period.
There were small irregular disturbances in both the UD and EW deformation evolutions over time at the HCW relic. These could be the effects related to human activities because the HCW was located in the downtown area of Heraklion city, some parts of the HCW were even reconstructed into modern buildings, playgrounds, or gardens viewed from Google Maps (Figure 4). The spatial-temporal evolutions of deformations were the combinations and reflectors to multiple driven forces [61,62], including earthquakes nearby and far away, the levels of surface water or groundwater, human activities, temperature changes, tectonic movements, local long-term subsidence due to the vicinity to the coastline, and so on. In this case, earthquakes played dominating roles that controlled the general trends of the HCW deformations.
Moreover, there were still residue phase noises in the GACOS corrected and the spatial-temporal filtered interferograms. Fully eliminating the multiple sources of noise remains to be a complex and challenging task [63,64]. Current methodologies could reduce or mitigate the effects of APS and improve the quality of the interferometric measurements among which the GACOS method has been proven to be successful and effective [63,65]. Moreover, the GACOS corrected residue phase components were more suitable to meet the conventional assumption of the APS behaviors in time and space domains [50], making the spatial-temporal filtering option more reliable. These all demonstrated that we had done our best to deal with the APS in this case study under the current condition of the state-of-the-art technology level.
The driven forces that cause the deformations of the HCW could be complicated, among which, we assume that the shake related to the earthquakes played an important role. It can be inferred from Figures 10 and 11 that, the maximum accumulative deformations were no more than 0.2 cm and 0.35 cm in the EW direction and UD direction, respectively. The potential risk of collapse or crack induced by these deformations should be taken into serious consideration in heritage conservation. By investigating and obtaining the amount and direction of the earthquake-induced deformation, the corresponding protection actions should be undertaken. Although the Mw 5.9 and the earthquake sequence of M 3.7, which was close to the HCW, did not cause significant damage to it, the earthquakes did introduce deformations to the wall. The potential risk should be taken care of before there would be further destruction. We also expect further multiple methods such as GNSS, field surveys, modeling, etc., to jointly work on the monitoring and protection of the HCW.

Conclusions
Aiming at promoting the heritage impact assessment, the purpose of this study was to investigate the relationships between the Mw 5.9 earthquake and the HCW that was 25 km away from the epicenter. The results showed that the co-seismic deformation is much larger than the post-seismic deformation. The InSAR time-series method revealed longer-term evolution since it recorded the deformation history over the observed period, which showed the correlation of deformation with the seismic-driven factors. The nearfield co-seismic defamation feature was mainly subsidence, the horizontal EW co-seismic deformation was half of the magnitude as the vertical deformation. The APS in the uncorrected interferogram was effectively mitigated by the GACOS correction. The MSBAS time-series analysis showed approximately the same deformation trends between the HCW and the epicenter area of the Mw 5.9 Arkalochori earthquake, indicating the Arkalochori earthquake did influence the deformations of the HCW. The sudden trend interruptions of both the UD and EW deformations, that occurred on 22 January 2022, indicated the existence of an earthquake sequence that happened even nearer on 16-18 January 2022, showing the deformations of the HCW were potentially influenced by another earthquake. The conclusions were made based on the post-seismic analysis using the MSBAS method. To directly answer the aforementioned 3 questions (3Ws): (1) Whether there were any influences between the Arkalochori earthquake and the HCW? Yes, not only the Arkalochori earthquake but also another earthquake sequence nearby contributed to the deformations of the HCW. (2) Where of the HCW suffered most? The closer to the epicenter, the larger the deformation occurred. (3) What kind of relationship there was? The HCW showed approximately the same deformation trends in both UD and EW directions as the epicenter area of the Arkalochori earthquake.
Space technologies link earthquake research and heritage protection together. The study provides new insight for investigating the interactive relations between earthquakes and cultural heritages and could be extended to more case studies such as ancient buildings, tombs, and inscriptions on precipices made from different materials.  Data Availability Statement: Sentinel-1 SAR images can be quarried and downloaded from European Space Agency (ESA) and the Alaska Satellite Facility (ASF) free of charge. The Precise Orbit Ephemerides of the data are freely available at (https://scihub.copernicus.eu/gnss/#/home (accessed on 16 October 2022)). The GACOS data are provided from (http://www.gacos.net/ (accessed on 16 October 2022)) freely.