Investigating Subsidence in the Bursa Plain, Turkey, Using Ascending and Descending Sentinel-1 Satellite Data

Abstract

We characterize and monitor subsidence of the Bursa Plain (southern Marmara region of Turkey), which has been interpreted as resulting from tectonic motions in the region. We quantify the subsidence using Interferometric Synthetic Aperture Radar (InSAR) time-series analysis. The Stanford Method for Persistent Scatterers InSAR package (StaMPS) is employed to process series of Sentinel 1 A-B radar images acquired between 2014 and 2017 along both ascending and descending orbits. The vertical velocity field obtained after decomposition of line-of-sight velocity fields on the two tracks reveals that the Bursa plain is subsiding at rates up to 25 mm/yr. The most prominent subsidence signal in the basin forms an east-west elongated ellipse of deformation in the east, and is bounded by a Quaternary alluvial plain undergoing average vertical subsidence at ~10 mm/yr. Another localized subsidence signal is located 5 km north of the city, following the Bursa alluvial fan, and is subsiding at velocities up to 25 mm/yr. The comparison between temporal variations of the subsiding surface displacements and variations of the water pressure head in the aquifer allows estimation of the compressibility of the aquifer, $\alpha$. It falls in the range of $0.5\times {10}^{-6}-2\times {10}^{-6}$ Pa⁻¹, which corresponds to typical values for clay and sand sediments. We find a clear correlation between subsidence patterns and the lithology, suggesting a strong lithological control over subsidence. In addition, the maximum rate of ground subsidence occurs where agricultural activity relies on groundwater exploitation. The InSAR time series within the observation period is well correlated with changes in the depth of the ground water. These observations indicate that the recent acceleration of subsidence is mainly due to anthropogenic activities rather than tectonic motion.

1. Introduction

Subsidence is one of the most diverse forms of ground vertical displacement, varying from local collapses to regional-scale sinking. Depending on the deformation regimes of the subsidence, vertical displacement may be associated with horizontal deformation that may have significant damaging effects. Land subsidence in the Bursa Plain in northwestern Turkey was first documented using four ALOS PALSAR (L-band) images acquired between 2007 and 2010 [1]. The authors investigated the surface deformation using the differential Interferometric Synthetic Aperture Radar (InSAR) technique and concluded that Bursa plain is subsiding at a rate of 10.9 cm/yr. They suggested that this subsidence is driven by local tectonic activity, as the Bursa basin likely corresponds to a pull-apart basin geometry, caused by strike-slip faulting. However, it is well known that subsidence may also occur due to the compaction of unconsolidated sediments that results from exploitation of groundwater [2,3,4] and due to urbanization [5,6]. It is therefore important to determine the mechanism of subsidence and its spatiotemporal characteristics because it allows the development of mitigation strategies that could minimize the risk of damage to surface structures and infrastructures [7,8,9,10].

Ground subsidence is a common phenomenon in urban areas built on thick, unconsolidated loose soil layers consisting of clay, silt, and peat formations [11,12,13], as is the case in the Bursa basin. Subsidence has been correlated to strata lithology, formation, and structure [14,15,16,17]. Dewatering of loose soil layers in which the water provides structural support produces volumetric contraction and surface subsidence [18]. Consequently, the depth of groundwater exploitation strongly controls the magnitude of land subsidence.

Since the 1960s, a considerable amount of literature has been published on land subsidence using various methods of investigations, including conventional practices and advanced ground sensing systems, depending on the spatial extent of the anticipated deformation. Data are interpreted to characterize the deformation mechanism and their environmental consequences [19]. These conventional techniques include optical leveling [20], GPS (Global Positioning System) surveys [21,22], groundwater monitoring [23], and tape-extensometers [24,25,26]. Advanced techniques include airborne and space based observtions, such as InSAR [27,28] and LiDAR [29], which can provide high-resolution imagery of the subsiding rates. InSAR is an observation technique that uses phase differences in reflected radar signals acquired at different times to detect and monitor different forms of ground deformations over wide areas (up to 400 × 400 km) with potentially high resolution and precision depending on the sensors [30,31,32,33]. Due to the development of new satellite technologies and advancing processing techniques, InSAR has been successfully used to investigate deformation resulting from subsidence due to the withdrawal of groundwater or other fluids [4,34,35,36], soil consolidation [37,38], subsurface mining or collapse of old mines [39,40], permafrost thawing [41], tectonic activity [42], and soluble earth materials [43]. Besides, InSAR has also been used to infer aquifer hydraulic properties, such as compressibility and the storage coefficient, as well as the mechanics of aquifer systems [2,4,44,45,46,47].

