Localized Subsidence Zones in Gävle City Detected by Sentinel-1 PSI and Leveling Data

Among different sets of constraints and hazards that have to be considered in the management of cities and land use, land surface subsidence is one of the important issues that can lead to many problems, and its economic consequences cannot be ignored. In this study, the ground surface deformation of Gävle city in Sweden is investigated using the Persistent Scatterer Interferometry (PSI) technique as well as analyzing the historical leveling data. The PSI technique is used to map the location of hazard zones and their ongoing subsidence rate. Two ascending and descending Sentinel-1 datasets, collected between January 2015 and May 2020, covering the Gävle city, were processed and analyzed. In addition, a long record of a leveling dataset, covering the period from 1974 to 2019, was used to detect the rate of subsidence in some locations which were not reported before. Our PSI analysis reveals that the center of Gävle is relatively stable with minor deformation ranged between −2 ± 0.5 mm/yr to +2 ± 0.5 mm/yr in vertical and east–west components. However, the land surface toward the northeast of the city is relatively subsiding with a higher annual rate of up to −6 ± 0.46 mm/yr. The comparison at sparse locations shows a close agreement between the subsidence rates obtained from precise leveling and PSI results. The regional quaternary deposits map was overlaid with PSI results and it shows the subsidence areas are mostly located in zones where the subsurface layer is marked by artificial fill materials. The knowledge of the spatio-temporal extents of land surface subsidence for undergoing urban areas can help to develop and establish models to mitigate hazards associated with such land settlement.


