Countrywide Monitoring of Ground Deformation Using InSAR Time Series: A Case Study from Qatar

: Over the past few decades the country of Qatar has been one of the fastest growing econ-omies in the Middle East; it has witnessed a rapid increase in its population, growth of its urban centers, and development of its natural resources. These anthropogenic activities compounded with natural forcings (e.g., climate change) will most likely introduce environmental effects that should be assessed. In this manuscript, we identify and assess one of these effects, namely, ground deformation over the entire country of Qatar. We use the Small Baseline Subset (SBAS) InSAR time series approach in conjunction with ALOS Palsar-1 (Jan. 2007 to Mar. 2011) and Sentinel-1 (Mar. 2017 to Dec. 2019) synthetic aperture radar (SAR) datasets to assess ground deformation and conduct spatial and temporal correlations between the observed deformation with relevant datasets to identify the controlling factors. The findings indicate: (1) the deformation products revealed areas of subsidence and uplift with high vertical velocities of up to 35 mm/yr; (2) the deformation rates were consistent with those extracted from the continuously operating reference GPS stations of Qatar; (3) many inland and coastal sabkhas (salt flats) showed evidence for uplift (up to 35 mm/yr) due to the continuous evaporation of the saline waters within the sabkhas and the deposition of the evaporites in the surficial and near-surficial sabkha sediments; (4) the increased precipitation during Sentinel-1 period compared to the ALOS Palsar-1 period led to a rise in groundwater levels and an increase in the areas occupied by surface water within the sabkhas, which in turn increased the rate of deposition of the evaporitic sediments; (5) high subsidence rates (up to 14 mm/yr) were detected over landfills and dumpsites, caused by mechanical compaction and biochemical processes; and (6) the deformation rates over areas surrounding known sinkhole locations were low (+/−2 mm/yr). We suggest that this study can pave the way to similar countrywide studies over the remaining Arabian Peninsula countries and to the development of a ground motion monitoring system for the entire Arabian Peninsula.


Introduction
Ground deformation has been reported from many areas around the world, and in many cases the observed land deformation was associated with losses of properties and life. Recent advances in Interferometric Synthetic Aperture Radar (InSAR) time series methodologies (e.g., Persistent Scatter Interferometry [PSI] [1,2], Small Baseline Subset [SBAS] [3,4] and the combination of persistent and distributed scatterers [DS] [5,6]) enabled the cost-effective detection and monitoring of ground deformation caused by anthropogenic and natural factors over large areas worldwide. Example applications include monitoring ground deformation in urban areas where large-scale infrastructure projects were conducted [7,8], subsidence in deltas [9,10], sinkhole-related deformation and collapse [11][12][13], and subsidence due to excessive groundwater extraction from fossil aquifers in semiarid and arid environments [14][15][16][17].
Recently, InSAR-based ground deformation monitoring systems have been set up and planned across an entire country (e.g., Norway [18], Germany [19], The Netherlands [20], Denmark [21], Italy [22,23]) or even a continent (e.g., Europe [24,25]). The main purpose of these systems has been to distribute InSAR-derived deformation products, but the interpretation of observed deformations was left largely to the users of these platforms. To the best of our knowledge none of these countrywide radar interferometric studies have been conducted over developing countries or arid areas.
In this study, we assess the ongoing and earlier ground deformation across the country of Qatar (11,571 km 2 ) and investigate its driving forces. Qatar has been subject to rapid economic development and population increase since 2005. The population was 744,029 in 2004, increased by some 300% in 2015 (2,404,776), and by 360% (2,715,919) in 2020 [26]. This rapid increase in population, coupled with the implementation of ambitious development plans and projects, was associated with progressive urbanization, land-use/landcover (LULC) change, and groundwater level variations related to intense groundwater extraction and injection operations. These anthropogenic activities compounded with natural forcings (e.g., climate change) will most likely introduce environmental effects that should be assessed.
We used Sentinel-1 data to investigate: (1) the recent and ongoing deformation (Mar. 2017 to Dec. 2019) over the entire country; (2) the progression of deformation through time by comparing the recent/ongoing deformation over the examined period to earlier deformation (Jan. 2007 to Mar. 2011) that was extracted from Advanced Land Observation Satellite ALOS Phased Array L-band Synthetic Aperture Radar (Palsar-1) data; and (3) the factors controlling the ongoing and earlier deformation by conducting spatial correlation of the observed deformation with static (geological and LULC maps, distribution of sinkholes and dumpsites) and temporal datasets (high resolution optical satellite imagery, precipitation data, groundwater levels, surface water area within the sabkhas, and digital elevation models [DEMs]).
Conducting a countrywide deformation monitoring exercise over Qatar is advantageous in many respects. Firstly, the project is conducted jointly with the Ministry of Municipality and Environment (MME) of Qatar, the agency responsible for collecting the detailed environmental datasets that were instrumental in deciphering the factors controlling the observed deformation. Secondly, Qatar is an arid area, with minimal vegetation cover, making it ideal for effective radar interferometric applications [3]. Thirdly, Qatar's fast growing economy and population over the past two decades have led to large infrastructure development projects, including developing a metro network in the Doha area, building manmade islands, and developing a new airport on reclaimed land. Monitoring the land deformation associated with natural and anthropogenic factors can assist in the identification and implementation of appropriate remedial processes in a timely fashion and could serve as a replicable model in similar settings worldwide.