The present study focuses on the surface deformation associated with groundwater extraction, and its temporal evolution in the Bursa Plain to understand how surface subsidence evolves in response to the water table fluctuations and changing soil properties. We investigate this evolution by integrating SAR interferometry and geological and hydrogeological data to shed insight on the underlying processes governing subsidence. We have restricted our study to the Bursa Plain located in the east of the Bursa basin along the northern slope of the Uludag Mount (Figure 1 and Figure 2a). We use PS-InSAR time series with Sentinel-1 data acquired between 2014 and 2017. The observed subsidence signal is interpreted by comparison with spatio-temporal piezometric level variations in the aquifer, and analyzed in terms of local geological conditions. We first describe the tectonic and geological setting of the study area, and other morphological characteristics. After summarizing the Sentinel-1 SAR data used and the processing methodologies employed, we present the mean velocity fields obtained from different viewing geometries, and their vertical components. Finally, we compare vertical InSAR velocity fields with hydrological and geological data to constrain the compressibility of the underlying aquifer that compacted due to a reduction of the water level decrease. In addition, we interpret Landsat 8 optical images to better understand the relationship between subsidence and the spatial pattern of alluvial fans, bare lands, urban and agricultural areas, and other land-use patterns.

2. Study Area

Located on the northern slope of Mount Uludag (2543 m a.s.l.) in north-west Turkey, the city of Bursa is one of the most industrialized and populated cities of Turkey, with a population of 2.9 million inhabitants [48]. The Bursa plain is located 5 km to the north of the city and bordered by the Mount Uludag to the south and Mount Katirli to the north (Figure 1b). The plain has a high agricultural potential due to its rich underground aquifer and surface water resources, and its fertile soils [49]. Consequently, the plain exerts considerable impact upon Turkey’s national economy, which has led to rapid population inflow over the past few decades. This determined a rapid expansion of industrial and residential zones into the fertile Bursa plain, previously mainly for farming activities.

From a tectonic point of view, the southern branch of the North Anatolian Fault (NAF) in the Marmara region controls the long-term tectonic evolution of the Bursa basin. The segmentation of the NAF in the region results in the formation of extentional basins (with a roughly N-S extension direction), at the origin of the Bursa Plain in the eastern part of the E-W trending Bursa-Gönen Depression, named Bursa basin hereafter. This basin is bounded to the north and the south by the Bandirma-Mudanya and Uludag-Sularya mountain ranges, respectively (Figure 1b) [50]. These highlands consist of metamorphic and granitic basement rocks and are dominant morphological features along the southern branch of the NAF [51]. Tectonic models suggest that the depressed areas opened as pull-apart basins and are filled with Upper Pleistocene fluvial terrace and alluvial fan materials (Figure 2b) [50]. The alluvium has high transmissivity in comparison with the Neogene deposits surrounding the Nilufer Valley [52]

The city of Bursa is located within the area with the highest probability of strong earthquakes in Turkey (i.e., first-degree seismic hazard zone). The seismicity of the region is controlled mainly by active faults, such as the Gemlik Fault, Geyve-Iznik Fault Yenisehir Fault, Bursa Fault, and the Inonu-Eskisehir Fault Zone [53]. Bursa also hosts thermal springs that follow major fracture zones [54]. The Quaternary sediments thicken to the south and exceeds 300 m in the Bursa basin [55,56], suggesting that the deposition pattern is tectonically controlled [50].

3. Datasets and Methodology

3.1. InSAR Datasets

We use C-band (~6 cm wavelength) Sentinel 1 A/B TOPS images to map the spatial distribution and rate of subsidence in the Bursa Basin. The data used in the present study consist of two datasets with short temporal sampling (Figure 3). 105 images were acquired on ascending Track 58 between 21 October 2014 and 28 December 2017, and 96 images were acquired on descending Track 138 between 8 November 2014 and 4 December 2017. The study area is fully covered by all tracks (Figure 1b).

3.2. Methodology

3.2.1. PS-InSAR Processing Methodology

