Analysis of Secular Ground Motions in Istanbul from a Long-Term InSAR Time-Series ( 1992 – 2017 )

The identification and measurement of ground deformations in urban areas is of great importance for determining the vulnerable parts of the cities that are prone to geohazards, which is a crucial element of both sustainable urban planning and hazard mitigation. Interferometric synthetic aperture radar (InSAR) time series analysis is a very powerful tool for the operational mapping of ground deformation related to urban subsidence and landslide phenomena. With an analysis spanning almost 25 years of satellite radar observations, we compute an InSAR time series of data from multiple satellites (European Remote Sensing satellites ERS-1 and ERS-2, Envisat, Sentinel-1A, and its twin sensor Sentinel-1B) in order to investigate the spatial extent and rate of ground deformation in the megacity of Istanbul. By combining the various multi-track InSAR datasets (291 images in total) and analysing persistent scatterers (PS-InSAR), we present mean velocity maps of ground surface displacement in selected areas of Istanbul. We identify several sites along the terrestrial and coastal regions of Istanbul that underwent vertical ground subsidence at varying rates, from 5 ± 1.2 mm/yr to 15 ± 2.1 mm/yr. The results reveal that the most distinctive subsidence patterns are associated with both anthropogenic factors and relatively weak lithologies along the Haramirede valley in particular, where the observed subsidence is up to 10 ± 2 mm/yr. We show that subsidence has been occurring along the Ayamama river stream at a rate of up to 10 ± 1.8 mm/yr since 1992, and has also been slowing down over time following the restoration of the river and stream system. We also identify subsidence at a rate of 8 ± 1.2 mm/yr along the coastal region of Istanbul, which we associate with land reclamation, as well as a very localised subsidence at a rate of 15 ± 2.3 mm/yr starting in 2016 around one of the highest skyscrapers of Istanbul, which was built in 2010.


Introduction
Very rapid social and economic transformation in recent decades caused a huge rural-to-urban migration all over the world, which has fueled urban growth.This massive population shift has brought many complex challenges together with regard to sustainable development and natural disaster management.Istanbul is the largest city in Turkey, with a population of approximately 14 million inhabitants, and one of the most rapidly growing cities in Europe [1].According to the Istanbul Transportation Master Plan (ITMP), when taking into account the consequence of Turkey's economic growth in the last two decades and the large amount of immigration, projections indicate that the population will overcome 20 million inhabitants in 2023 [2].This rapid population growth poses major threats to the city's development when considering its vulnerability to natural disasters such as earthquakes, landslides, and floods, due to heavy and unplanned urbanization practices.Besides, the short distance (~20 km) of the main active branch of the North Anatolian Fault to the city poses a major threat to Istanbul [3].
Spaceborne interferometric synthetic aperture radar (InSAR) is a powerful remote sensing tool that enables observations of Earth's surface day and night under all weather conditions with high precision.Over the past decades, the method has been widely exploited in order to measure and monitor ground deformation induced by earthquakes [4][5][6][7] volcanoes [8], the withdrawal of ground water or other fluids [9,10], soil consolidation [11,12], mining [13], landslides [14,15], permafrost melting [16], ground subsidence [17,18], land reclamation [19], and sinkholes [20].Previous studies have shown the capacity of InSAR methods to measure and map land subsidence due to various anthropogenic factors, including ground water extraction in megacities such as Tianjin [21,22], Shanghai [23], Mexico City [24], and lithological factors in Bandung basin, Indonesia [25].Among InSAR techniques, persistent scatterers InSAR (PS-InSAR) was developed to tackle limitations related to temporal and geometrical decorrelation and atmospheric effects [26][27][28][29][30][31].It enables monitoring the temporal evolution of the ground motion by exploiting multiple SAR images acquired over the same area.It uses the radar return signal reflected from persistent scatterers (PS, pointwise phase-stable targets) such as rooftops, large rock outcrops, bridges, or motorways, where the spatial density of such PS is high [27].PS-InSAR analyses provide a time series of PS displacements and average surface velocities by searching a motion model that is relative to a reference point, and assumed to be motionless.
In the megacity of Istanbul, PS-InSAR time series analyses have allowed the monitoring of ground motions induced by anthropogenic activities [32], lithological features [33], and tectonic activities [7,34].However, there has been little discussion about the long-term temporal evolution of the ground motion and its possible causes.This study presents a PS-InSAR analysis of the secular ground motion in the urbanized metropolitan area of Istanbul.The processed InSAR data (Figure 1) span nearly 25 years, from 1992 to 2017, with multi-sensor images acquired along ascending and descending orbits.Most of the surface motion anomalies that we identify are associated with ground subsidence that has been induced by various factors, including natural compaction, and anthropogenic activities.These subsidence anomalies are carefully measured and analyzed from the perspective of urbanization and the assessment of geohazards for the city of Istanbul.The causes of subsidence in cities are diverse, and include factors such as lithology (i.e., rock type), variations in soil moisture content, groundwater exploitation, and overburden loads associated with human activity.In Istanbul, considering the proximity of several segments of the active North Anatolian Fault (NAF) in the Sea of Marmara, which have remained unbroken since 1776, the characterisation of subsidence susceptibility for Istanbul is crucial with regard to hazard mitigation and urban planning, as it can identify the vulnerable parts of the region that are prone to possible future earthquake damage.Thus, the main goal of the present study was to use a long-term PS-InSAR time series to: (1) quantify subsidence phenomena and discuss associated causes such as lithology-controlled natural compaction and anthropogenic activities, and (2) monitor the temporal evolution of the subsidence.

Background of the Study Areas
The present study primarily focusses on six areas where an InSAR time series enabled the detection of anomalous ground displacement: the region of Haramidere, where a series of active natural landslides had been previously mapped (circle 1 in Figures 1 and 2), and the Ayamama floodplain (circle 2 in Figure 1), a geological setting that is discussed in Section 2.1 (Figure 2), as well as several local subsidence areas that are related to anthropogenic activities (circles 3 to 6 in Figure 1).

Geological Setting of Study Areas 1 and 2
The Haramidere and the Ayamama streams are located near the boundary between the Istranca Methamorphic Unit (Paleozoic-Mesozoic) and Istanbul Unit (Paleozoic).These units are covered by the Eocene sedimentary sequence of the Thrace basin [36].The Paleozoic metamorphic basement of the study area consists in the east in a thick turbiditic sandstone-shale sequence of the Carboniferous age, while the Eocene cover is made of limestone, marl, and claystone units, which are transgressive on Çatalca metamorphic units in the west (Figure 2).