Introduction
The management of cities and land use are challenging activities, where different sets of constraints, hazards, and preventive decisions have to be considered. Among the hazards, land surface subsidence and its subsequent effects on building and infrastructure can be considered as a major geo-hazard in many cities around the world. The principal causes of such deformation can be attributed to an aquifer-system compaction associated with groundwater/oil depletion, drainage of organic soils, underground mining, hydrocompaction, natural compaction, sinkholes, stress ended by a load of construction and foundation type, and thawing permafrost [1][2][3][4].
Globally, numerous cities and regions undergoing development are characterized by different rates of land movement phenomena, such as in Hong Kong [5], London [6], Las Vegas [7], and Tehran [3]. Therefore, knowledge of large-scale land surface subsidence could be an important factor in the long-term planning and managing of cities.
for different purposes, such as urban spatial planning and development activities, environmental management, geo-hazard monitoring, and assessment efforts. For instance, the PanGeo project (https://www.eurogeosurveys.org/projects/pangeo/), which was a geo-hazards European project, used satellite radar data integrated with local geological information to provide geo-hazard information for the largest 52 cities across Europe, including Stockholm and Gothenburg in Sweden. Therefore, the integration of the radar satellite data with the local in situ data (e.g., geology and precise leveling data) would be of great importance for the geo-hazard mapping.
Generally, land subsidence impacts in urban regions are quite numerous and can be classified into infrastructural, environmental, economic, and social impacts. In Gävle, which is one of the old coastal growing cities in Sweden, the subsidence impacts can be seen in forms of cracking and damage of buildings and infrastructure, especially in the northern part of the city, which negatively will affect the rehabilitation and maintenance cost and the quality of the living environment. For instance, according to the latest surveys and reports made by the Gävle municipality (www.gavle.se), some streets in the northern part of the city (i.e., Norra Strandgatan) were repaired due to the subsidence problems. Furthermore, visible related problems probably due to subsidence of the buildings in the Gävle Strand area, have been reported by, e.g., WSP and SWECO Companies. Additionally, subsidence for 13 benchmark points in the area of Alderhomen, around Gävle beach, ranging between 4 and 12 mm, have been reported by Lindén and Rosenqvist [41], using leveling data collected between 2006 and 2011. Hedberg and Ottekrall [42], have reported similar problem for three different buildings located in the city center using leveling data between 1985 and 2011.
The city of Gävle crosses the Gävleån River that formed a large delta at the outlet in the Gulf of Gävle. The sub-surface geology map of the city shows areas marked with artificial fill that have been used for urban development activities. According to the Department of the Environment, Transport, and the Regions, and Environment Agency in the UK (DETR/EA, 1998) "the fill material at any site is commonly an admixture of organic, chemical and inert material which can lead to serious problems in the built environment". Thus, the ongoing and expected development of the city according to the municipality of Gävle, along the coastal line, within the formed delta and the artificial fill area, is prone to subsidence hazard and needs further investigation with, e.g., satellite remote sensing (InSAR) for deformation mapping. According to the Gävle municipality, there are plans for numerous large construction projects; e.g., Näringen, which is one of the largest urban transformations in the municipality and targeted to become one of Europe's most sustainable districts with 6000 new houses and space of 450,000 m 2 of gross area businesses. The targeted area is located around and within the identified geo-hazard zone areas according to the previous studies. Moreover, the availability of long precise leveling observation records that last for about 45 years in monitoring the stability of some locations in the city, in addition to other information from remote sensing data would help and motivate us to further study and map the localized deformation in the city of Gävle.
The aim of this study is to generate a precise deformation map of Gävle city, to detect and highlight the hazard zones, which relatively deform with a higher rate in this city, and correlate the results with subsurface geology. In this study, Sentinel-1 satellite radar datasets provided by the European Space Agency (ESA) have been processed and analyzed using the PSI technique, which is essentially an improved version of the standard PS-InSAR algorithm, developed by Perissin and Wang [43] and implemented in the SARPROZ software. For this purpose, two freely available datasets of Sentinel-1 radar images in ascending and descending geometries were used, that covered more than a five-year period starting from January 2015 to May 2020. Available long records of precise leveling data, back to early 70s, observed and processed by the municipality and the University of Gävle, were utilized to examine the obtained PSI results for the Sentinel-1 observation period as well as to analyze the long-term stability of some locations in the city which were not reported before.

Study Area
The study area covers about 24 km 2 of Gävle city which is located by the Baltic Sea in east-central Sweden, at about 60 degrees north and 17 degrees east (see Figure 1). Gävle city, which is the capital of Gävleborg County is the oldest in the historical Norrland, with a population of about 102,000 in 2019, which makes it the 13th largest city in Sweden. According to the quaternary deposit map of the study area (Figure 2), provided by the Geological Survey of Sweden (SGU), the surface geology of the city mainly consists of post-glacial sand which is dominant and covers most of the city center and extends toward the west. A glaciofluvial sedimental strip crosses the city from southwest to the northeast. Glacial clay and Post-glacial clay layers located toward the south, the east, and partially to the north of the city. In addition, some areas to the northeast are of artificial fill material. The two banks of the outlet part of Gävleån River (i.e., about 6 km across the city) that crosses the city from the west to the east toward the Baltic Sea, consist of artificial fill material and have been constructed. The fill areas were part of the original route of the river canal and extend toward the north from the riverbank. Furthermore, the outer areas of the city are almost sandy till layers type. According to the SGU (https://www.sgu.se/) the estimated soil depth in the city ranges between 3 to more than 50 m in some areas along the coastline, and between 10 and 20 m in the artificial fill areas. Figure 1 shows the location map of the city, urban zone, and the area covered by leveling measurements (the red square in Figure 1b).

PSI Method
According to Crosetto et al. [19], the comprehensive DInSAR phase components can be read as below, if a digital elevation model (DEM) of the imaged scene is available (i.e., the topographic term can be simulated and subtracted from the interferometric phase ∆ϕ int yielding the so-called DInSAR phase components ∆ϕ D−int ): where ∆ϕ Disp is the displacement phase term in the LOS direction, ∆ϕ Top−res represents the residual topographic error components, ∆ϕ atm is the atmospheric phase components at an acquisition time of each image (Master and Slave), and it can be reduced by applying stacking technique [44], ∆ϕ orb is the orbital error phase components of each image (Master and Slave), and it can be estimated using precise orbits [45], ∆ϕ noise is the phase noise which can be minimized with filtering and multi-look techniques [45,46]. The last term of Equation (1) represent the general 2π ambiguity (integer value n) is a result of errors in the phase unwrapping process, which is necessary to solve for this ambiguity, i.e., the DInSAR phases are bounded in the range (−π, π) [11,47]. The goal of any DInSAR technique is to estimate the displacement phase components ∆ϕ Disp from the ∆ϕ D−int , which requires its separation from the other components in Equation (1).
To overcome the limitations of the deformation monitoring using the DInSAR (i.e., temporal and geometrical decorrelation and atmospheric disturbances), Ferretti et al. [13,14] developed the first PSI technique based only on a subset of pixel analysis. It is called permanent scatterers (PS) and assumed to possess stable phase values during the time series of the whole used datasets of the same area. From these datasets, a single master is selected considering a high coherence in all formed interferograms over the whole datasets. Temporal sampling of the data should be regular to some extent. Prior to the PS estimation, a DEM of the area should be used to remove the topographic phase components [19]. Accuracy of about 20 m is sufficient for the used DEM [14].
The PS is the strong reflecting objects that are dominant in a cell and it can be identified based on the amplitude stability index and the temporal coherence [48]. Stable pixel selection depends on the use of amplitude dispersion index [49], which can be computed as in Ferretti et al. [50]: where σ A and m A are the standard deviations and the mean of the amplitude, respectively. Therefore, a pixel with high reflection has a low D A . Temporal coherence explains the mean difference between the observed and modeled phase for each pixel and interferogram [9]. Over the selected PS points, the atmospheric phase components should be estimated for every interferogram, where the interpolated map is called the atmospheric phase screen (APS). According to Colesanti et al. [51], at least 25 images are needed in addition to a PS density of about 5 to 10 PS/km 2 to achieve proper estimation and removal of the APS effects. Furthermore, a deformation model has to be introduced prior to the velocity estimation e.g., linear, nonlinear, and seasonal models [52]. It is worth mentioning here that, the PSI technique [13,14] was followed by several contributions and new approaches such as PSI procedure for crustal displacement in non-urban environments [23], the adaptation of the GPS LAMBDA method to PSI [53], the SqueeSAR algorithm, which is an extension of the original PS-InSAR algorithm [54], and the PSI approach introduced by Perissin and Wang [43], which is called the QPS algorithm. The PSI method by Perissin and Wang (i.e., QPS) exploits partially coherent targets, and it has been incorporated in the SARPROZ InSAR processor software which is used in this study.

Sentinel-1 Data and Analysis
In this study, freely available C-band Sentinel-1 A and B data from the European Space Agency (ESA) were downloaded (https://scihub.copernicus.eu/dhus/#/home) and used to detect and quantify the ground surface deformation of Gävle city, using the PSI method [18,51]. Two datasets consisting of 41 and 50 Single Look Complex (SLC) images on ascending and descending geometries were processed and analyzed over the study area covering periods from January 2015 to May 2020 and from June 2015 to May 2020, respectively (see Table 1). The temporal resolution for the acquired images was about one month, with some exceptions due to the unavailability of some images with certain characteristics or some images were removed due to the winter season which may reduce image coherence. For the data processing, SARPROZ software [49] was used, which could determine the LOS rate and its associated uncertainty for each individual PS point. The SLC images were co-registered to a single master. Then, the 50 m grid local DEM, provided by the National Land Survey of Sweden (Lantmäteriet), and the precise orbits for each image were used to remove the topographic and flat earth components, respectively. Reference points were selected among the selected Persistent Scatter Candidates (PSCs), which are relatively unaffected by deformation. For the PS point selections, a quality factor called amplitude stability index is used and defined as 1 − D A in the software.
The first and second sets of data i.e., the 41 ascending and 50 descending images, were used to analyze the land surface of the study area. We used the images dated by 2016-11-06 (for ascending) and 2017-11-13 (for descending) as master images, respectively, according to certain characteristics including the length of the baseline, where a high stack coherence is achieved.
The APS effect was estimated and removed in the location of the PSCs, which were selected based on amplitude stability indexes of 0.7 and 0.8 as quality parameters for the ascending and descending datasets, respectively. Temporal coherence of 0.60 was used to select all persistent scatterers (PS) points. For each dataset, the line of sight (LOS) velocity and displacement time series were estimated utilizing a linear model. To ensure reliable results, the two reference points (located in Sätra, Figures 3 and 4) with zero displacements, were selected in the same very stable area and away from the known deformed zones. Furthermore, the estimated LOS displacement velocities were decomposed to the east-west and vertical components by combining ascending and descending results using the following Equations (Equations (3) and (4)) presented by Milillo et al. [55] and implemented and used in SARPROZ software.
where θ and d are the incidence angle and the estimated displacement velocity, respectively, in the ascending (Asc) and descending (Des) geometry. AD pair (i.e., Ascending Descending pair) method was implemented by using 10 m as a maximum offset distance (in planar and height) between the selected pair.

Levelling Data and Analysis
A long record of leveling measurements has been conducted and processed by the municipality of the city of Gävle and partially by the University of Gävle, covering the period from 1974 to 2019. The data were used to examine the obtained PSI results as well as to analyze the ground stability of some locations in the study area. The record consists of four short separated leveling lines ranging between 100 and 500 m and has been established by the municipality to monitor four different buildings located in the center of the city (see Figure 1c) within the vicinity of the Gävleån River. Each line consists of several metal pegs that mounted in the building foundation and associated with a certain benchmark (BM). A double run leveling procedure has been performed using different types of leveling instruments including Leica DNA03, NA2, and NA720. The measurements are referenced to the official national height system RH 2000 in Sweden. Periods of observation records for each building varies and ranges between 6 years and 34 years. For instance, the levelling records for Building-1 ranged from 1985 to 2019, Building-2 ranged between 2000 and 2019, Building-3 ranged between 1976 and 1982, and Building-4 between 1974 and 1988 (Figure 1c). In addition, the temporal resolution for repeated leveling observations also varies from months to several years (this will be discussed in Figure 6). By using regression analysis, the annual rates of vertical changes have been computed at each peg.

Geological Data and Analysis
According to the National Geotechnical Institute "Statens Geotekniska Institute" (SGI) assessment report of Gävle municipality [56], the dominant soil along the coastal stretch of the Gävle municipality is almost blocky and large-block moraine, while areas with sand and fine sediments, such as clay and silt, are found in areas near to watercourses and small lakes as well as in sea bays (cf. Figure 2). Therefore, based on the soil condition, the report highlighted five areas in the region, attributed to the hazard of landslide, race, and erosion [56]. The identified five areas almost consist of silt and/or clay as well as fill that probably underpinned of fine sediment and adjacent to the coastal lakes' watercourses. The fill term in this study is used to describe ground that has been formed by human activities in which natural soils or sometimes man-made materials are deposited, in contrast to natural soils which has its origin in geological processes [57]. At the outlet of the Gulf of Gävle, the Gävleån River, which crosses the city, formed a large delta, which has been covered using fill material accordingly and used for urban development. In addition, the associated small river arms are also channeled or have been filled.
By combining the geological information of Gävle city with the resulted PSI deformation map, it is possible to assess the hazard of land movement in the city. The combined ascending and descending PSI results (i.e., vertical components) and their rate of movement were analyzed based on the provided 1:25000-1:100000 quaternary deposit map by the SGU (Figure 2). The PS points were overlaid on the geological map where the rate of displacement can be evaluated accordingly (cf. Figure 7).

Results and Discussions
In this section, we report our results of the estimated land subsidence rate over the city of Gävle and the relation between the obtained PSI results with the available leveling data and the geological information.

PSI Results
The two Sentinel-1 datasets which consist of ascending and descending images were analyzed over the city, where the PSI technique was applied to detect relative ground deformation, using the processing criteria in Section 3.2. By using 0.6 as a temporal coherence mask, a total of 7265 and 9314 PS points were obtained using ascending and descending geometries, respectively, by employing SARPROZ software.
Generally, the obtained ascending and descending results show good agreement in the LOS direction relative to selected references, which confirm the stability of the two references that have been selected in the same location (see Figures 3 and 4). It is worth noting that, the selection of the reference points has to be very careful due to its effects on all PS points. The selected two reference points (ID 5776 and ID 1976) have zero velocity and 0.94 and 0.99 temporal coherence, respectively, and are located in a very stable part of the city (Sätra area) that consist of types of sandy till layers. Figure 3 shows all generated ascending PS points, the main deformed areas (Area-1 in the white square), and five different small areas (i.e., A1-A5) where the PS rates have been averaged to generate time series for the selected five areas using 5 to 10 PS points (see Figure 5a), in addition to the location of the reference point. Point selection for averaging is based on the closeness of the points to each other (i.e., within 100 m radius) and having the minimum coherence of 0.75. The identified hazard zone Area-1, which is located toward the northeast of the city center, shows the maximum negative LOS movements in the city, based on the estimated PS points in the area which reaches up to −6 mm/yr. Some sparse areas and points showed a notable subsidence rate ranging between −1 mm/yr and −3 mm/yr, e.g., in the railway lines (area A5) to the northwest of Area-1 (Figure 3). The rest of the land surface in the city is almost relatively stable, where the maximum movement is ranged between about −2 mm/yr and +2 mm/yr. The small and sparse positive LOS movements may be attributed to ground heave (i.e., when the soil is saturated by water or accumulated due to construction activities).  Figure 4 shows all the generated descending PS points and the area with maximum subsidence rate (Area-1), and similarly the five different small areas (i.e., A1-A5), in addition to the location of the reference point. In both cases (i.e., ascending and descending results), the density of the obtained PS results is relatively poor (about one PS for each 2500 m 2 ), based on the medium spatial resolution of the SLC Sentinel-1 radar data which is 5 m × 25 m. In contrast, more dense PS points would be expected when using high-resolution images e.g., X-band images, which offer a resolution up to 1.1 m × 0.6 m as for TerraSAR-X [58]. Figure 5 shows the time series and the accumulative LOS relative movement in millimeter for the five ascending and descending small areas (i.e., A1-A5), and the reference points ID 5776 and ID 1976, respectively. There are a lack of data during periods of 31 October 2015 to 15 June 2016 (about 7 months) and between 12 January 2018 to 18 April 2018 for the descending dataset, in addition to different periods in the ascending set, ranged between 2 and 5 months (Figure 5), despite the short revisit time of Sentinel-1 sensors which is 6-12 days. The results show a continuous subsiding rate with a notable deviation in e.g., January, and February which can be attributed to the winter season. The reference points (ID 5776 and ID 1976) demonstrate a very stable and consistent time series with a minor deviation in some months e.g., February.

Leveling Data And the PSI Results
By utilizing the leveling records associated with the four buildings (Figure 1c), the annual rate of change at each leveling peg has been computed using the regression analysis method. Figure 6 illustrates the time series and the computed rate of change at a selected peg for each building. The selection of the peg point is based on its closeness to the generated PS points on or very near to the building to facilitate the validation process. It is worth noting that, the leveling records of the Buildings 1 and 2 (Figure 6a,b) extend for a long period (1985 to 2019) and overlap with the satellite radar data (2015 to 2020). In contrast, the leveling records of the Buildings 3 and 4 (Figure 6c,d) ended at the years 1982 and 1988, respectively, i.e., there was no overlap between the leveling records and the SAR data. However, even that old leveling records for Buildings 3 and 4 in that period was used here to infer if the present-day subsidence rate from PSI is comparable with the one obtained from old leveling records at different time period. For validation purposes, the available leveling information data were used to confirm the obtained PSI results, where the generated ascending and descending PS points on or close to the buildings and the leveling pegs were considered.
The LOS movement rates of all generated PS points from the PSI analysis are relative to one stable reference point, ID 5776 for ascending and ID 1976 for descending, while the computed displacement rates of the metal pegs, mounted in the buildings are referenced to specific benchmarks (BM), and these benchmarks could be moving with respect to the stable reference of the PSI analysis. To overcome this relative problem, the trend rate in the nearest PS points to each peg in the building will be computed relative to the nearest PS point to the leveling benchmark. Then, a comparison between the computed leveling rate (i.e., between the BM and the peg) and the relative rate between the PS points that represent the BM and peg (i.e., between the two PS points near to the BM and to the peg) can be performed. It is worth noting also, that the PSI results are in the LOS direction, while the leveling rate results in the vertical direction. To compensate for such differences, the PSI results (i.e., the displacement rate and the accumulative displacement) have been converted to vertical and horizontal (only possible for east-west) movements from the LOS direction by combining the ascending and descending datasets according to Milillo et al. [54] (see Figures 7 and 8 and Equations (3) and (4)). Table 2 shows the track type i.e., ascending and descending pair (AD), point number (ID), the relative vertical displacement rate in mm/yr, the vertical accumulative displacement in mm, and the coherence of the generated combined PS points on or near to the four buildings using about five years' radar data. In addition, the corresponding computed annual rate of the leveling measurements (Pre. Lev) in mm/yr for the four buildings, using different records ranged between 1974 and 2019, have been presented.   Table 2. Comparison between the computed precise leveling rates (Pre. Lev) using four different leveling records relative to four different BMs, and the relative vertical rate of the generated combined PS points in the four validation buildings using the ascending and descending combination of the 41 and 50 SAR images, respectively, collected between 2015 and 2020. According to Table 2, the relative vertical PS results of point (ID 498) on Building-1 shows close agreement with the computed leveling rate at Point 7 (nearest combined PS point to the leveling point) i.e., −0.9 mm/yr vs. −1.2 mm/yr, respectively. In Building-2, the relative vertical results of point (ID 430) shows very close agreement in terms of the rate and the accumulative displacement with the leveling data at Point 6. In Building-3, despite the relativity short and old record of the leveling data (from 1976 to 1982) when compared to radar data (from 2015 to 2020), a lesser agreement is revealed in terms of the displacement rate between the leveling records at point 6 and the relative vertical combined PSI result at point (ID 396), which is −2.0 mm/yr and −0.6 mm/yr, respectively. In Building-4, the vertical relative combined result shows almost complete stability (i.e., zero rates) with a very minor cumulative displacement for the 5 years (−1.6 mm), while the historical leveling record (from 1974 to 1988) shows relative displacement rate of −1.8 mm/yr and accumulative displacement of −31.0 mm for a total of 14 years (i.e., there is no agreement). The result reveals that Building-4 was experiencing a high relative subsidence rate during the years after construction phase, due probably to the load and foundation type, and reached complete stability during the last years.
The above comparison between the obtained relative PSI results in the vertical direction, using combined ascending and descending data for about 5 years, and the precise leveling information record for different periods at the four buildings reveals close agreement in Buildings 1 and 2.

Subsidence and Subsurface Geology Effect
By correlating the geological information of Gävle city with the resulted combined (ascending and descending) PSI deformation map of the land surface, it is possible to assess the vertical land movement in the city. The combined measured PS points and their rate of the vertical movement were analyzed based on the provided 1:25000-1:100000 quaternary deposit map by the SGU (Figure 2). The combined PS points were overlaid on the geological map where the high rate of displacements can be illustrated accordingly ( Figure 7). Moreover, the resultant horizontal movement in the east-west direction has been overlaid on the topographic map of the city (Figure 8).
Based on our results, the Area-1 (see Figures 3, 4 and 7) shows the PS points with maximum displacement rates that reach up to −6 mm/yr. It is clearly visible from the corresponding geological map that the subsiding area belongs to the artificial fill zone that covers the banks of the outlet of Gävleån River and extends toward the north [59]. Some sparse PS points with a relatively high rate of subsidence are located on areas marked with clay to silt deposits (pale yellow solid zone in Figure 7) in the southeast of the city. The land surface of the central part of the city is relatively stable and it is almost located on the post-glacial sand deposits. However, some sparse PS points, mostly along and around the river, show a westward movement about 1-2 mm/yr relative to the selected reference point (Figure 8). This can be further investigated in future by analyzing the higher resolution radar satellite images, e.g., TerraSAR-X, which provides deformation maps with denser PS points.
Field investigation reveals that most of the construction buildings in areas A1 and A4, those located in the artificial fill Area-1, are of light type (e.g., warehouses and workshops). Meanwhile, area A2, which is the newest build area in the city, contains multi-story buildings. Recently, visible related problems probably due to subsidence of the buildings in area A2 has been reported to construction consultancy companies. Additionally, it is worth noting that, the fill area to the northeast of Area-1, which contains a large number of huge storage tanks, shows very stable PS results, with minor vertical movement of order of −1.0 to −2.9 mm/yr concentrated in the railway line in the vicinity of the oil storage tanks area, while the storage tanks area shows vertical movements of order less than −1.0 mm/yr. The storage tank foundations are probably built with special care to the subsurface geology type in the area. Furthermore, uneven subsidence evidence has been reported in the church building in the city center (Heliga Trefaldighetskyrkan) to the east of Building-1, where the very visible cracks can be attributed to different reasons such as building foundation settlement or age-related settlement or the close location to Gävleån River. However, due to the resolution of the PSI results, it was not possible to get enough PS points on this building that could be used for evaluation.

Conclusions
In this paper, we have carried out a comprehensive study to generate an accurate deformation map of Gävle city, to detect and highlight the hazard zones. For this purpose, two datasets of Sentinel-1 radar images in ascending and descending geometries, provided by the ESA, and covered the period of January 2015 to May 2020, were processed and analyzed using the PSI method and combined ascending descending datasets. Available long records of precise leveling data were utilized to validate the PSI results. Furthermore, a correlation analysis was performed between the PSI results and the available geological maps for the study area.
Our results show that land surface in Gävle city is relatively stable with exceptions of some localized subsiding zones. There is no obvious pattern of deformation observed for the whole or a big part of the city. Only localized deformation zones were detected by our analysis that have been subsiding in the last five years. The correlation of our results with subsurface geology types, i.e., the artificial fill, suggests that the subsidence occurs due to relatively weak subsurface layers which either are affected by loading of new constructions or due to hydrocompaction. A further analysis is suggested to analyze the hydrogeological and geotechnical measurements and correlate them with InSAR results to better understand the ongoing land deformation process in this city.
The detected subsiding zones show a maximum displacement rate that reaches up to −6 ± 0.46 mm/yr in the LOS direction. To the northwest of Area-1 (i.e., in Gävle railway station's marshalling yard), the area A5 is experiencing average subsidence ranging from −1 ± 0.47 mm/yr to −3 ± 0.53 mm/yr. The rest of the city, including the city center, is almost relatively stable with minor deformation rates ranged between −2 ± 0.5 mm/yr to +2 ± 0.5 mm/yr in vertical and east-west components. A comparison between the computed precise leveling rates using the available four different leveling records relative to four different BMs, and the relative vertical rate of the generated PS points in the four validation buildings, shows close agreement in Buildings 1 and 2. According to the correlation analysis between the obtained PSI results and the geological information of the city, the achieved result is highly correlated with the quaternary deposit map, e.g., in the detected hazard zone that is reported as an artificial fill area (see Figure 7). This conclusion can be considered by the Gävle municipality for their future projects, e.g., in Gävle harbor.
Knowledge of the spatial and temporal extents of land surface subsidence for cities that are vulnerable and undergoing development can help to develop and establish measures to mitigate hazards associated with such land settlement.
Author Contributions: N.A.A.G. wrote the manuscript and processed and analyzed the data. M.B. designed the study and helped to write and prepare the manuscript. He also contributed to historical and precise leveling data processing. F.N. contributed to supervision and writing the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding:
We are immensely grateful to J. Gust. Richert Stiftelse foundation (Project number 2019-00485) for their financial support for this study.