All interferograms were generated using the latest version of the open source (GNU General Public License) InSAR processing system called “Generic Mapping Tools Synthetic Aperture Radar (GMTSAR)” [61]. To correct the topographic contribution to the radar phase, we used the Shuttle Radar Topography Mission (SRTM) 3-arcsecond digital elevation model [62]. All interferograms were computed based on a single master network for Permanent Scatterer InSAR (PS-InSAR) analysis. The choice of the master images for both tracks (blue dots on Figure 3) minimizes the spatial and temporal baselines. The single master stacks of interferograms were processed using the software package, StaMPS [63,64]. StaMPS allows the identification of PS points, using both amplitude and phase information. In the first step, the initial selection of PS points is performed based on their noise characteristics, using amplitude analysis. The amplitude dispersion criterion is defined by ${\mathrm{D}}_{\mathrm{Amp}}=\left({\mathsf{\sigma }}_{\mathrm{Amp}}/{\mathrm{m}}_{\mathrm{Amp}}\right)$, where ${\mathsf{\sigma }}_{\mathrm{Amp}}$ and ${\mathrm{m}}_{\mathrm{Amp}}$ are the standard and mean deviation of the amplitude in time, respectively [65]. We selected a threshold value of 0.4 for ${\mathrm{D}}_{\mathrm{Amp}}$, which minimizes the random amplitude variability and eliminates highly decorrelated pixels in some areas covered with vegetation, agricultural fields, or snow. Once coherent PS points have been selected based on amplitude analysis, the PS probability is refined by phase analysis in a series of iterations. This process allows the detection of stable pixels even with low amplitude. Once the final selection of PS has been done, the residual topographic component can be removed. Then, 3D unwrapping of the PS phase is performed both spatially and temporally. This analysis enables retrieval of the average Line-Of-Sight (LOS) subsidence rate maps. To remove atmospheric effects, we used the freely available Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) developed by Bekaert et al. [66]. This toolbox uses ERA-Interim (ERA-I, European Center for Medium-Range Weather Forecasts) numerical weather model datasets [67].