Background of the Study Areas
The present study primarily focusses on six areas where an InSAR time series enabled the detection of anomalous ground displacement: the region of Haramidere, where a series of active natural landslides had been previously mapped (circle 1 in Figures 1 and 2), and the Ayamama floodplain (circle 2 in Figure 1), a geological setting that is discussed in Section 2.1 (Figure 2), as well as several local subsidence areas that are related to anthropogenic activities (circles 3 to 6 in Figure 1).

Geological Setting of Study Areas 1 and 2
The Haramidere and the Ayamama streams are located near the boundary between the Istranca Methamorphic Unit (Paleozoic-Mesozoic) and Istanbul Unit (Paleozoic).These units are covered by the Eocene sedimentary sequence of the Thrace basin [36].The Paleozoic metamorphic basement of the study area consists in the east in a thick turbiditic sandstone-shale sequence of the Carboniferous age, while the Eocene cover is made of limestone, marl, and claystone units, which are transgressive on Çatalca metamorphic units in the west (Figure 2).

Study Area 1: Haramidere Valley and Avcilar Neighborhood
The Haramidere valley is located between Buyukcekmece Lake and Kucukcekmece Lake in the Avcilar Peninsula.It is located about 15 km north of the NAF, which cuts across the Sea of Marmara (Figure 2).Although the Avcilar district is located at a distance of 120 km west of the epicenter of the 17 August 1999 Izmit earthquake, it was the only area in the Istanbul metropolitan region that suffered extensive damage [40].During this earthquake, 27 buildings collapsed, 2076 other buildings were heavily damaged, and 273 casualties were reported in Avcilar [41,42].The maximum ground acceleration that was measured on soft sediments was 0.25 g, which is six times higher than the peak ground acceleration recorded on the bedrock in the center of Istanbul [40].This difference is the result of the amplification of seismic waves in surficial layers with soft lithology [40,[42][43][44][45][46].Despite a low background seismicity, which suggested no active faulting in the region [47], destructive and widespread damage during the 1999 Izmit earthquake drew considerable attention on this area, requiring the reassessment of active faulting.In order to investigate any relationship between the faults and damage observed in Avcilar and the vicinity, fault-mapping was refined from field studies and seismic reflection analyses [37, [48][49][50].These studies showed that the presence of secondary faults might be an important driving factor for the localised seismic amplification [7,51].
Landslides have been identified as another important geohazard in the suburb of Avcilar for many years [38,52,53].Recent events in that area often result from the reactivation of old landslides due to the overloading of the existing landslides by new constructions.The investigation of the old landslides in more detail is very important in order to anticipate measures that could be taken to avoid future possible damage in the urbanized center of the Avcilar peninsula.A total of 391 landslides were mapped (Figure 2) in the region [38].Approximately half of all of the landslides were distributed between Buyukcekmece and the Haramidere valley, which are important local landforms in Figure 2. Simplified geological and structural map of study areas 1 and 2 (circles 1 and 2 as in Figure 1).Numerous active landslides (dark green patches) were mapped between the Küçükçekmece and Büyükçekmece lakes [simplified from Ergintav et al., .The inset map shows a figure area north of the Sea of Marmara, and the main segments of the NAF in red [35].

Study Area 1: Haramidere Valley and Avcilar Neighborhood
The Haramidere valley is located between Buyukcekmece Lake and Kucukcekmece Lake in the Avcilar Peninsula.It is located about 15 km north of the NAF, which cuts across the Sea of Marmara (Figure 2).Although the Avcilar district is located at a distance of 120 km west of the epicenter of the 17 August 1999 Izmit earthquake, it was the only area in the Istanbul metropolitan region that suffered extensive damage [40].During this earthquake, 27 buildings collapsed, 2076 other buildings were heavily damaged, and 273 casualties were reported in Avcilar [41,42].The maximum ground acceleration that was measured on soft sediments was 0.25 g, which is six times higher than the peak ground acceleration recorded on the bedrock in the center of Istanbul [40].This difference is the result of the amplification of seismic waves in surficial layers with soft lithology [40,[42][43][44][45][46].Despite a low background seismicity, which suggested no active faulting in the region [47], destructive and widespread damage during the 1999 Izmit earthquake drew considerable attention on this area, requiring the reassessment of active faulting.In order to investigate any relationship between the faults and damage observed in Avcilar and the vicinity, fault-mapping was refined from field studies and seismic reflection analyses [37, [48][49][50].These studies showed that the presence of secondary faults might be an important driving factor for the localised seismic amplification [7,51].
Landslides have been identified as another important geohazard in the suburb of Avcilar for many years [38,52,53].Recent events in that area often result from the reactivation of old landslides due to the overloading of the existing landslides by new constructions.The investigation of the old landslides in more detail is very important in order to anticipate measures that could be taken to avoid future possible damage in the urbanized center of the Avcilar peninsula.A total of 391 landslides were mapped (Figure 2) in the region [38].Approximately half of all of the landslides were distributed between Buyukcekmece and the Haramidere valley, which are important local landforms in the region (Figure 2).Duman et al. [52,53] used field geological observations to argue that the parameters that control landslide initiation are shallow groundwater levels, lithology, and liquefaction.

Study Area 2: Ayamama River Valley
The Ayamama valley is located in the western part of Istanbul, east of the Haramidere valley (Figure 2).The river flows north-south from the eastern part of the Basaksehir district, and towards the Sea of Marmara in the Bakirkoy district [54].It runs through the heavily urbanised and highly populated area of the European side of Istanbul.The lower parts of its basin show various land use patterns and a high density of population [55].The Ayamama stream is known to produce seasonal flooding.According to the municipality of Istanbul, sedimentation and illegal urbanisation in the riverbeds have decreased floodplain capacity, which has subsequently caused periodic floods and overflows.One of the worst flooding events that affected the region was on 9 September 2009, which caused 31 casualties and 50 injuries.In the city's development plan, the Ayamama river basin and its surrounding zone had been set aside for recreational areas.However, after an amendment was made to develop it into residential area in 1997, industrial and residential land use rapidly increased along the axis of the stream [56].

Study Areas 3 to 6: Subsidence Events Associated with Land Reclamation and Urbanisation
Istanbul lies on both sides of the Bosphorus Strait ( İstanbul Bogazı), and has been subjected to heavy and unplanned urbanisation.This rapid urbanisation gave rise to land reclamation along the coastal areas of the Marmara Sea in order to provide new recreational areas to compensate for the destruction of green areas [57].Nearly 2.6 square kilometers of land have been gained from the Istanbul coast by filling up the sea since 2000 [58].We have focussed on the two recreational areas in the Yenikapi neighborhood and the Maltepe district for the investigation of the ground deformation related to land reclamation by the PS-InSAR method (sites 4 and 6).Another two local subsidence phenomena have been observed along the banks of the Golden Horn (Haliç, site 3), and around a skyscraper located in the Levent neighborhood of Istanbul (site 5).

Datasets
The SAR data used in the present study consisted of 291 C-band (5-6 GHz, ~6 cm wavelength) images acquired with various sensors, including 51 (ERS 1/2) images on two overlapping ascending tracks spanning from 1992 to 2001, 52 ENVISAT images on two overlapping ascending tracks spanning from 2003 to 2011, and 188 Sentinel 1A/B TOPS (Terrain Observation with Progressive Scans in azimuth) images in one ascending and two descending orbits acquired between 2014-2017 (Figure 3, see Table 1 for details).The Istanbul metropolitan area is entirely covered by all of the tracks.These multiorbit/sensor datasets enable us to examine the consistency between different sets of observations by inter-comparisons.For the Sentinel data, the period when the two satellites 1A/1B were operational is indicated in orange (before this period, only satellite 1A was operational).

Methodology
Single look complex images and interferograms of the Envisat/ASAR and ERS archives were generated using the ROI_PAC software [59] and the Delft Object-oriented Radar Interferometric Software (DORIS), respectively [60].All of the interferograms of the SENTINEL 1A/B satellites datasets were generated using the latest version of the software "Generic Mapping Tools Synthetic Aperture Radar (GMTSAR)" [61] using the Shuttle Radar Topography Mission (SRTM) 3-arcsecond digital elevation model for correcting topographic contributions to the radar phase.All of the interferograms were generated based on a single master network in order to use them in the PS-InSAR analysis.The choice of the master images for each track (red dots on Figure 3) were made so as to minimize the spatial and temporal baselines.The stacks of interferograms were processed with the PS-InSAR approach using the software package STAMPS (permanent scatterers InSAR analysis), which allows the identification of radar benchmarks by selecting pixels on the basis of their noise characteristics [62,63].For the selection of permanent scatterers, STAMPS takes into account the statistical relationship between amplitude and phase stability, which is quantified by the ratio between the standard deviation and the mean of the amplitude through time for each pixel (amplitude dispersion index).In the present study, we selected a threshold value of 0.42 for this amplitude dispersion index [27], which minimises the random amplitude variability, and eliminates highly decorrelated pixels in some areas covered with vegetation, agricultural fields, or snow.For the Sentinel data, the period when the two satellites 1A/1B were operational is indicated in orange (before this period, only satellite 1A was operational).

Methodology
Single look complex images and interferograms of the Envisat/ASAR and ERS archives were generated using the ROI_PAC software [59] and the Delft Object-oriented Radar Interferometric Software (DORIS), respectively [60].All of the interferograms of the SENTINEL 1A/B satellites datasets were generated using the latest version of the software "Generic Mapping Tools Synthetic Aperture Radar (GMTSAR)" [61] using the Shuttle Radar Topography Mission (SRTM) 3-arcsecond digital elevation model for correcting topographic contributions to the radar phase.All of the interferograms were generated based on a single master network in order to use them in the PS-InSAR analysis.The choice of the master images for each track (red dots on Figure 3) were made so as to minimize the spatial and temporal baselines.The stacks of interferograms were processed with the PS-InSAR approach using the software package STAMPS (permanent scatterers InSAR analysis), which allows the identification of radar benchmarks by selecting pixels on the basis of their noise characteristics [62,63].For the selection of permanent scatterers, STAMPS takes into account the statistical relationship between amplitude and phase stability, which is quantified by the ratio between the standard deviation and the mean of the amplitude through time for each pixel (amplitude dispersion index).In the present study, we selected a threshold value of 0.42 for this amplitude dispersion index [27], which minimises the random amplitude variability, and eliminates highly decorrelated pixels in some areas covered with vegetation, agricultural fields, or snow.

InSAR-Derived Land Subsidence Maps
Figure 4 shows the mean line-of-sight (LOS) velocity fields that were calculated from a PS-InSAR time series on the metropolitan city of Istanbul.A main deformation feature, which is common to all of the tracks, is the NNW-SSE elongated area in the southwest region of Istanbul, along the Haramidere valley on the Avcilar peninsula (Figure 4, circle 1 in a, site 1, as described in Sections 2.1 and 2.2).The associated motion is away from the satellite for both the descending and ascending tracks, suggesting a dominant subsidence signal, peaking at up to 10 mm/year in the line-of-sight.We identified another likely subsiding area along the Ayamama stream valley.This anomaly is elongated in a WNW-ESE direction, and is located along the western side of the Kucukcekmece Lake (Figure 4a circle 2, named as site 2).The area labeled with circle 3 in Figure 4a (named as site 3) covers the primary inlet of the Bosphorus, and is called Golden Horn urban waterway.This Golden Horn's bank is also likely subject to subsidence.Other potential subsidence anomalies caused by the settlement of the reclamation in the coastal areas were identified along the northern coast of the Marmara Sea (circles 4 and 6 in Figure 4a corresponding to sites 4 and 6).Lastly, we point to a local subsidence signal at a rather fast rate that was observed around a skyscraper (Figure 4a, circle and site 5) located in the Levent neighbourhood of Istanbul.The spatial and temporal variations of these different deformation patterns are discussed along with their possible underlying causes in the following sections.
Haramidere valley on the Avcilar peninsula (Figure 4, circle 1 in a, site 1, as described in Sections 2.1 and 2.2).The associated motion is away from the satellite for both the descending and ascending tracks, suggesting a dominant subsidence signal, peaking at up to 10 mm/year in the line-of-sight.We identified another likely subsiding area along the Ayamama stream valley.This anomaly is elongated in a WNW-ESE direction, and is located along the western side of the Kucukcekmece Lake (Figure 4a circle 2, named as site 2).The area labeled with circle 3 in Figure 4a (named as site 3) covers the primary inlet of the Bosphorus, and is called Golden Horn urban waterway.This Golden Horn's bank is also likely subject to subsidence.Other potential subsidence anomalies caused by the settlement of the reclamation in the coastal areas were identified along the northern coast of the Marmara Sea (circles 4 and 6 in Figure 4a corresponding to sites 4 and 6).Lastly, we point to a local subsidence signal at a rather fast rate that was observed around a skyscraper (Figure 4a, circle and site 5) located in the Levent neighbourhood of Istanbul.The spatial and temporal variations of these different deformation patterns are discussed along with their possible underlying causes in the following sections.

Self-Consistency between InSAR Measurements
We first quantitatively assessed the consistency of InSAR mean velocity measurements calculated from different datasets, across different time periods, for the widest areas of subsidence in the Haramidere and Avcilar valleys.An inter-comparison of the line-of-sight displacement rates in these areas was performed using a two-step procedure.In the first step, the mean velocity fields were downsampled by a factor of 10 in order to minimise errors within the geolocalisations (spatial mismatch between measurements).In a second step, we selected all of the pixels in the localised deforming areas of Haramidere Valley and Avcilar neighbourhood that were common to a pair of overlapping tracks, and compared the distribution of their subsidence rates for each pair of tracks (Figure 5).The consistency of velocities was good overall (mean correlation = 0.62).The observed differences could have originated from various factors, such as: (1) different incidence angles depending on the satellite, as they were not taken into account in the inter-comparison (Table 1), (2) different temporal coverages, (3) geolocalization uncertainty, (4) InSAR processing errors, and (5) seasonal effects [64].
We first quantitatively assessed the consistency of InSAR mean velocity measurements calculated from different datasets, across different time periods, for the widest areas of subsidence in the Haramidere and Avcilar valleys.An inter-comparison of the line-of-sight displacement rates in these areas was performed using a two-step procedure.In the first step, the mean velocity fields were downsampled by a factor of 10 in order to minimise errors within the geolocalisations (spatial mismatch between measurements).In a second step, we selected all of the pixels in the localised deforming areas of Haramidere Valley and Avcilar neighbourhood that were common to a pair of overlapping tracks, and compared the distribution of their subsidence rates for each pair of tracks (Figure 5).The consistency of velocities was good overall (mean correlation = 0.62).The observed differences could have originated from various factors, such as: (1) different incidence angles depending on the satellite, as they were not taken into account in the inter-comparison (Table 1), (2) different temporal coverages, (3) geolocalization uncertainty, (4) InSAR processing errors, and (5) seasonal effects [64].

Site 1: Haramidere Valley
Close-up views of the elongated pattern of subsidence in the Haramidere valley are shown on 6 for all tracks.In order to quantify the vertical, subsidence component of the motion, we decomposed the mean PS-InSAR line-of-sight velocity fields into east-west and vertical components using the method described by Samieie-Esfahany et al. [65].We only used the velocity fields calculated from the Sentinel 1A/B images (two tracks, 58 [Asc] and 138 [Dsc]) that covered the same time interval to

Site 1: Haramidere Valley
Close-up views of the elongated pattern of subsidence in the Haramidere valley are shown on 6 for all tracks.In order to quantify the vertical, subsidence component of the motion, we decomposed the mean PS-InSAR line-of-sight velocity fields into east-west and vertical components using the method described by Samieie-Esfahany et al. [65].We only used the velocity fields calculated from the Sentinel 1A/B images (two tracks, 58 [Asc] and 138 [Dsc]) that covered the same time interval to calculate this decomposition.We assumed that there was no north-south displacement, due to the low sensitivity of the LOS data to this component of displacement.Doing so, we reduced the number of unknown variables for each permanent scatterer point to two displacement components, d down , the vertical displacement (positive downward), and d west , the horizontal displacement in the east-west direction (positive toward west).In a first step, we resampled the mean line-of-sight velocities for the ascending and descending tracks onto a 0.0005 • × 0.0005 • regular grid (approximately 10 m spacing).We used the nearest neighbour procedure in the resampling of the persistent scatterer pixels that were within 30 m of the center of each grid nodal point.In a second step, we selected all of the pixels that exist in both the ascending and descending tracks.For further interpretation of the displacements, we referenced the two tracks using reference points located in an area that was assumed to be stable (circle in the NE part of Figure 6).In a third step, the decomposition of the line-of-sight velocity fields into east-west and vertical components was calculated, taking into account the local incidence angle of the satellite view (Figure 7).
low sensitivity of the LOS data to this component of displacement.Doing so, we reduced the number of unknown variables for each permanent scatterer point to two displacement components, ddown, the vertical displacement (positive downward), and dwest, the horizontal displacement in the east-west direction (positive toward west).In a first step, we resampled the mean line-of-sight velocities for the ascending and descending tracks onto a 0.0005° × 0.0005° regular grid (approximately 10 m spacing).We used the nearest neighbour procedure in the resampling of the persistent scatterer pixels that were within 30 m of the center of each grid nodal point.In a second step, we selected all of the pixels that exist in both the ascending and descending tracks.For further interpretation of the displacements, we referenced the two tracks using reference points located in an area that was assumed to be stable (circle in the NE part of Figure 6).In a third step, the decomposition of the lineof-sight velocity fields into east-west and vertical components was calculated, taking into account the local incidence angle of the satellite view (Figure 7).8).It is referenced to the mean value of the PS-InSAR points within the circle labeled R, which is considered a stable area.The DSC and ASC labels are for the descending and ascending orbits, respectively.
The vertical displacement field shows that subsidence occurs on both banks of the Haramidere valley, and follows an elongated area in the valley in a northwest-southeast direction (Figure 7c).The maximum subsidence is centered on the valley (Figure 7f).This region has a long, slow-moving landslides history, and is located in an area with shallow water level, poor soil conditions, and weak lithology, which are all parameters that are considered as favoring landslides [17,40,66].The subsidence that we observed coincided overall with previously mapped active landslide zones.However, the contours of these mapped landslides did not match precisely with the areas that had the highest subsidence in our vertical velocity map.The horizontal component in the east-west direction, which was derived by decomposition, showed a horizontal movement in the opposite direction on both banks of the valley, toward the valley center (Figures 7b,e).The sign change in the east-west velocity was also centered on the valley axis, as the maximum subsidence.We concluded The mean velocity value of the persistent scatterers (PS)-InSAR points within the solid black circle in the center of the maps has been used to illustrate the temporal evolution of the subsidence associated with weak lithology (Figure 8).It is referenced to the mean value of the PS-InSAR points within the circle labeled R, which is considered a stable area.The DSC and ASC labels are for the descending and ascending orbits, respectively.
The vertical displacement field shows that subsidence occurs on both banks of the Haramidere valley, and follows an elongated area in the valley in a northwest-southeast direction (Figure 7c).The maximum subsidence is centered on the valley (Figure 7f).This region has a long, slow-moving landslides history, and is located in an area with shallow water level, poor soil conditions, and weak lithology, which are all parameters that are considered as favoring landslides [17,40,66].The subsidence that we observed coincided overall with previously mapped active landslide zones.However, the contours of these mapped landslides did not match precisely with the areas that had the highest subsidence in our vertical velocity map.The horizontal component in the east-west direction, which was derived by decomposition, showed a horizontal movement in the opposite direction on both banks of the valley, toward the valley center (Figure 7b,e).The sign change in the east-west velocity was also centered on the valley axis, as the maximum subsidence.We concluded that the observed signal was consistent, rather, with the subsidence of the soil/surface on the valley slope (Figure 7d), moving downslope due to landslide or soil creep gravity.This downslope movement had both horizontal and vertical components, which were well captured by the InSAR data.
that the observed signal was consistent, rather, with the subsidence of the soil/surface on the valley slope (Figure 7d), moving downslope due to landslide or soil creep gravity.This downslope movement had both horizontal and vertical components, which were well captured by the InSAR data.In order to analyse the subsidence temporal evolution (Figure 8), we selected permanent scatterer points located in an area that was previously detected as undergoing landslide activity, and where the subsidence rate was among the highest observed in the present study.For the sake of consistency between the datasets acquired from different viewing geometries, these selected PS points were from an area where the horizontal velocity was considered negligible (see the circle in center of Figure 6), so that the line-of-sight velocities were converted into vertical velocities by a simple geometrical equation.The date of the first SAR image used here was taken as the reference time of the time series.As seen in Figure 8, the three datasets used for the Haramidere valley and Avcilar area had different starting dates and temporal coverage.For comparison, we set one reference time as 26 May 1992 for the three datasets, and the time series that were mapped from the Envisat and Sentinel datasets were shifted with a constant, which was calculated by assuming that the site was undergoing subsidence with a constant rate (Figure 8).
Although the subsidence rates that were calculated from the ERS and Envisat datasets were consistent with each other and matched well, the subsidence rates that were obtained from the Sentinel datasets were slightly smaller.Two reasons might lie behind this difference.The first reason is related to the relatively short duration of the SENTINEL 1 A/B time series, which could alter the accuracy of the rate estimate.The second reason may be the retardation of the settlement due to a long-term decay of the soil consolidation rate related to ground water extractions.Such an In order to analyse the subsidence temporal evolution (Figure 8), we selected permanent scatterer points located in an area that was previously detected as undergoing landslide activity, and where the subsidence rate was among the highest observed in the present study.For the sake of consistency between the datasets acquired from different viewing geometries, these selected PS points were from an area where the horizontal velocity was considered negligible (see the circle in center of Figure 6), so that the line-of-sight velocities were converted into vertical velocities by a simple geometrical equation.The date of the first SAR image used here was taken as the reference time of the time series.As seen in Figure 8, the three datasets used for the Haramidere valley and Avcilar area had different starting dates and temporal coverage.For comparison, we set one reference time as 26 May 1992 for the three datasets, and the time series that were mapped from the Envisat and Sentinel datasets were shifted with a constant, which was calculated by assuming that the site was undergoing subsidence with a constant rate (Figure 8).
Although the subsidence rates that were calculated from the ERS and Envisat datasets were consistent with each other and matched well, the subsidence rates that were obtained from the Sentinel datasets were slightly smaller.Two reasons might lie behind this difference.The first reason is related to the relatively short duration of the SENTINEL 1 A/B time series, which could alter the accuracy of the rate estimate.The second reason may be the retardation of the settlement due to a long-term decay of the soil consolidation rate related to ground water extractions.Such an exponential decay of ground subsidence was proposed to explain an InSAR time series on the Great Salt Lake in Utah [67].
Ground subsidence in the Avcilar peninsula has been previously reported at a mean rate of 6 mm/yr using an InSAR analysis of ERS-1 and ERS-2 satellite images taken between 1992-1999 [17], and at a rate of 10 mm/yr using ERS 1/2 and ENVISAT satellite images taken between 1992-2010 [64].These authors concluded that the spatial coverage of the land subsidence in this area, which was overall consistent with ours, was associated with partially saturated and unconsolidated shallow layers of soil formation with a relatively weak lithological profile.The results of the present study thus support the observation of land subsidence in Avcilar, at similar to higher rates.Such types of lithology-controlled subsidence have also been observed in the Bandung basin on the island of Java in Indonesia, using an ALOS-1 dataset, although at much higher rates (up to 12 cm/yr), along the boundaries between consolidated rocks and unconsolidated sediments [25].
Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 18 exponential decay of ground subsidence was proposed to explain an InSAR time series on the Great Salt Lake in Utah [67].
Ground subsidence in the Avcilar peninsula has been previously reported at a mean rate of 6 mm/yr using an InSAR analysis of ERS-1 and ERS-2 satellite images taken between 1992-1999 [17], and at a rate of 10 mm/yr using ERS 1/2 and ENVISAT satellite images taken between 1992-2010 [64].These authors concluded that the spatial coverage of the land subsidence in this area, which was overall consistent with ours, was associated with partially saturated and unconsolidated shallow layers of soil formation with a relatively weak lithological profile.The results of the present study thus support the observation of land subsidence in Avcilar, at similar to higher rates.Such types of lithology-controlled subsidence have also been observed in the Bandung basin on the island of Java in Indonesia, using an ALOS-1 dataset, although at much higher rates (up to 12 cm/yr), along the boundaries between consolidated rocks and unconsolidated sediments [25].

Site 2: Ayamama River
Another subsidence zone lies along the dense settlements of the Ayamama stream banks in the western part of Istanbul (site 2, Figure 4a).Here, the river has the typical morphology of a delta, and the streambed is mostly composed of young alluvial deposits, varying in thickness in the range 3-10 m [68].The subsidence rates we measured were about 6 mm/yr at maximum in the line-of-sight, corresponding to a maximum vertical subsidence rate of about 10 mm/yr (Figure 9).The area affected by subsidence shrunk gradually (Figure 9) during the observation period, which we interpreted as resulting from the flood prevention and remediation project of the river and stream system, which started in 2008.The subsidence pattern along the Ayamama stream was firstly reported by Walter et al. [33] with a more limited dataset, who suggested that the ongoing land subsidence in the region is related to the sediment consolidation process, and that there might be direct or indirect consequences of the destructive flood events, such as the one in 2009.This is thus consistent with our interpretation of the recovery of the region following the restoration of the river and stream system.(referenced to points in area labeled R in Figure 6).

Site 2: Ayamama River
Another subsidence zone lies along the dense settlements of the Ayamama stream banks in the western part of Istanbul (site 2, Figure 4a).Here, the river has the typical morphology of a delta, and the streambed is mostly composed of young alluvial deposits, varying in thickness in the range 3-10 m [68].The subsidence rates we measured were about 6 mm/yr at maximum in the line-of-sight, corresponding to a maximum vertical subsidence rate of about 10 mm/yr (Figure 9).The area affected by subsidence shrunk gradually (Figure 9) during the observation period, which we interpreted as resulting from the flood prevention and remediation project of the river and stream system, which started in 2008.The subsidence pattern along the Ayamama stream was firstly reported by Walter et al. [33] with a more limited dataset, who suggested that the ongoing land subsidence in the region is related to the sediment consolidation process, and that there might be direct or indirect consequences of the destructive flood events, such as the one in 2009.This is thus consistent with our interpretation of the recovery of the region following the restoration of the river and stream system.

Sites 4, 6, and 3
We also identified other sites along the terrestrial and coastal region of Istanbul undergoing ground subsidence.The subsidence that we related to land reclamation in two recreational areas in the Yenikapi neighbourhood (Figure 10a, site 4 in Figure 4a) and the Maltepe district (Figure 10d, site 6 in Figure 4a), which were constructed in 2014, was observed at a vertical rate of 10 mm/yr.This subsidence was presumably related to the compaction induced by the primary consolidation process of the alluvial clay beneath the reclamation zone.The subsidence rates that were measured in both reclamation areas were likely dependent on the physical characteristic and thickness of the underlying alluvial deposit and used matrix of the filling material [68,69].Another local subsidence area related to a similar phenomena was observed along the shores of the Golden Horn (Haliç) (Figure 10b, site 3 in Figure 4a).A significant part of this subsidence was located on reclamation areas that have been transformed into parks and recreational facilities along both banks of the Golden Horn.Besides that, the shorelines in this area have also undergone significant urban changes within the frame of the renovation project of all of the waterfronts in Istanbul.During the renovation, sediments made of loose clay deposits were removed from the shallowest parts of the Golden Horn, which might have triggered subsidence in the nearby waterfront areas.The subsidence pattern along the Golden Horn that we observed was consistent with the coastal subsidence that was previously described in the Figures 5 and 7 of [31], which, used a dataset of high-resolution TerraSAR-X SAR images covering the period 2010-2012.In their analyses of subsidence evolution over the urbanised

Sites 4, 6, and 3
We also identified other sites along the terrestrial and coastal region of Istanbul undergoing ground subsidence.The subsidence that we related to land reclamation in two recreational areas in the Yenikapi neighbourhood (Figure 10a, site 4 in Figure 4a) and the Maltepe district (Figure 10d, site 6 in Figure 4a), which were constructed in 2014, was observed at a vertical rate of 10 mm/yr.This subsidence was presumably related to the compaction induced by the primary consolidation process of the alluvial clay beneath the reclamation zone.The subsidence rates that were measured in both reclamation areas were likely dependent on the physical characteristic and thickness of the underlying alluvial deposit and used matrix of the filling material [68,69].Another local subsidence area related to a similar phenomena was observed along the shores of the Golden Horn (Haliç) (Figure 10b, site 3 in Figure 4a).A significant part of this subsidence was located on reclamation areas that have been transformed into parks and recreational facilities along both banks of the Golden Horn.Besides that, the shorelines in this area have also undergone significant urban changes within the frame of the renovation project of all of the waterfronts in Istanbul.During the renovation, sediments made of loose clay deposits were removed from the shallowest parts of the Golden Horn, which might have triggered subsidence in the nearby waterfront areas.The subsidence pattern along the Golden Horn that we observed was consistent with the coastal subsidence that was previously described in the Figures 5 and 7 of [31], which, used a dataset of high-resolution TerraSAR-X SAR images covering the period 2010-2012.In their analyses of subsidence evolution over the urbanised region of Istanbul, these authors similarly concluded that the patterns of settlement along the Golden Horn shores were caused by anthropogenic factors arising from unsustainable urban development.
Remote Sens. 2018, 10, x FOR PEER REVIEW 13 of 18 region of Istanbul, these authors similarly concluded that the patterns of settlement along the Golden Horn shores were caused by anthropogenic factors arising from unsustainable urban development.

Site 5
High skyscrapers might suffer rapid settlement and declination, and cause local ground surface subsidence due to the consolidation of the underlying soft soil deposits.In this study, we detected a very localised subsidence pattern around a skyscraper built in 2010, and other high-rise buildings in the Levent neighbourhood of Istanbul starting from 2016, with an average subsidence rate of about

Site 5
High skyscrapers might suffer rapid settlement and declination, and cause local ground surface subsidence due to the consolidation of the underlying soft soil deposits.In this study, we detected a very localised subsidence pattern around a skyscraper built in 2010, and other high-rise buildings in the Levent neighbourhood of Istanbul starting from 2016, with an average subsidence rate of about 15 mm/yr (Figure 10c).A time series analysis of the point targets surrounding the skyscraper showed a very rapid increase of subsidence rate during the first four months of 2016, which might be related to the groundwater lowering around the foundation of the building, which may cause a downward movement of the surface surrounding the building.

Conclusions
Istanbul has been subject to intense industrialisation and population increase, especially since the 1960s, which is causing very rapid urbanisation and heavy land-use changes.We identified several sites along the terrestrial and coastal region of Istanbul that have been undergoing vertical ground subsidence at rates ranging from 5 mm/yr to 15 mm/yr.In the present study, a PS-InSAR time series analysis was performed using 291 C-band SAR images in order to characterise these subsidence phenomena by combining multi-track/sensor InSAR datasets and provide insights into the potential hazards induced by local soil conditions and human activities.Using the PS-InSAR technique, enough time-coherent pixels were obtained over six different sites in Istanbul.The most extended, clearly visible subsidence signals were detected over the western part of Istanbul consistently in all of the frames.The spatiotemporal variability of the ground displacement was measured over the last 25 years in Avcilar, which suffered extensive damage during the 1999 Izmit earthquake and along the Haramidere valley, where a long-lasting landslide history has been reported.The time series analysis in this region revealed that the Haramidere valley banks are undergoing ground vertical subsidence at a rate of 10 ± 2 mm/yr.Another subsidence area has been reported along the Ayamama river banks, which is made of shallow alluvial deposits due to the local substratum consolidation process, at a rate of 10 ± 1.8 mm/yr, and the surface area affected by subsidence has been shrinking gradually following the restoration of the river and stream system, which started in 2008.The inter-comparison of PS-InSAR measurements from different satellite sensors for the western part of Istanbul showed that the correlation between the mean velocity fields were high between 1992-2011, and decreased with time due to spatiotemporal changes in the ground deformation.Sentinel-1A/B InSAR measurements acquired between 2014-2017 showed that reclaimed lands in both the European (Yenikapi reclamation area) and Asian (Maltepe reclamation area) coastlines of Istanbul underwent significant subsidence of up to 8 ± 1.3 mm/yr as a result of the primary consolidation process of the alluvial clay beneath the filling material.Lastly, a very localised subsidence pattern was detected around a skyscraper, with average subsidence rate of 15 ± 1.2 mm/yr.
On the whole, we can conclude that during the interseismic period, human-driven changes produced a more significant control on the coastal subsidence in Istanbul than natural factors.In future studies, such high-resolution SAR data over the very dense urban areas of Istanbul could help continuously monitor the urbanised areas suffering from land subsidence.With this purpose, future plans include the processing of high-resolution X-band sensors COSMO-Skymed and TerraSAR-X datasets in order to better quantify the surface settlement and reveal the underlying causes of these settlements, and thus provide more complete data to public organisations that are in charge of sustainable urban policies and hazard mitigation.

Figure 1 .
Figure 1.Study area and satellite synthetic aperture radar data coverage used in the present study.The shaded topography is given by the Shuttle Radar Topography Mission (SRTM) along the North Anatolian Fault (NAF) in the Sea of Marmara, and major faults are drawn in red [35].Rectangles labeled with sensor and track numbers indicate the coverage of the SAR images that were used in the present study.The red and black arrows indicate satellite's line-of-sight look and flight directions, respectively.Circles with numbers show the study regions in the paper (details are given in the text).

Figure 1 .
Figure 1.Study area and satellite synthetic aperture radar data coverage used in the present study.The shaded topography is given by the Shuttle Radar Topography Mission (SRTM) along the North Anatolian Fault (NAF) in the Sea of Marmara, and major faults are drawn in red [35].Rectangles labeled with sensor and track numbers indicate the coverage of the SAR images that were used in the present study.The red and black arrows indicate satellite's line-of-sight look and flight directions, respectively.Circles with numbers show the study regions in the paper (details are given in the text).

Figure 2 .
Figure 2. Simplified geological and structural map of study areas 1 and 2 (circles 1 and 2 as in Figure 1).Numerous active landslides (dark green patches) were mapped between the Küçükçekmece and Büyükçekmece lakes [simplified from Ergintav et al., Duman et al. and Ozgul et al. [37-39].The inset map shows a figure area north of the Sea of Marmara, and the main segments of the NAF in red [35].

Figure 3 .
Figure 3. Baseline versus time plots for the seven tracks used in this study.The red dots indicate the master image used as a reference for each track.For the Sentinel data, the period when the two satellites 1A/1B were operational is indicated in orange (before this period, only satellite 1A was operational).

Figure 3 .
Figure 3. Baseline versus time plots for the seven tracks used in this study.The red dots indicate the master image used as a reference for each track.For the Sentinel data, the period when the two satellites 1A/1B were operational is indicated in orange (before this period, only satellite 1A was operational).

Figure 4 .
Figure 4. Averaged line-of-sight velocity maps of the Istanbul metropolitan area from an interferometric synthetic aperture radar (InSAR) time series analysis, with varying time spans depending on the sensor.Negative velocities (cold colors) represent the displacement of the ground toward the satellite, and positive velocities (warm colors) indicate the displacement away from the satellite.Red lines in the Sea of Marmara indicate the submarine branches of the NAF.Average lineof-sight velocity (a) for Sentinel ascending track 58.The solid black circles labeled from 1 to 6 indicate the locations of the subsidence anomalies that are discussed in the present study.1: Haramidere River, 2: Ayamama Stream Valley, 3: Golden Horn, 4: Yenikapi reclamation area, 5: Skyscrapers in Levent neighbourhood, 6: Maltepe reclamation area; (b) for SENTINEL 1A/B descending track 36; (c) for Sentinel descending track 138; (d) for Envisat descending track 107; (e) for ERS descending track 336; (f) for ENVISAT descending track 336 and (g) for ERS descending track 107.

Figure 4 .
Figure 4. Averaged line-of-sight velocity maps of the Istanbul metropolitan area from an interferometric synthetic aperture radar (InSAR) time series analysis, with varying time spans depending on the sensor.Negative velocities (cold colors) represent the displacement of the ground toward the satellite, and positive velocities (warm colors) indicate the displacement away from the satellite.Red lines in the Sea of Marmara indicate the submarine branches of the NAF.Average line-of-sight velocity (a) for Sentinel ascending track 58.The solid black circles labeled from 1 to 6 indicate the locations of the subsidence anomalies that are discussed in the present study.1: Haramidere River, 2: Ayamama Stream Valley, 3: Golden Horn, 4: Yenikapi reclamation area, 5: Skyscrapers in Levent neighbourhood, 6: Maltepe reclamation area; (b) for SENTINEL 1A/B descending track 36; (c) for Sentinel descending track 138; (d) for Envisat descending track 107; (e) for ERS descending track 336; (f) for ENVISAT descending track 336 and (g) for ERS descending track 107.

Figure 5 .
Figure 5. Quantitative comparison of the line-of-sight displacement rates between all of the tracks used in the present study.The upper-right triangle matrix shows pairwise correlations of different tracks with correlation values and color intensities (blue and red indicate low and high correlation, respectively).In the lower-left triangle, the black dots denote the points that can be extracted on both tracks.SEN and ENV in the panel represent SENTINEL and ENVISAT, respectively.

Figure 5 .
Figure 5. Quantitative comparison of the line-of-sight displacement rates between all of the tracks used in the present study.The upper-right triangle matrix shows pairwise correlations of different tracks with correlation values and color intensities (blue and red indicate low and high correlation, respectively).In the lower-left triangle, the black dots denote the points that can be extracted on both tracks.SEN and ENV in the panel represent SENTINEL and ENVISAT, respectively.

Figure 6 .
Figure 6.Zoom views of displacement rates in the Haramidere Valley and its surrounding area.The mean velocity value of the persistent scatterers (PS)-InSAR points within the solid black circle in the center of the maps has been used to illustrate the temporal evolution of the subsidence associated with weak lithology (Figure8).It is referenced to the mean value of the PS-InSAR points within the circle labeled R, which is considered a stable area.The DSC and ASC labels are for the descending and ascending orbits, respectively.

Figure 6 .
Figure 6.Zoom views of displacement rates in the Haramidere Valley and its surrounding area.The mean velocity value of the persistent scatterers (PS)-InSAR points within the solid black circle in the center of the maps has been used to illustrate the temporal evolution of the subsidence associated with weak lithology (Figure8).It is referenced to the mean value of the PS-InSAR points within the circle labeled R, which is considered a stable area.The DSC and ASC labels are for the descending and ascending orbits, respectively.

Figure 7 .
Figure 7. Decomposition of horizontal and vertical components of ground displacement using only S-1 datasets.(a) The shaded topography was given by the Shuttle Radar Topography Mission (SRTM) along the Avcilar region.Fault lines were simplified from Ergintav et al. [37]; (b) Vertical component.Patches with thick dark boundaries correspond to the landslides that were identified geological maps, as shown on Figure 2 [simplified from Duman et al. and Ozgul et al. [38,39]; (c) Horizontal component in the east-west direction.(d)Valley-perpendicular elevation profile extracted from (a); (e,f) are horizontal and vertical velocity profiles extracted from (b,c) respectively.

Figure 7 .
Figure 7. Decomposition of horizontal and vertical components of ground displacement using only S-1 datasets.(a) The shaded topography was given by the Shuttle Radar Topography Mission (SRTM) along the Avcilar region.Fault lines were simplified from Ergintav et al. [37]; (b) Vertical component.Patches with thick dark boundaries correspond to the landslides that were identified geological maps, as shown on Figure 2 [simplified from Duman et al. and Ozgul et al. [38,39]; (c) Horizontal component in the east-west direction.(d)Valley-perpendicular elevation profile extracted from (a); (e,f) are horizontal and vertical velocity profiles extracted from (b,c) respectively.

Figure 8 .
Figure 8.Time series of the vertical displacement at the selected PS points circled in center of Figure6(referenced to points in area labeled R in Figure6).

Figure 8 .
Figure 8.Time series of the vertical displacement at the selected PS points circled in center of Figure6(referenced to points in area labeled R in Figure6).

Figure 9 .
Figure 9. Spatiotemporal characteristics of subsidence along the Ayamama river valley study area.(a) Mean displacement rates in line-of-sight (LOS) obtained from a track 336 ERS dataset.The black dashed lines indicate two profiles, one in an area with an active subsidence (profile labeled 1-4), and the other one with the same length in an area considered as stable and used as a reference (profile labeled 1′-4′).The inset map indicates the temporal and spatial pattern of subsidence for the region, from 1992 to 2017; (b-h) Rates along the profiles 1-4 (black) and 1′-4′ (red) taken for each track; (i) Temporal evolution of the coastal subsidence of selected points around point 3 in Figure 9a.

Figure 9 .
Figure 9. Spatiotemporal characteristics of subsidence along the Ayamama river valley study area.(a) Mean displacement rates in line-of-sight (LOS) obtained from a track 336 ERS dataset.The black dashed lines indicate two profiles, one in an area with an active subsidence (profile labeled 1-4), and the other one with the same length in an area considered as stable and used as a reference (profile labeled 1 -4 ).The inset map indicates the temporal and spatial pattern of subsidence for the region, from 1992 to 2017; (b-h) Rates along the profiles 1-4 (black) and 1 -4 (red) taken for each track; (i) Temporal evolution of the coastal subsidence of selected points around point 3 in Figure 9a.

Figure 10 .
Figure 10.Vertical velocities obtained by the decomposition of mean velocity fields of Sentinel 1 data (T58 ascending track and T138 descending track) superimposed on a Google Earth image of Istanbul, and the relevant time series of the vertical displacement.Black, red, and blue triangles represent the ascending T58, descending T36, and descending T138 tracks, respectively.(a) Yenikapi coastal and land reclamation area (circle 4 in Figure 4a).The color scale represents the vertical displacement of the surface; (b) Golden Horn area (circle 3 in Figure 4a); (c) Highly urbanised area of Istanbul, with subsiding persistent scatterer points clustered around the highest skyscraper of Istanbul (circle 5 in Figure 4a); (d) Maltepe reclamation zone (circle 6 in Figure 4a).

Figure 10 .
Figure 10.Vertical velocities obtained by the decomposition of mean velocity fields of Sentinel 1 data (T58 ascending track and T138 descending track) superimposed on a Google Earth image of Istanbul, and the relevant time series of the vertical displacement.Black, red, and blue triangles represent the ascending T58, descending T36, and descending T138 tracks, respectively.(a) Yenikapi coastal and land reclamation area (circle 4 in Figure 4a).The color scale represents the vertical displacement of the surface; (b) Golden Horn area (circle 3 in Figure 4a); (c) Highly urbanised area of Istanbul, with subsiding persistent scatterer points clustered around the highest skyscraper of Istanbul (circle 5 in Figure 4a); (d) Maltepe reclamation zone (circle 6 in Figure 4a).

Table 1 .
Characteristics of each processed track.

Table 1 .
Characteristics of each processed track.
1 SENTINEL 1 A/B TOPS; 2 Density of permanent scatterers (PS) in the urban areas.