Study Area
The State of Qatar is located on the eastern side of the Arabian shelf extending into the Gulf as a 160 km long and a 60-90 km wide peninsula along a broad, north-south trending fold, the Qatar arc [27] (Figure 1a). The relief of the Qatar Peninsula is, in general, low. Analysis of a lidar-based DEM reveals that the elevation ranges from −2.7 to 102 m in reference to the Qatar National Height Datum (QNHD), and averages 24.5 m. Many of the areas with elevations below sea level are currently occupied by inland sabkhas (e.g., Dukhan (Sabkha 1; area: 69.7 km 2 )) and coastal sabkhas (e.g., Mamlaha (Sabkha 2; area: 14.4 km 2 ), Afjat Wadi Abu Al Ghurban (Sabkha 3; area: 15.9 km 2 ), and Afjat Saadoon (Sabkha 4; area: 8.8 km 2 )) ( Figure 1). The elevation of the Dukhan sabkha ranges from −2.7 to 2 m, and for the coastal sabkhas it ranges from −1.3 to 1 m in reference to QNHD. The Dukhan Sabkha, the largest inland salt flat (70 km 2 ) in the Qatar Peninsula, is located within an upturned L-shaped depression, where highly saline groundwater is found at shallow depths ranging from 2 to 120 cm [28]. The sources of the groundwater within the sabkhas are thought to be discharge of continental groundwater, evaporative pumping of seawater, and limited runoff [29]. The shallow water table, high temperatures, high rate of evaporation, and saturation of the sabkha sediments by brine promote the formation of evaporites within the capillary zone and on the surface as well [30].
Approximately 80% of Qatar's land surface is covered by the chalky dolomite and limestone of the Eocene Dammam Formation (Figure 1a). In southern Qatar, the Dammam Formation is unconformably overlain by 40-80 m thick limestone, dolomite, and the evaporite rocks of the Miocene Dam Formation [27]. Quaternary deposits are largely found on the coastal and near-coastal areas and are composed of sabkhas, and unconsolidated sediments (fine-grained alluvium, eolian sand, beach deposits, etc.; Figure 1a).
Throughout the Pleistocene the climate alternated between wet and dry periods [17,31], and it is during these wet periods that the karst was formed [32]. The karst formations of central Qatar have been associated with the subsurface dissolution of carbonate and sulfate deposits during Pleistocene wet climatic periods, whereas the pitted karst terrain in northern Qatar was associated with differential dissolution along joint-flow drainage [32]. The discovery of many cavities during the implementation of massive infrastructure and construction projects in recent years was one of the main reasons behind launching a new national initiative for subsurface geological and geotechnical mapping [27,33]. Shallow groundwater systems are found in the karstic Dammam and the Upper Umm El Radhuma and Russ aquifers, and the majority of sinkholes were reported from the Upper Dammam formation (Figure 1b) [32].
Qatar has been injecting treated wastewater (28% of the total treated wastewater, or 63 million m 3 , in 2017) into aquifers to replenish its scarce groundwater resources [34]. Changing groundwater levels due to treated wastewater injection and also extractions (dewatering) for large construction projects raised concerns about the extension of already formed sinkholes and the creation of new ones within urban environments in Qatar. The increase in urbanization and construction activities also created an increase in solid waste production. Solid waste in Qatar is currently being stored in two landfills (Umm Al Afai (L1) and Mesaieed (L2)), two dumpsites (Rawdat Rashed (D1) and Umm Thanyatain (D2)), and a solid waste processing plant (P1; Figure 1c

Data and Methodology
The adopted approach involved three major steps. These were: (1) the extraction of deformation rates across Qatar throughout the ALOS Palsar-1 and Sentinel-1 periods (Step I), (2) the collection and analysis of relevant auxiliary climatic (precipitation), geodetic (GPS stations), and hydrologic datasets (temporal groundwater levels; Step II), and (3) the correlation of the extracted rates of deformation with all of the relevant auxiliary datasets in search of causal effects (Step III).

SAR Datasets and SBAS Time Series Analysis (Step I)
SAR data from two different instruments were used to investigate the magnitude and nature of ground deformation: ALOS Palsar-1 (L-band with a 23.62 cm wavelength) and Sentinel-1 (C-band with a 5.56 cm wavelength). ALOS Palsar The SAR instrument with its side looking geometry combines the projections of 3dimensional ground deformation into the line-of-sight (LOS) direction. It is favorable to combine ascending and descending datasets from the same SAR instrument to retrieve 3dimentional deformation [36]. If combining ascending and descending datasets of the same instrument is not possible and datasets from different SAR missions are being utilized together, it is preferred to use the data in the same viewing geometry to minimize the differences due to the look direction, especially where deformation has a significant horizontal component [37]. In our study area, the ALOS Palsar-1 archives have [16][17][18][19][20][21][22] scenes available in ascending mode and only one in descending mode. This forced us to use the ascending Palsar-1 dataset only. Ascending and descending Sentinel-1 data were available over the study area and were both analyzed. The ascending Sentinel-1 dataset produced a considerably higher number of scenes/interferograms with strong atmospheric artifacts, phase jumps and unwrapping errors which prevented us from acquiring reliable deformation estimates, whereas the descending dataset provided more reliable and spatially complete results. Given the above-mentioned constraints we chose to work with descending Sentinel-1 and ascending ALOS Plasar-1 datasets and to back-project the LOS deformation into the vertical direction to enhance the compatibility of the estimated deformation rates from two datasets. Back-projection was performed by dividing the estimated LOS deformation by the cosine of the incidence angle [38]. In doing so, we assumed that the deformation is dominated by the vertical direction, a reasonable assumption given the flat topography of the study area, and the lack of any reported lateral tectonic deformational features.  Figure 2b-d, are images included in, and those excluded from, the final SBAS process (green and red dots, respectively), the main primary image (yellow dot), and the interferometric pairs (two dots joined by a line). * SBAS networks are plotted in Figure 2.
A multitemporal differential SAR interferometry approach was adopted in this study to enable the measurement of surface deformation velocity and its evolution through time. The technique is based on the isolation of deformation phase component between the two SAR acquisitions. In this study, we employed the SBAS method [3]. The SBAS method relies on (1) creating a stack of interferograms with small temporal and orbital separation (baseline) to limit the temporal and spatial decorrelation phenomena and (2) filtering out the topographic artifacts and atmospheric phase component using spatial and temporal information available in the processed data [3]. The key steps followed in the processing of Sentinel-1 and ALOS Palsar-1 datasets are provided in the following sections.

Sentinel-1 SBAS Time Series Analysis
The SBAS time series analysis method based on a linear deformation model was employed [3], and was implemented using SARscape software version 5.5.2 [39]. Neighboring same-track Single Look Complex (SLC) Sentinel-1 frames (Track 137; Frames 504 and 509) were mosaicked to cover the entire country of Qatar. ESA's Sentinel-1 precise orbits were used to refine the orbital information. The images were reduced in size to cover only the Qatar Peninsula (11,580 km 2 ). SLC image stacks were coregistered to a single reference image. A multilooking ratio of 24:6 (range:azimuth; ~90 m ground resolution) was applied to reduce the noise and the computational time. Differential interferograms were then created for image pairs having small temporal (< 365 days) and perpendicular baselines (< 130 m). The ALOS 30 m DEM was utilized to remove the topographic phase component. The Goldstein phase filtering method was applied to improve the phase unwrapping [40]. Phase unwrapping was performed for the pixels with coherence values higher than 0.3 using 3D Delaunay triangulation method [41]. The interferograms were refined by removing the residual phase gradients due to orbital errors, atmospheric signal or long wavelength deformation signal by applying second order polynomial estimation method using 265 ground control points (GCP) [38]. The phase unwrapping process was reapplied for the refined interferograms. Mean velocity field and possible topographic phase residuals were calculated in the first step of SBAS inversion with a linear deformation model. Topographic phase residuals were subtracted from wrapped interferograms. In the second inversion step, time series of deformation were estimated by the singular value decomposition (SVD) inversion [3]. The atmospheric phase component was estimated and removed by low-pass spatial and high-pass temporal filtering with filter sizes of 1200 m and 365 days, respectively [42]. The velocity and time series of deformation were retrieved and geocoded in geographic coordinates. In this step, different threshold values of several quality indices (temporal coherence, velocity precision, and height precision) were evaluated to identify and mask out pixels with low quality [39]. The selected threshold values for these indices were manually identified by step-by-step visual evaluation of the generated velocity maps and selection of the threshold values that produced reasonable coverage and minimized noiselike velocity patterns. The optimal threshold values identified for temporal coherence, velocity precision, and height precision were 0.35, 4 mm/yr, and 8 m, respectively. Following the completion of the initial SBAS processing with 71 scenes, the intensity and coherence images, wrapped and unwrapped interferograms, and velocity estimates were inspected for strong atmospheric artifacts, radio frequency interference (RFI) from ground sources, and unwrapping errors. A total of 34 out of 71 scenes were identified as being affected by one or more of the above-mentioned sources of error. Nine scenes exhibited RFI and produced interferograms with strong atmospheric artifacts, 13 scenes were affected by RFI, and 12 produced interferograms with strong atmospheric contributions. Particularly, scenes acquired after September 2019, suffered from frequent RFI artifacts. The final SBAS results were retrieved by reprocessing the data stack after omission of these 34 scenes. The SBAS network was created from 37 scenes and was composed of 177 interferograms. The mean temporal and perpendicular baselines for Sentinel-1 SBAS network were 82 days and 47 m, respectively ( Table 1). As described earlier, estimated LOS deformation was back-projected into the vertical direction by dividing by the cosine of the incidence angle with the assumption of zero lateral movement [38].

ALOS Palsar-1 SBAS Time Series Analysis
The Qatar Peninsula is covered by two ALOS Palsar-1 Tracks, 574 (Frames: 470, 480, 490, and 500) and 575 (Frames: 480, 490, and 500; Table 1 and Figure 2). Each ALOS Palsar-1 frame was analyzed individually using the SBAS procedures outlined in the previous section. Table 1 lists the number of scenes, number of interferograms, mean temporal and perpendicular baselines for each frame. In the development of ALOS Palsar-1 SBAS networks, we utilized the interferograms with small and long perpendicular and temporal baselines given that L-band interferometry is less affected by geometric and temporal decorrelation. A few of the long baseline interferograms displayed low coherence and were omitted from the SBAS processing. For ALOS Palsar-1 data, the multilooking ratio was kept at 6:29 and 12:29 (range:azimuth; ~90 m ground resolution) for the FBD and FBS imaging modes, respectively. The number of GCPs used for interferogram refinement ranged from 100 to 160 per ALOS Palsar-1 frame. For each frame, the LOS velocity and time series of deformation were retrieved and back-projected into the vertical direction by dividing by the cosine of the incidence angle assuming a purely vertical deformation field [38]. The derived vertical deformation products for each frame were geocoded in geographic coordinates then mosaicked to cover the entire Qatar Peninsula.

GPS Data
The InSAR ground deformation rates were compared to those extracted from Qatar's Continuously Operating Reference Station Network (CORS-Qatar) run by MME since 2009. GPS data was processed using precise point positioning (PPP) module of GIPSY-OASIS II v6.4 [42] developed by JPL at NASA. JPL final precise orbit and clock estimates for GPS constellation were used in the analysis. Three-component (east-west, northsouth, and up-down) daily position data was retrieved for the period 2009 to 2018. Linear vertical velocities were estimated at five GPS stations from daily position data after removal of the seasonal cycle. GPS-derived vertical velocities were compared with InSAR-based LOS velocities projected onto vertical direction assuming a purely vertical deformation field. The SBAS velocities over the GPS stations were extracted by averaging those of the SBAS results within a circle (radius: 150 m) centered over the station for comparison.

Precipitation Data
The satellite-based daily precipitation data [43] was retrieved from the Giovanni online data system, developed and maintained by the NASA GES DISC [44]. Annual precipitation was calculated from 1 July to next year's 30 June and centered on 1 January to include the rainy season, November to May [45].

Groundwater Level Data
The groundwater level data provided by Qatar Public Works Authority, Project CP761/1 Shallow Groundwater Monitoring Phase II is reported in meters in reference to Qatar National Height Datum (QNHD; ~1.3 m above Chart Datum [46]). Data reported before mid-2017 were collected manually, approximately a month apart, but the more recent observations were collected more frequently (daily) via a data logger system. Temporal variations in groundwater levels of the wells proximal to investigated sabkha areas were examined in relation to the observed deformation over the sabkhas in search of causal effects.

Surface Water Area within the Sabkhas
Temporal variations in surface water coverage within the sabkhas (Sabkha 1, Sabkha 2, and Sabkha 4) were calculated in Google Earth Engine [47] using Landsat-based monthly global surface water dataset (1984-2019) [48]. No surface water was detected within Sabkha 3 throughout the observation periods.

Spatiotemporal Correlation in a GIS Environment (Step III)
Deformation products extracted via InSAR time-series analyses based on ALOS Palsar-1 and Sentinel-1 datasets were included in a GIS database along with other auxiliary static (geologic and LULC maps, sinkholes, and depressions) and temporal datasets (highresolution optical satellite imagery, groundwater level, and surface water area within the sabkhas) for spatiotemporal correlations and data visualization.

Results and Discussion
The InSAR ground deformation maps extracted from ALOS Palsar-1 (January 2007-March 2011) and Sentinel-1 (March 2017-December 2019) datasets revealed that the majority of Qatar shows little or no deformation, yet subsidence and uplift with high vertical velocities of up to 35 mm/yr were observed locally (Figure 3a,b,d and e). The highest deformation rates were observed over sabkhas (salt flats) and landfills. The nature of deformation over sabkhas during the ALOS Palsar-1 and the Sentinel-1 periods was similarin this case uplift-yet different in magnitude, as was the deformation over the landfills. Landfills experienced subsidence during the ALOS Palsar-1 and the Sentinel-1 periods, but the subsidence rates were different. We counted a total of 10 inland and coastal sabkhas of various sizes within the study area, the most prominent of which are the inland sabkha of Dukhan and the coastal sabkhas of Mamlaha, Afjat Wadi Abu Al Ghurban, and Afjat Saadoon (Figure 3). High subsidence rates (up to 14 mm/yr) were observed over all major landfills (e.g., Umm Al Afai, Mesaieed), dumpsites (Rawdat Rashed and Umm Thanyatain), and the solid waste processing plant during the investigation period of both the Sentinel-1 and the ALOS Palsar-1 missions (Figure 3).
The histograms of deformation rates extracted from ALOS Palsar-1 and Sentinel-1 datasets over Qatar show similar Gaussian distributions with standard deviations of 1.7 and 2.7 mm/yr, respectively (Figure 3d). Figure 3e shows deformation rates along transect XX' (Figure 3c) and the associated root mean square error (RMSE). The transect XX' crosscuts the Dukhan Sabkha (Sabkha 1) undergoing uplift and the Umm Al Afai landfill (L1) undergoing subsidence. Examination of Figure 3e reveals larger RMSE values for deformation rates from ALOS Palsar-1 dataset compared to those from the Sentinel-1 dataset. Data from five GPS stations were used to validate SBAS results. All five GPS stations showed very low vertical deformation rates (<1 mm/yr) (Figure 3f). GPS-derived vertical deformation rates were consistent with the vertical deformation rates calculated from ALOS Palsar-1 and Sentinel-1 datasets (Figure 3g). The difference between GPS-and SBAS-derived velocities at all five GPS stations were under 3 mm/yr (mean: 1.4 mm/yr; standard deviation: 1.2 mm/yr).  (Figure 4h) within the sabkhas reveals an average uplift rate of 17 mm/yr for the four sabkhas during the Sentinel-1 period, approximately four times larger than the average rate (4 mm/yr) during the ALOS Palsar-1 period. Furthermore, the observed increase in the uplift rate during Sentinel-1 period was more prominent over the coastal sabkhas (Sabkhas 2, 3 and 4) compared to the inland sabkha (Sabkha 1). With the exception of Sabkhas 2 and 3 during ALOS Palsar-1 period, seasonal variations in deformation rates over the investigated sabkhas were not observed (Figure 5e). One explanation for the observed uplift throughout the ALOS Palsar-1 and Sentinel-1 investigation periods is the continuous deposition of evaporites in the surficial and near-surficial sabkha sediments due to evaporation of saline waters within the sabkhas. The coastal sabkhas are probably more affected by climate change-related sea level rise compared to the inland sabkhas. A rise in sea level will introduce more solutes in the neighboring coastal sabkha settings, and promote salt deposition, which in turn causes uplift within these sabkhas. Next, we investigate the validity of this hypothesis and the reason(s) behind the observed increase in uplift rates during the Sentinel-1 period.

Uplift in Sabkhas (Salt Flats)
The annual total precipitation, ground water levels, and monthly surface water area within the sabkhas show a noticeable increase in the recent years 2016-2020 (wet period) including the Sentinel  (Figure 5d), and in some wells periodic increases in the ground water levels are observed during rainy seasons (Figure 5c).
The above-mentioned observations are consistent with the following hypothesis: the increased precipitation during wet periods leads to a rise in groundwater levels to surface or near-surface levels, an increase in areas occupied by surface water within sabkhas, which in turn increases the rate of deposition of surficial evaporitic sediments. This hypothesis, if true, can explain the higher rates of uplift during the Sentinel-1 period compared to the ALOS Palsar-1 period (Figures 4 and 5). A similar uplift rate up to 15 mm/year was reported from salt flats of Atacama Desert in Chile using ENVISAT data and was attributed to accretion of evaporites within the capillary zone [49].
.  Figure 4a. (h) Boxplot for the deformation rates within the sabkha boundaries during the ALOS Palsar-1 and Sentinel-1 periods. The boundaries of the sabkhas were modified from the geological map of Qatar [35]. Time series of cumulative deformation averaged over each sabkha is given in Figure 5. Abbreviations are defined in caption of Figure 1c.  Additional insights into the evolution of the landfill over time and further clarifications of the above-mentioned observations could be made through comparisons between the DEMs at locations TS4 (northern section) and TS5 and TS6 (southern section) on Figure  6d. At location TS4 (northern section), a 6 m elevation difference is observed between SRTM and ALOS DEMs, but none between the ALOS and lidar DEMs, which indicates negligible, if any, filling operations since 2011. The subsidence rate at TS4 was found to be higher during the earlier ALOS Palsar-1 period (14 mm/yr) compared to the rate during the more recent Sentinel-1 period (9 mm/yr) indicating that the compaction rate decreased over time in the older sections of the landfill. At location TS6 (southern section), a 3 m elevation difference is observed between the ALOS and lidar DEMs, indicating more recent additions of fill material. At location TS6 (the most southern section), a high subsidence rate (up to 10 mm/yr) was observed during the more recent Sentinel-1 period compared to the near-stable conditions (subsidence of 2 mm/yr) during the ALOS Palsar-1 period. The southern section at large was coherent during the Sentinel-1 (subsidence rate at TS5 of up to 11 mm/yr), but not during the ALOS Palsar-1, period; this could be explained by active filling in this section during ALOS Palsar-1, but not during the Sentinel-1, period.
The observed subsidence is here attributed to mechanical compression and biochemical processes. The former refers to the compression of air-filled pores within the waste products, a continuous and long-term process caused by the compression of waste products, whereas the latter (biochemical processes) refer to the decomposition of organic matter [51]. Similar observations (subsidence over landfills) were reported from many landfills [52,53] worldwide.

Sinkholes
During our field trip in 2018, we visited three sinkholes (Dahl Al Hamam, Dahl Duhail, and Dahl Musfer): the first two were within the Doha city. Our field investigations in Dahl Al Hamam and Dahl Duhail showed clearly that the sinkhole in both areas occurred over caves that extend laterally beyond the boundaries of the sinkhole. Our observations are consistent with findings from earlier geophysical (gravity) investigations that revealed the presence of additional voids and subsurface extensions in and around the Dahl Al Hamam sinkhole [54]. The salinity of the water in the cave and its apparent tidal interconnection was interpreted as evidence for the presence of a wider network of interconnected caves extending some 4.5 km towards the Gulf [33,54]. These observations together with the ongoing anthropogenic activities that are changing groundwater levels such as wastewater injection and dewatering for large construction projects raised concerns about the stability of the sinkholes and cavities within the urban environments of Qatar and prompted us to investigate the stability of these sinkholes, especially the ones in urban environments.
Since the size of all the investigated sinkholes was smaller than the ground resolution of the InSAR products (90 m), deformation rate and time series were extracted over 12 individual pixels, each covering one of the investigated sinkholes. In Figure 7 and Table  2, the SBAS analysis revealed that all 12 sinkholes studied did not show apparent deformation. The Sentinel-1 data for the individual sinkholes showed less scatter than the ALOS Palsar-1 data. For example, the RMS of deformation estimates at the individual sinkholes ranged from 0.1 to 4 mm extracted from Sentinel-1 dataset compared to 1.0 to 12 mm derived from ALOS Palsar-1 dataset. We attribute this difference to lower sensitivity of L-band to deformation compared to C-band and limited number of scenes in the ALOS Palsar-1 data stacks. In any case, the deformation rates for both periods were within the uncertainty levels (+/−2 mm/yr) and analysis of cumulative deformation did not show any evidence for acceleration of deformation over any of the investigated sites. It is noteworthy that the deformation velocity and time series extracted from the pixels surrounding the sinkholes were also low and very similar to the estimated velocities over the central pixels reported in Table 2.  Although we did not observe significant deformation over any of the investigated sinkholes, these results should be considered with caution for a number of reasons. These include the use of coarse ground resolution (90 m) products and back-projecting the LOS displacement into the vertical direction assuming the absence of horizontal displacement; this assumption could be invalid for sinkholes still in formation or for expanding sinkholes that are experiencing significant horizontal deformational component. Furthermore, precursory deformation events leading to sinkhole formation can last weeks to years with progressively increasing subsidence rates reaching up to 1.5 mm/day prior to failure [12,13]. These short term, nonlinear precursory deformation events are unlikely to be detected by the SBAS method, particularly when a linear deformation model is adopted. The increase in the frequency of the SAR data acquisitions, the use of high spatial resolution deformation products, the adoption of nonlinear models, and the combined use of ascending and descending datasets can potentially improve the success of detecting the precursory deformation leading to sinkhole formation and expansion of existing sinkholes.

Limitations
Our approach has its limitations as well. We estimated deformation from the descending C-band Sentinel-1 dataset and the ascending L-band ALOS Palsar-1 dataset over two different time periods. This could introduce uncertainties in our attempts to compare the deformation rates estimated from both datasets given the differences in C-band and L-band sensitivity to deformation, look direction (ascending versus descending), and local incidence angle. To mitigate the differences in the look direction and incidence angle, we back-projected LOS deformation into the vertical direction before comparing deformation rates from the two datasets. This procedure allows comparisons to be made where horizontal deformation is negligible. We believe that the majority of LOS deformation over the study area is dominated by the vertical deformation given the flat topography of the study area, and the lack of any reported lateral tectonic deformational features.
Although GPS observations at five stations showed vertical velocities comparable to InSAR-derived vertical velocities, none of those stations were located over the areas experiencing high deformation rates (e.g., sabkhas and landfills/dumpsites), and no other independent field observations (e.g., GPS and leveling) were available over these areas to validate the observed high rates of deformation. Small-scale deformation (e.g., over sinkholes or buildings) are likely to have gone undetected in our analysis due to the coarse ground resolution (90 m) selected for the SBAS analysis; for sinkhole and building-scale deformation monitoring, methods employing persistent and distributed scatterers at finer ground resolution are recommended. C-band Sentinel-1 scenes over Qatar, particularly after September 2019, suffered from frequent RFI artifacts. Repeating the InSAR time-series analysis on smaller areas not affected by RFI artifacts will enable the utilization of larger stacks of Sentinel-1 scenes. Combining ascending and descending geometries, where available, should also be considered to retrieve 3-dimentional ground motion.

Conclusions
Qatar has been subject to rapid economic development and population increase since 2005. This rapid increase in population, coupled with the implementation of ambitious development plans and projects, was associated with progressive urbanization, LULC change, and groundwater level variations related to intense groundwater extraction and injection operations. In this study, we investigated whether these anthropogenic activities, together with natural forcings were responsible for the observed land deformation. We conducted SBAS time series analysis using SAR images (L-band ALOS Palsar-1 and Cband Sentinel-1) to generate the first map of ground deformation across the Qatar Peninsula.
We identified several sites that have been experiencing vertical ground deformation at rates of up to 35 mm/yr. The factors controlling the observed deformation were assessed by conducting spatiotemporal correlation between the deformation and the auxiliary static (geologic and LULC maps, distribution of sinkholes and landfill/dumpsites) and temporal (high resolution optical satellite imagery, precipitation, groundwater levels, surface water area within sabkhas) datasets.
The most extended, and clearly visible deformation is observed over active sabkha areas in southeastern and southwestern Qatar. The uplift in the sabkha areas (up to 35 mm/yr) is attributed to continuous deposition of evaporites in the surficial and near-surficial sabkha sediments due to the evaporation of saline waters within the sabkhas. High subsidence rates (up to 14 mm/yr) were observed over the Umm Al-Afai landfill and were attributed to mechanical compression and biochemical processes. The correlation and integration of radar-based deformation rates with elevation profiles extracted from DEM data sources enabled the unraveling of the evolution of the landfill through time, namely, its areal extent, and the nature and timing of the filling. The ongoing anthropogenic activities do not seem to be exacerbating the development of sinkholes. Areas in and around the known sinkholes did not show significant and/or accelerated deformation.
We presented the first countrywide deformation analysis over the country of Qatar in the Arabian Peninsula. We utilized current and earlier radar datasets and investigated controlling factors of the observed deformation by using relevant auxiliary datasets. We suggest that this study can pave the way to the establishment of similar countrywide studies over the remaining Arabian Peninsula countries and perhaps lead to the development of a ground motion monitoring system for the entire Arabian Peninsula in the future.