The results obtained from the InSAR time series were compared with the groundwater level pressure head variations obtained from piezometric measurements. The water level data were collected in two wells (#58581 and #62821), whose locations are shown in Figure 2b, and used to demonstrate the correlation in time between land subsidence and the level of water in the aquifer. The piezometric heads in the study area were measured monthly during the InSAR observation time, from October 2014 to December 2017.

3.2.2. Decomposition of Sentinel-1 Data into 2D Displacement Rates

The subsidence signal in the Bursa plain is covered entirely by Sentinel-1 data from two different orbits and therefore only two different viewing geometries are available. Consequently, we can only retrieve two of the three components of the actual deformation vector. We decomposed the mean PS-InSAR line-of-sight velocity fields into east-west and vertical components. We assumed that the north-south displacement was negligible because the near polar orbits of the SAR satellites produce a low sensitivity to the north-south component of displacement [68]. In a first step, we resampled the mean line-of-sight velocities for the ascending and descending tracks onto a 100 × 100 m regular grid, which was used as input for displacement decomposition. We used a nearest neighbor procedure to resample the persistent scatterer pixels that are within 200 m of the center of each grid nodal point. In a second step, we selected all the pixels that exist in both the ascending and descending tracks. Before decomposition, we transformed the InSAR mean velocity fields of both tracks into the same reference frame using a reference area considered as a stable area (black star in Figure 2b and Figure 4). In a last step, the line-of-sight velocity fields were decomposed into two components: The horizontal component along the east-west direction (${\mathrm{d}}_{\mathrm{hor}}$) and the vertical component (${\mathrm{d}}_{\mathrm{ver}}$) were computed, considering the local incidence angle of the satellite view by solving the following equation [69]:

$$\left(\begin{array}{c}{\mathrm{d}}_{\mathrm{asc}}\\ {\mathrm{d}}_{\mathrm{dsc}}\end{array}\right)=\left(\begin{array}{cc}{\cos \mathsf{\theta }}_{\mathrm{asc}}& -{\cos \mathsf{\alpha }}_{\mathrm{asc}}{\sin \mathsf{\theta }}_{\mathrm{asc}}\\ {\cos \mathsf{\theta }}_{\mathrm{dsc}}& -{\cos \mathsf{\alpha }}_{\mathrm{dsc}}{\sin \mathsf{\theta }}_{\mathrm{dsc}}\end{array}\right)\left(\begin{array}{c}{\mathrm{d}}_{\mathrm{ver}}\\ {\mathrm{d}}_{\mathrm{hor}}\end{array}\right)$$

(1)

where ${\mathsf{\theta }}_{\mathrm{asc}}$ and ${\mathsf{\theta }}_{\mathrm{dsc}}$ represent the local incidence angles and ${\mathsf{\alpha }}_{\mathrm{asc}}$ and ${\mathsf{\alpha }}_{\mathrm{dsc}}$ are the satellite heading angles in the ascending and descending modes, respectively.

4. Results

InSAR-Derived Land Subsidence Maps

Figure 4 shows the InSAR-derived mean line-of-sight velocity fields calculated from the PS-InSAR time series analysis for the Bursa Plain. There is no ground-truth data (i.e., data collected at the ground surface) available in the study area. Consequently, the mean phase values for both tracks were initially used for referencing the velocity fields in the time series analysis. The area in the south of the subsiding basin located within the shaded polygon representing the extension of the Bursa city was chosen as the reference area because it is considered as a stable area (black star in Figure 3).

The velocity fields of both tracks exhibit similar spatial patterns. This similarity suggests that the LOS velocities show mostly vertical deformation with negligible horizontal motion and atmospheric effects. The most prominent deformation pattern in the eastern part of the basin, northeast of the Bursa city, is thus interpreted as subsidence. In order to assess the consistency between ascending and descending orbits, we performed an inter-comparison of the line of-sight displacement rates (Section 5.1). Slight differences in the locations of the subsidence patterns are also discussed below. While near-urban areas (marked by shaded polygon) have a high density of PS points, such points are sparse in the basin due to dense vegetation coverage and higher deformation rates.

The spatial distribution of the subsidence rates for the period of 2014–2017, obtained after decomposition of ascending and descending LOS velocity maps and corresponding to the component, is shown in Figure 5. Figure 5c shows that the city of Bursa, which is settled on the northwest slope of the Uludag Mountain, is not affected by the subsidence. The maximum subsidence rate at the center of the depression is ~25 mm/yr. The shape and the location of the subsidence pattern coincide with the extension of the agricultural zone, which relies heavily on groundwater exploitation.

5. Discussion

5.1. Self-Consistency Checking between the Ascending and Descending InSAR Observations

In order to show the self-consistency between InSAR measurements from two different viewing geometries, we first compared InSAR mean velocity measurements calculated from the two different Sentinel 1A/B datasets for the subsidence area in the Bursa plain and its surrounding region. Figure 6 shows the correlation between both datasets within the region of interest. The inter-comparison of the two line-of-sight displacement rate maps in this area was performed using a two-step procedure. In the first step, the mean velocity fields were resampled to a grid of 100 × 100 m to minimize the error of geo-localization (spatial mismatch between measurements). In the second step, we selected all pixels in the localized deforming areas that are common to both tracks. Then, we compared the distribution of the subsidence rates for both tracks. We found a root-mean-square value of 1.53 mm/yr and the coefficient of determination (R²) was close to 0.85. These values indicate that the two subsidence patterns (mean velocity) obtained from the different datasets are correlated with each other dominated by vertical motion. This agreement demonstrates that the deformation obtained from both datasets exhibits similar linear motions between 2014 and 2017, lending confidence to the subsidence patterns and amplitude determined by these datasets.

Comparison of the deformation rates between ascending and descending InSAR measurements of the same or different sensors is a common practice for quantitatively checking the consistency of measurements. The observed differences could stem from insufficient acquisitions, georeferencing uncertainty, seasonal biases, long temporal gaps, different temporal coverage, and, most importantly, different looking geometries. In order to show the impact of the looking geometry on the observed differences of the land subsidence induced by a single point source, we simulated a simple subsidence depression (Figure 7a). When this synthetic subsidence signal is observed with an incidence angle of 37.5°, considered as a mean incidence angle for Sentinel 1, the observed LOS deformation pattern will differ depending on the looking geometry of the sensor, either with a heading angle of −13.5° on ascending orbits, or −166.5° for descending orbits (Figure 7b,c). The maximum subsidence value moves away from the center along the azimuth look directions between the ascending and descending tracks. This idealized subsidence pattern gives a large dispersion and low correlation between the ascending and descending tracks as observed for our datasets.

5.2. Water Table Variations and InSAR Time Series

Using piezometric measurements, we can better constrain the hydro-geological processes that produce the observed subsidence. Agricultural activity that relies significantly on irrigation contributes to the groundwater depletion in aquifer systems of fine-grained sediments [71]. Ground-truth monitoring, satellites [72], InSAR [73], GPS [74], and numerical modeling have been used to understand the spatiotemporal variations of groundwater depletion. Many studies have shown that InSAR-derived subsidence patterns agree with water budgets [71,75,76]. Exploitation of aquifers decreases the groundwater volume, causing land subsidence that can be detected by InSAR. Given the lithological and hydrogeological data throughout the subsiding Bursa basin, we can constrain the amount of the groundwater depletion from the inversion of InSAR-derived land subsidence.

The Bursa ground water basin has a surface area of 208 km². To investigate the correlation between groundwater exploitation and surface deformation, we compared the temporal variations of our measured LOS displacements with temporal variations of the water table measured in two wells during the period of InSAR measurements (2014–2017). We only used LOS measurement points that fell within a 1 km radius around the wells. As seen in Figure 8, the decreasing trend observed on InSAR displacement time series and water level data show a very good agreement. The pressure head variations of both wells show ~10 and ~20 m of water level decrease starting in mid-2017, coinciding with an acceleration of the subsidence in InSAR data.

5.3. Compressibility of the Aquifer in the Bursa Plain

We used the subsidence measured by InSAR data and the level of water in the two wells (Figure 8) to estimate the compressibility of the aquifer that compacted due to the water level decrease. The compressibility of an aquifer is defined as:

$$\mathsf{\alpha }=-\frac{{\mathrm{dV}}_{\mathrm{T}}}{{\mathrm{d}\mathsf{\sigma }}_{\mathrm{e}}·{\mathrm{V}}_{\mathrm{T}}}$$

(2)

where ${\mathrm{V}}_{\mathrm{T}}$ is the total volume of the aquifer, ${\mathrm{dV}}_{\mathrm{T}}$ its temporal evolution, and ${\mathrm{d}\mathsf{\sigma }}_{\mathrm{e}}$ is the change in effective pressure in the aquifer [77]. We considered an increment of drawdown, ${\mathrm{dV}}_{\mathrm{w}}$, due to the pumping of water in the aquifer. This drawdown induces an increment of subsidence at the surface, ${\mathrm{dS}}_{\mathrm{u}}$, related to a compaction of the reservoir such that:

$$\frac{{\mathrm{dV}}_{\mathrm{T}}}{{\mathrm{V}}_{\mathrm{T}}}=-\frac{{\mathrm{dS}}_{\mathrm{u}}}{\mathrm{b}},$$

(3)

where $\mathrm{b}$ is the thickness of the aquifer, taken here as equal to 50 m. The change in effective pressure is:

$${\mathrm{d}\mathsf{\sigma }}_{\mathrm{e}}={\mathsf{\rho }}_{\mathrm{w}}{\mathrm{gdV}}_{\mathrm{w}},$$

(4)

where ${\mathsf{\rho }}_{\mathrm{w}}=1000$ kg·m⁻³ and g = 9.81 m·s⁻² are the density of water and the constant of gravity, respectively. The compressibility can be calculated as:

$$\mathsf{\alpha }=\frac{{\mathrm{dS}}_{\mathrm{u}}}{{\mathrm{b}\mathsf{\rho }\mathrm{gdV}}_{\mathrm{w}}}$$

(5)

Applied to the Bursa basin and considering a unit surface area of 1 m², Figure 8 shows that the drawdown increased in the range of $9<{\mathrm{dV}}_{\mathrm{w}}<18\mathrm{m}$ during the year of 2017 and the surface subsidence increased in the range of $0.035<{\mathrm{dS}}_{\mathrm{u}}<0.045\mathrm{m}$. As a result, the compressibility is in the range of $0.5·{10}^{-6}<\mathsf{\alpha }<2·{10}^{-6}$ Pa⁻¹, corresponding to typical values for clay and sand sediments [77,78]. These values agree with the alluvium sedimentary filling of the Bursa basin (see geological map in Figure 2).

5.4. Lithological and Tectonic Control over Land Subsidence

Alluvial deposits in the Bursa plain area in which we observed subsidence are made of Neogene and Quaternary deposits [51]. The thickness of these deposits varies from 50 to 300 m [55]. They are composed of predominantly soft clay and loose silt and sand [50]. Such highly compressible lithological units are more vulnerable to land subsidence due to groundwater extraction. We estimated a maximum vertical subsidence rate of 25 mm/yr for the 2014–2017 period in the Bursa plain. The spatial distribution of subsidence is characterized by an east-west oriented, elongated, elliptical shape of deformation along the axis of the Bursa plain, in areas where agricultural activities are well developed (Figure 5c). The vertical mean velocity field shows the subsidence to be on the order of 10 mm/yr in the part of the alluvial plain that extends for approximately 200 km² in the northern part of the Bursa city center. Within this subsidence zone, an oblong subsidence lobe (red in Figure 5a) is subsiding the fastest, sinking at the fastest observed rate of 25 mm/yr. The overall subsidence pattern in the study area shows very clear correlation with the shape of the basin, agricultural land use, and geological units of the study area. These correlations suggest that the land deformation is partially controlled by lithology and human activities. The fact that the pattern of subsidence mimics the shape of the basin may suggest that the alluvial thickness of the basin impacts the subsidence rate.

Similar to our results, Kutoglu et al. [1] have detected a subsidence signal in the Bursa plain between 2007 and 2010 and interpreted this signal as due to movements of the active faults that in the long term, have created the Bursa basin. However, the subsidence rate of 10.9 cm/yr proposed by Kutoglu et al. [1] is too high to be attributed to tectonic movements, such as shallow (surface creep) or lower crust deep creep [58]. As shown by GPS velocities in Figure 1b, the N-S horizontal velocity difference between the northern and southern shoulders of the Bursa basin that could be caused by aseismic slip on the basin boundary faults is negligible or too small to be detected in a few decades. As seen in Figure 9, there are no differences in the InSAR velocity across both the southern and northern boundary faults, confirming that there is no shallow creep on them. A deep creep would cause subsidence in the entire basin, not only in the alluvial plain. Therefore, our new observations and the strong correlation between the subsidence rate and ground water level suggest that most of the current subsidence detected in the past years is most likely due to the decline of the ground water table, rather than tectonic movements. An extremely high subsidence rate for the period of 2007–2010 was probably due to either droughts or more intensive pumping of the aquifer.

We have detected the very localized subsidence signal in the eastern tip of the basin overlapping with the agricultural area results from water pumping activities, but we do not have sufficient data to characterize this local signal.

6. Conclusions and Perspectives

The population of Bursa is increasing due to the city’s growing economy. This growing population is causing rapid urbanization, heavy land use, and land cover change. Such a rapid land development is forcing the land use to change from agricultural to more urban activities on the city’s outskirt and hinterland. We have detected ground subsidence in the northern part of Bursa at rates ranging from 5 to ~25 mm/yr in an area of 200 km² from the time series analysis of Sentinel 1-A/B images collected between 2014 and 2017. The use of both ascending and descending tracks, showing consistent results, is critical here to better quantification of the subsidence rate and further intercomparison with other SAR sensors’ data will continue to improve the subsidence monitoring.

In order to investigate the possible correlation between detected InSAR-derived land subsidence and water level change, we analyzed piezometric records for the same 2014–2017 period. These records show a water table drop starting in mid 2017 that correlates with a subsidence acceleration in the InSAR-derived time series. The coincidence of the InSAR subsidence spatial extent with the alluvial deposits suggests a hydrogeological control of the Bursa plain on the observed subsidence. The eastern part of the subsiding area is sinking at a maximum rate of ~25 mm/yr, where agricultural activities rely significantly on groundwater extraction. The spatial pattern of subsidence reveals that the urban environment is not yet affected by significant subsidence. Finally, no velocity gradient associated with active faults is presently observed.

All these observations indicate that the subsidence in Bursa plain is due to excessive groundwater extraction, not tectonic motions. This study demonstrates the capability of InSAR time series, combined with geological and hydrogeological data, to detect and monitor subsidence and identify potential control factors in the Bursa Province. Further work, complementary to InSAR, would be required to better monitor and asses the hazard related to land subsidence, using continuous GPS time-series in combination with in-situ measurements, including piezometric and geophysical data, such as borehole extensometers. In order to predict the future development of the land subsidence along the basin and differentiate the tectonic and anthropogenic contribution to the subsidence signal, further research should be performed to investigate elastic and hydraulic characteristics of the aquifer system. Because the entire Bursa basin is affected by land subsidence, and elaborate and orderly use of land is necessary to secure the infrastructures of the urban areas.