Land Subsidence Induced by Rapid Urbanization in Arid Environments: A Remote Sensing-Based Investigation

Abstract

The rapid increase in the population of many of the older major cities within the countries of the Saharan-Arabian Desert is steering vast and disorganized urban expansion and in many cases introducing adverse environmental impacts such as soil erosion, rise in groundwater levels, and contamination of shallow aquifers, as well as development of deformational features including land subsidence. Using the rapidly growing city of Riyadh (1992: 467 km²; 2018: 980 km²), the capital of the Kingdom of Saudi Arabia as a test site, we utilized Small Baseline Subset (SBAS) interferometric analyses of 2016 to 2018 Sentinel-1 images together with multi-temporal high-resolution images viewable on Google Earth, GPS, field, land use land cover (LULC), and geological data to assess the distribution and rates of land subsidence and their causal effects. Three main causes of subsidence were identified and assessed: (1) discharge of wastewater effluents from septic systems in newly urbanized areas that lead to an increase in soil moisture, rise in groundwater levels, waterlogging, and wetting and hydrocompaction of dry alluvium loose sediments causing land subsidence (up to −20 mm/y) in wadis and lowlands; (2) the subsurface dissolution of karst formation by wastewater effluents and the collapse of voids and cavities at depth under stresses introduced by heavy construction machinery, causing sagging and land subsidence (up to −5 mm/y); and (3) leveling, compaction, and degradation of municipal and building waste materials in organized landfills and disorganized dump sites that resulted in significant land subsidence (up to −21 mm/y) and differential settling that could jeopardize the stability of structures erected over these sites. Our findings highlight the potential use of the advocated integrated approach to assess the nature and extent of land deformation associated with rapid urban growth in arid lands, and to identify areas most impacted for the purpose of directing and prioritizing remediation efforts.

1. Introduction

Only 2% of the World population lived in urban centers in 1800, but 200 years later, these urban centers hosted some 50% of the World population and it is projected that by year 2030, 61% of the world population will reside in these centers [1]. Much of urban population growth is driven by migration from rural areas to benefit from the advantages offered by urban centers including, but not limited to, enhanced job opportunities, education, and health services [1].

Rapid urban growth has serious environmental repercussions including deterioration of air quality, degraded visibility, creation of urban heat islands, and elevated day and night temperatures [2,3,4]. Additional environmental impacts include inadequate water supplies, lack of sanitary landfills, and industrial pollution.

Many of the rapidly growing urban centers are witnessing high rates of subsidence. Using GPS, leveling, and interferometric synthetic aperture radar (InSAR) results, subsidence rates of about 5–10 cm/y and up to about 20 cm/y were reported from the cities of Jakarta, Bandung, and Semarang in Indonesia and were attributed to groundwater extraction, natural consolidation of alluvium soil, load of constructions, and tectonic activities [5]. In Hanoi city, in Vietnam, a mean vertical subsidence rate of 6.8 cm/y was extracted over three new urban districts using a multi-temporal InSAR analysis of ALOS images and was attributed to excessive groundwater extraction [6]. Similar findings were reported from Houston, Mexico City, Osaka, San Jose, Shanghai, Tokyo, and Venice, where subsidence was caused largely by pumping of ground water from unconsolidated sediment [7].

Less attention has been paid to assess the nature and magnitude of land deformation associated with rapid growth of urban centers in arid desert settings. Over the past few decades, many of the major cities across Saharan Africa and Arabia witnessed rapid expansion that continues to be fueled by rising populations and by the economic and cultural opportunities provided by these mega urban centers. Riyadh (To avoid confusion in this paper between the province of Riyadh and the city of Riyadh, we will refer to the former as “Riyadh Province” and the latter simply as “Riyadh.”), the capital of the Kingdom of Saudi Arabia (KSA), is one of these cities. The population of Riyadh increased from 0.4 million in 1970 to 4.5 million in 2005, and its area increased as well (24 km² in 1970 to 700 km² in 2005) [8]. Similar observations could be made for many other megacities in Saharan Africa (e.g., Cairo, and Khartoum) and Arabia (e.g., Baghdad, and Damascus). For example, Cairo, the capital of Egypt, had a population of 13 million in 1999 and 18 million in 2015, and its area increased from 423 km² in 1999 to 656 km² in 2015 [9]. Baghdad, in Iraq, had a population of a 5.3 million in 2008 and approached 7 million in 2019, as its area increased from 622 km² in 2008 to 727 km² in 2019 [10].

The economic development of the KSA has led to rapid urban expansion, with about 83% of population living in urban areas in 2015 compared to only 21% in 1950 [11,12]. In particular, large provinces such as Riyadh and Makkah in central and western Arabia, respectively (Figure 1), have witnessed the most dramatic urban expansion—almost a threefold increase in urban land use between 1992 and 2013 (e.g., 2340 to 6097 km² and 1616 to 3461 km² in Riyadh and Makkah provinces, respectively [13]). Riyadh, the largest urban center within the Riyadh Province is the study area of this research.

The rapid urban expansion and the recent discovery of deformation features (e.g., large sinkholes) in Riyadh [17] have drawn more attention to understanding the spatial distribution and mechanism of formation of deformation features in the city and to examining the role of human activities in their origin. This rapid urbanization was, in many sections of Riyadh, apparently disorganized. The majority of households had individual septic tanks in which the effluents were released in the ground probably causing dissolution of carbonate bedrock, sink hole formation and land subsidence, groundwater contamination, rise in groundwater levels, and discharge in lowlands. The significant increase of urban areas at the expense of rural lands over a short time period in Riyadh, have led to the development of poorly managed and randomly distributed dump sites in some areas. Natural hazards associated with urban encroachment over abandoned dump sites have been reported from numerous locations worldwide, including groundwater contamination and leakage of leachate in the city of Patras, Greece [18], and in Delhi city [19], and fatal slope failure in Quezon City, Philippines [20], and in Koshe landfill, Ethiopia [21]. Additionally, settlements established over abandoned dump sites are expected to cause land subsidence due to urban loading and waste degradation [22]. This study will assess the extent to which practices associated with rapid urbanization (e.g., release of effluents, poorly managed dump sites) could have caused land deformation (subsidence) in Riyadh.

Land subsidence is a slow phenomenon that refers to the gradual lowering of the surface of the land due to a reduction of the volume of subsurface materials and an increase in subsurface voids by either natural or anthropogenic processes [23]. Land subsidence is often related to urban expansion [24,25], sediment compaction [26], landslides [27], large-scale construction [28], mining activities, and mechanical compression and biochemical processes in dump sites [29]. Traditionally, land subsidence and associated deformation features has been examined using ground-based geodetic measurements including differential leveling, GPS, and LIDAR surveys [30]. Satellite-based measurements using interferometric synthetic aperture radar (InSAR) techniques provide a cost-effective and accurate tool for investigating land subsidence over large-scale areas [31]. InSAR techniques have been used widely to study land subsidence related to groundwater abstraction in different regions worldwide, such as northern California and southern Nevada [32,33,34], central and northern KSA [35,36], and in Beijing, China [37]. The Small Baseline Subset (SBAS) technique has emerged as an advanced method of InSAR time series analysis with algorithms that rely on small-baseline differential SAR interferograms to increase the volume of data for analysis [38,39].

In this study, we integrate field, remote sensing (e.g., InSAR and satellite imagery), land use and landcover (LULC), and hydrogeological data to accomplish the following tasks: (1) estimate the temporal and spatial deformation rates across Riyadh using SBAS techniques on Sentinel-1 data over a two year period (2016 to 2018), and (2) examine natural and anthropogenic factors that possibly control the observed deformation through spatial and temporal correlation of the observed deformation with relevant datasets (distribution of buildings, roads, dump sites, lithological units, wadis, and lowlands) in a GIS environment.

2. The Climatic, Geologic, Hydrogeologic, and Urban Settings of Riyadh

Riyadh extends in a more or less north–south direction and is largely built in the lowlands bordered on the west by the Tuwaiq Mountains and on the east by Hith Mountains (Figure 2). Figure 2 is a 3-D rendering of a Landsat 8 false color multispectral band image (spatial resolution: 30 m; acquisition date: 2020) draped over an ALOS World 3D digital elevation model (DEM; resolution 30 m). The Tuwaiq Mountains are bisected by southwest-northeast trending drainage networks (e.g., Wadi Mahadiyah, Wadi Wubay, Wadi Laban, and Wadi Namar; Figure 2) that drain into a main wadi (Wadi Hanifa). The latter (Wadi Hanifa) crosscuts the city in a northwest–southeast direction.

The Comprehensive Riyadh Strategic Plan (CRSP) was established to organize the expansion of the city with all its components (urban, environmental, economic, transportation, housing, and infrastructure). Four urban growth phases were set in the strategic plan between 2010 and 2030, with a start date and a defined extension for each of the four phases. Phase 1, 2, 3, and 4 start in 2010, 2015, 2020, and 2025, respectively, bringing the total area of Riyadh in these four phases to 2398 km², 2763 km², 3055 km², and 3131 km², respectively (Figure 1b; [16]). In the selection of the study area (Figure 1a,b; area outlined by red polygon), we visually identified (from Landsat 8 multispectral band image acquired in 2020) contiguous urban areas occupied by buildings and infrastructure.

The KSA in general including the city and province of Riyadh is characterized by desert climate with significant variations between day and night temperatures [40]. The annual average temperature is around 28 °C, while the maximum temperature in the summer can reach 50 °C and the mean precipitation in the Riyadh region is 87 mm/y [41].

The uplift associated with the Red Sea opening led to the exposure of the basement complex and the overlying Paleozoic sedimentary sequences on either side of the rift [14,42]. In central KSA, the outcrops of these formations, which range in age from Triassic to Miocene, are exposed near the Arabian-Nubian shield and dip steeply eastward towards the Arabian Gulf (Figure 1c; [14,43]). These formations include, from older to younger: (1) the shale, limestone, and sandstone of the Dhruma Formation; (2) the limestone and lime mudstone of the Tuwaiq Mountain Formation; (3) the muddy and lagoonal carbonates of the Hanifa Formation; (4) the lime mudstones, dolomite and sandstone of the Jubaila Formation; (5) the limestone, dolomite, and anhydrite of the Arab Formation; (6) the Hith anhydrite Formation; (7) the Sulaiy limestone Formation; (8) the Yamama limestone Formation; (9) the shale and dolomite of the Buwaib Formation; and (10) the Kharj Formation and alluvium deposits containing fragments of limestone, marl, shale, sandstone, chert, gravel, and gypsum. These stratigraphic units are covered locally by surficial deposits of gravels, sand, silt, and clays.

Karst-related phenomena have been observed and reported from large areas in the central and eastern parts of the KSA in carbonate and evaporite formations such as the Arab, Hith, Sulaiy, Umm ar Radhuma, Rus, Dammam, and Dam formations [44,45,46]. Most of the karstic features are controlled by joint systems [45] and were developed during previous wet climatic periods (Pliocene or Pleistocene; [44]). Riyadh is located within a dissolution-induced subsidence depression where thick Ca-sulfate rocks (up to 130 m) and carbonates dominate the subsurface [14]. The Riyadh urban area may have been the location of interstratal dissolution of karst units, leading to sagging and collapse of the overlying formations and the development of deformation features [47].

The faults and the folds reported in the area surrounding Riyadh indicate that they trend mostly in the northwest–southeast, northeast–southwest, west-northwest–east-northeast, and east–west directions (Figure 1). These faults and folds occur as deep-seated or shallow-to-surficial structures that have been initiated during the Late Cretaceous and experienced continued deformation until the Eocene [15,48].

Riyadh is underlain by the Mega Aquifer System (MAS), which is subdivided into two major aquifer systems, the Upper and Lower MAS. The significant aquifers for the Upper MAS are the Wasia-Biyadh and Umm er Radhuma; for the Lower MAS, the Minjur Aquifer is most significant. The Upper and Lower MAS are separated by the Upper Jurassic to Lower Cretaceous anhydrite aquitard of the Hith Formation [49,50]. These multiple aquifers within the MAS are partially hydraulically connected by subvertical fault systems which act as “hydraulic conduits” allowing cross-formation flow between the various aquifers [36,49,50].

The majority of households had individual septic tanks in which the effluents were released in the ground and concentrated sewage was collected from each household, treated in decentralized wastewater treatment plants, and used for landscaping, and agricultural irrigation, and industry. The release of effluents from septic tanks led to a general rise in groundwater levels in Riyadh, ponding in some locations, and discharge of groundwater into lowlands and wadis [51,52]. Moreover, encroachment of urbanization on major valleys (wadis) was a common practice. These wadis collect precipitation over extensive areas and channel it through the main valleys as runoff and groundwater flow in the underlying aquifers. Despite the paucity of precipitation, these main wadis (e.g., Hanifah and Laban) commonly become the site of flash floods [53].

In recent years, access to basic public services including electricity, drinking water, and sanitation has become a goal of the Ninth National Development Plan (NDP) and the country’s Vision 2030. Decentralized sewage systems are being replaced by centralized sewage systems and septic tanks are being decommissioned. In 2015, there was 55 wastewater treatment plants with a capacity of 1144 million cubic meters per year, and 42 more were under construction with an added capacity of 799 million cubic meters per year [54]. These sanitation projects will, with time, slow down the adverse impacts of the disorganized urbanization during the earlier periods, including land subsidence.

The significant increase of urban areas at the expense of rural lands over a short time period in Riyadh, has led to an annual addition of solid wastes of 2871 × 10³ tons and the development of randomly distributed dump sites in some areas [55]. These solid wastes were dumped outside the Riyadh city limits, compressed by heavy machinery, and covered with thick soil [56]. With time, the city limits expanded and some of the urban settlements could have been established over areas once occupied by dump sites, exposing residents of these areas to various hazards associated with dump sites such as groundwater contamination, release of greenhouse gases, and slope stability and settlement problems. Throughout the text, we investigate the extent to which the geologic, tectonic, climatic, and hydrologic settings and the practices associated rapid urbanization have contributed to land deformation in Riyadh.

3. Data and Methods

Our approach involves addressing two main tasks: estimating the subsidence rates within and surrounding Riyadh and understanding the formation mechanism of these deformation features. The first task—estimating of the temporal and spatial deformation rates over Riyadh and its surroundings—was addressed by extracting the radar line-of-sight (LOS) velocities using SBAS techniques on Sentinel-1 data (acquisition period: 2016–2018) and comparing these velocities to GPS data. The second task—investigating the natural and anthropogenic factors that possibly control the observed deformation—was conducted through spatial correlation of the observed deformation with relevant spatial datasets in a GIS environment such as the distribution of lithological units, wadis, and features displayed in temporal Google Earth images and LULC maps (e.g., buildings, roads, dump sites, and water-logged areas). Findings were corroborated by observations collected (2016–2020) by our team during our field campaigns and by information provided by the locals residing in the investigated areas.

3.1. Processing of InSAR Data and Extraction of Displacement Rates (Task I)

The SBAS algorithm is based on creating a stack of interferograms with small temporal and orbital baselines to reduce the temporal and spatial decorrelation phenomena, and compensating topographic artifacts and atmospheric phase component by spatial and temporal filtering using the information available in the processed data [57]. The European Space Agency (ESA)’s Sentinel-1A data were utilized, with a revisit time of 12 days in central Arabia. A series of 39 ascending Sentinel-1 scenes (

Table 1. Sentinel 1A SAR Acquisitions (ascending mode; Orbits: Path 72, Frame 73 and 77).

SR#	Mission	Orbit	Day	Month	Year	Comment
1	Sentinel (1A)	013969	16	11	2016	Scenes used
3	Sentinel (1A)	014319	10	12	2016
4	Sentinel (1A)	014494	22	12	2016
5	Sentinel (1A)	014669	03	01	2017
6	Sentinel (1A)	014844	15	01	2017
7	Sentinel (1A)	015019	27	01	2017
8	Sentinel (1A)	015194	08	02	2017
10	Sentinel (1A)	015544	04	03	2017
11	Sentinel (1A)	015719	16	03	2017
13	Sentinel (1A)	016069	09	04	2017
14	Sentinel (1A)	016244	21	04	2017
15	Sentinel (1A)	016594	15	05	2017
17	Sentinel (1A)	016944	08	06	2017
18	Sentinel (1A)	017119	20	06	2017
19	Sentinel (1A)	017294	02	07	2017
20	Sentinel (1A)	017469	14	07	2017
21	Sentinel (1A)	017819	07	08	2017
23	Sentinel (1A)	018169	31	08	2017
25	Sentinel (1A)	018519	24	09	2017
26	Sentinel (1A)	018694	06	10	2017
27	Sentinel (1A)	018869	18	10	2017
28	Sentinel (1A)	019044	30	10	2017
29	Sentinel (1A)	019219	11	11	2017
30	Sentinel (1A)	019394	23	11	2017
31	Sentinel (1A)	021144	23	03	2018
34	Sentinel (1A)	021669	28	04	2018
35	Sentinel (1A)	021844	10	05	2018
36	Sentinel (1A)	022019	22	05	2018
37	Sentinel (1A)	022194	03	06	2018
38	Sentinel (1A)	022369	15	06	2018
0	Sentinel (1A)	009944	14	02	2016	Scenes eliminated
2	Sentinel (1A)	014144	28	11	2016
9	Sentinel (1A)	015369	20	02	2017
12	Sentinel (1A)	015894	28	03	2017
16	Sentinel (1A)	016769	27	05	2017
22	Sentinel (1A)	017994	19	08	2017
24	Sentinel (1A)	018344	12	09	2017
32	Sentinel (1A)	021319	04	04	2018
33	Sentinel (1A)	021494	04	16	2018

) were selected and the SBAS radar interferometric technique [57] was applied using Sarscape software to create 255 interferograms (Figure 3). A subset of these scenes and interferograms were subsequently eliminated because they yielded interferograms with heavy atmospheric contributions and/or low coherence (red points; Figure 3). The remaining scenes (30 scenes; 16 November 2016–15 June 2018) and interferograms (159) were included in the remaining processing steps. The ALOS World 3D digital elevation model was used to remove topographical phase contributions. A multi-looking ratio of 5:1 (range:azimuth) was applied leading to square pixel of approximately 15 m ground resolution. Phase unwrapping was performed on the interferograms using Delaunay triangulation [58] for the pixels having coherence value higher than 0.35. Then, over 100 evenly distributed ground control points (GCPs) were selected and used for orbital refinement by estimating and removing the residual phase ramps on the unwrapped interferograms. The GCPs were selected from locations with high coherence and where no deformation is observed or expected. In the first inversion, based on a linear velocity model, possible topographic phase residuals and mean velocity field were calculated using the unwrapped phases. Topographic phase residuals were subtracted from the wrapped interferograms. In the second inversion step, the time series of displacement was estimated by the singular value decomposition (SVD) inversion approach using refined unwrapped interferograms [57]. The GCPs used for orbital refinement were also used as reference points in phase to displacement conversion step assuming that the deformation at the reference points were zero. Atmospheric phase components were estimated and removed by low-pass spatial (1200 m) and high-pass temporal (365 days) filtering [59]. LOS deformation velocities (Figure 4a;

Table 2. The distribution and rates of deformation of sites showing LOS subsidence exceeding 4 mm/y over Riyadh.

#	Site	Latitude	Longitude	Area	Subsidence * (mm/y)
1	W1	24°38′38.709″ N	46°31′51.384″ E	Wadi Wubay	−4 to −9
2	W2	24°41′44.044″ N	46°36′22.379″ E	Wadi Hanifah	−1 to −4
3	W3	24°39′32.734″ N	46°33′24.014″ E	Wadi Wubay	−2 to −8
4	W4	24°39′55.499″ N	46°34′18.458″ E	Wadi Wubay	−1 to −7
5	W5	24°39′5.564″ N	46°35′12.577″ E	Wadi Hanifah	−1 to −9
6	W6	24°38′27.016″ N	46°34′52.182″ E	Wadi Hanifah	−1 to −18
7	W7	24°38′41.857″ N	46°35′51.722″ E	Wadi Hanifah	−2 to −9
8	W8	24°36′21.642″ N	46°36′39.024″ E	Wadi Laban	−2 to −6
9	W9	24°35′55.154″ N	46°31′55.488″ E	Wadi Laban	−1 to −4
10	W10	24°34′15.424″ N	46°36′11.816″ E	Wadi Namar	−1 to −13
11	W11	24°35′20.698″ N	46°39′3.233″ E	Wadi Hanifah	−1 to −6
12	W12	24°33′40.369″ N	46°37′38.978″ E	Wadi Namar	−1 to −6
13	W13	24°35′44.25″ N	46°31′18.543″ E	Wadi Laban	−1 to −6
14	W14	24°50′58.932″ N	46°54′44.591″ E	Al Janadriyah	−3 to −5
15	L1	24°41′4.59″ N	46°53′8.083″ E	Khashm Al An	−1 to −7
16	L2	24°51′3.42″ N	46°37′52.625″ E	An Narjis	−2 to −4
17	L3	24°51′36.777″ N	46°37′49.09″ E	An Narjis	−1 to −4
18	L4	24°49′33.448″ N	46°39′44.518″ E	Al Yasmin	−1 to −7
19	L5	24°52′57.178″ N	46°38′13.64″ E	An Narjis	−3 to −5
20	L6	24°52′41.952″ N	46°37′22.343″ E	An Narjis	−3 to −5
21	L7	24°52′9.118″ N	46°38′11.129″ E	An Narjis	−3 to −6
22	L8	24°51′22.118″ N	46°36′46.65″ E	An Narjis	−1 to −7
23	L9	24°51′58.132″ N	46°37′30.828″ E	An Narjis	−3 to −5
24	L10	24°45′16.158″ N	46°37′56.401″ E	An Nakhil	−1 to −5
25	L11	24°33′44.12″ N	46°32′12.723″ E	Tuwaiq	−2 to −4
26	L12	24°33′43.813″ N	46°31′28.054″ E	Tuwaiq	−3 to −4
27	L13	24°32′53.479″ N	46°31′13.872″ E	Tuwaiq	−2 to −5
28	K1	24°51′41.453″ N	46°34′33.728″ E	Al Qirawan	−3 to −5
29	K2	24°51′46.361″ N	46°33′51.911″ E	Al Qirawan	−2 to −5
30	K3	24°51′8.338″ N	46°33′46.788″ E	Al Qirawan	−3 to −10
31	K4	24°50′23.807″ N	46°34′9.237″ E	Al Qirawan	−1 to −8
32	K5	24°49′50.841″ N	46°35′45.102″ E	Al Qirawan	−2 to −5
33	D1	24°46′0.131″ N	46°54′3.156″ E	Ar Rimayah	−1 to −6
34	D2	24°49′48.788″ N	46°43′1.198″ E	Qurtubah	−1 to − 4
35	LF1	24°31′24.093″ N	46°40′29.554″ E	Al Marwah	−1 to −21
36	LF2	24°29′44.341″ N	46°41′9.629″ E	Al Marwah	−1 to −11
37	LF3	24°37′32.879″ N	46°53′56.093″ E	As Sulay	−1 to −21
38	LF4	24°36′47.926″ N	46°47′8.627″ E	Al Manakh	−1 to −9

* Range of LOS velocities (mm/y) over each of the investigated sites.

) were retrieved and geocoded in geographic coordinates. In this step, pixels with the coherence value lower than 0.4 were masked. The deformation velocity (in mm/y) was contoured using a raster calculator in ArcMap for areas of subsidence and exported into the Google Earth domain to investigate the nature (urban, geologic, and hydrologic) of the setting undergoing subsidence and the factors causing the observed deformation through comparisons with temporal satellite imagery.

3.2. Processing of GPS Data to Validate InSAR-Derived Displacement Rates (Task I)

There are two GPS stations within the InSAR-based displacement mapping study. One of the stations is an IGS network station (SAS0; 46.639° E, 24.736° N). The daily solutions and vertical displacement velocity at this station is provided by Nevada geodetic laboratory [60] for the period between September 2019 and January 2021 (1.3 years). Vertical displacement velocity at this station is 2.28 ± 1.66 mm/y. Since the data acquisition period was not enough to acquire reliable long-term displacement velocity, this station was not included in the following discussion.

The second station (RY99, 46.694° E, 24.674° N; Figure 4) within the study area is part of Kingdom of Saudi Arabia Continuously Operating Reference Station (KSA-CORS) network operated by General Authority for Survey and Geospatial Information (CASGI). The time series of ~4.8 years covering the time-span 2015–2020 was processed using Bernese 5.2 [61] with the final IGS orbits, and satellite clocks and Earth orientation parameters. The processing included 18 permanent stations from different sources (International GNSS Services [IGS], EUREF Permanent GNSS Network [EPN], Scripps Orbit and Permanent Array Center [SOPAC], and University NAVSTAR Consortium [UNAVCO];

Table 3. GPS station position and velocities in east, north, and up components with uncertainties in mm/y relative to ITRF2014.

#	Stations	Latitude	Longitude	GPS Velocity mm/y			Standard Deviation mm/y
				Ve	Vn	Vu	$\sigma e$	$\sigma n$	$\sigma u$
1	ADIS	09°02′06.00″ N	38°45′57.59″ E	24.64	18.83	−3.20	0.27	0.17	0.48
2	ANKR	39°53′13.20″ N	32°45′28.80″ E	3.00	11.57	−4.43	0.59	0.22	0.79
3	BHR3	26°12′32.39″ N	50°36′28.79″ E	31.11	30.32	−0.18	0.10	0.08	0.31
4	BHR4	26°12′32.39″ N	50°36′28.79″ E	31.03	30.32	−0.10	0.09	0.06	0.57
5	BSHM	32°46′44.40″ N	35°01′22.80″ E	22.65	19.87	0.25	0.09	0.10	0.46
6	GRAZ	47°04′01.20″ N	15°29′34.80″ E	22.51	15.18	1.81	0.12	0.09	0.38
7	IISC	13°01′15.60″ N	77°34′11.99″ E	42.32	35.71	−1.19	0.16	0.11	0.52
8	ISBA	33°20′27.60″ N	44°26′16.80″ E	24.84	28.39	−0.36	0.15	0.19	0.57
9	MAL2	02°59′45.59″ S	40°11′38.40″ E	26.13	16.15	−0.42	0.12	0.15	0.62
10	MARS	43°16′44.40″ N	05°21′14.40″ E	20.08	15.85	−0.35	0.12	0.08	0.39
11	MAS1	27°45′50.39″ N	15°37′58.79″ W	17.13	16.92	−0.39	0.19	0.08	0.41
12	MBAR	00°36′03.60″ S	30°44′16.79″ E	25.58	17.91	0.62	0.11	0.17	0.49
13	RAMO	30°35′52.79″ N	34°45′46.79″ E	23.86	19.23	1.35	0.13	0.15	0.38
14	RY99	24°40′26.39″ N	46°41′38.40″ E	31.48	29.71	−0.15	0.08	0.19	0.34
15	SOFI	42°33′21.59″ N	23°23′41.99″ E	24.02	11.43	0.56	0.24	0.11	0.37
16	TEHN	35°41′49.20″ N	51°20′02.40″ E	26.82	19.14	−2.43	0.36	0.27	0.65
17	YEBE	40°31′29.99″ N	03°05′20.39″ W	19.27	15.96	0.75	0.16	0.06	0.28

), 13 of which are included in the ITRF2014 and were used as fiducial stations to attain a comprehensive configuration around the GPS station and for datum definition purposes. We applied the processing strategy described by [62,63]. The GPS data processing was carried out on a daily basis which resulted in a set of 1726 daily solutions in ITRF2014 reference frame [64]. The built-in combination tool (ADDNEQ2) in Bernese software was used for the estimation of station position and velocity for all the processed stations. An outlier detection/rejection process was applied followed by a second round of estimations for station position and velocity based on the clean solution. To estimate a realistic margin of standard deviations, we scaled the standard deviations by the ratio of the formal error and repeatability estimated by the combination step [63,65,66,67].

The LOS SBAS velocity over the GPS station was obtained over the pixel centered on the RY99 GPS station. The extracted GPS velocity in the vertical direction was projected to satellite geometry by multiplying by the cosine of the incidence angle at the GPS station assuming no, or negligible, lateral displacement and was compared to the LOS velocity extracted from the InSAR time series analysis.

3.3. Urban Land Use Changes between 1992 and 2018 (Task II)

Eight LULC maps for the years 1992, 1994, 1999, 2003, 2008, 2012, 2015, and 2018 extracted from European Space Agency Climate Change Initiative (ESA-CCI) land cover viewer datasets (http://maps.elie.ucl.ac.be/CCI/viewer/ accessed on 2 September 2020) were compared to estimate the rates of urban encroachment in Riyadh. These are multi-sensor datasets (i.e., ENVISAT MERIS, SPOT-VEGETATION, PROBA-V, and AVHRR images) that provide a 37 classes-unified global LULC coverage with 300 m spatial resolution [68]. These maps were validated by comparison to the GlobCover 2009 dataset [69]. Urban areas, with a class number of 190, was extracted as a separate thematic map for each of the investigated years. Results are plotted in Figure 5.

3.4. Spatial Analysis of Relevant Datasets (Task II)

All collected and processed data were co-registered to a unified geographic projection (datum: WGS-84; UTM Zone: N38) in a GIS environment. Geological, topographic, and land use and land cover maps together with high-resolution images viewable on Google Earth were used to monitor the environmental impacts associated with urban expansion. These include the rise in groundwater level and the increase in soil moisture, waterlogging, and vegetation cover. The integrated datasets were also used to identify areas most susceptible to these adverse impacts. The integration of observations from various products provided valuable clues pertaining to the formation mechanisms of observed land deformation in Riyadh.

4. Results

Inspection of Figure 5 reveals that the area occupied by urban development within our study area doubled between years 1992 (467 km²) and 2018 (980 km²). Temporal Google images show that this urban encroachment was largely on Riyadh’s western and northwestern sections over the foothills of Tuwaiq Mountains and the interleaving valleys and lowlands (Figure 2) and less so on the eastern and southern sections. A similar pattern is observed on the LOS land deformation SBAS map (Figure 4a); areas experiencing subsidence exceeding 4 mm are largely concentrated on the western and northwestern sides, of the city where 32 of the 38 identified subsidence sites are located. These observations suggest that the observed LOS SBAS land deformation is largely caused by anthropogenic activities related to urban expansion. We first compare the LOS SBAS velocities to GPS data to validate or calibrate the SBAS measurement then conduct spatial correlation of the observed LOS SBAS deformation with relevant temporal and spatial datasets in a GIS environment to assess the observed deformation and to identify factors controlling it.

4.1. Comparison of SBAS and GPS Velocities

The east and north velocity components of the GPS measurement at the RY99 GPS station reported in Table 3 represent the northeast movement of the Arabian Plate towards Zagros Mountains in Iran [70,71]. It is safe to assume that this lateral plate motion movement is constant within the study area given that the area is located in the center of the Arabian Plate. These constant displacements (lateral or vertical) affecting the entire extent of the processed SAR data with the same magnitude (rate) will not be detected by the InSAR method. Instead, the LOS InSAR velocities describe the observed relative deformation within the extent of the processed SAR data. Hence, the observed SBAS deformation accounts only for the relative deformation within the processed area.

We compared the LOS velocities from the SBAS analysis (1.29 ± 0.2 mm/y) during the investigated period to those vertical GPS velocities (−0.12 ± 0.27 mm/y) projected to satellite LOS direction derived from the Riyadh RY99 GPS station (Figure 6). The reported error for the LOS SBAS velocity over the GPS station was estimated using a statistical approach that accounts for trend errors in the time series [72,73]; specifically, Monte Carlo simulations were adopted by fitting trends and other terms for many (n = 10,000) synthetic datasets and the standard deviation of the extracted synthetic trends was interpreted as the trend error. The estimated statistical error is an underestimate of the overall uncertainties in LOS SBAS velocities because it does not fully account for the additional uncertainties related to SBAS processing (e.g., atmospheric artifacts, unwrapping errors, and residual topography). Those uncertainties could exceed 2 mm/y for the number of utilized Sentinel-1 scenes [74]. Thus, the radar-based velocities are not statistically different from GPS velocities within 2 sigma error and no calibration for the radar-based velocities was conducted.

4.2. Subsidence Over Wadis and Lowlands

Inspection of the deformational velocities displayed in Figure 4a reveals moderate to high subsidence rates (up to −21 mm/y) over highly localized areas in the western sections of the city. Inspection of digital elevation maps over these areas reveals that many of these areas are concentrated within wadis and lowlands in the western sections of the city, but much less so on the eastern side (Figure 4a,b). Examples of subsiding areas within wadis include sites W2, W5, W6, W7, and W11 in wadi Hanifah or in wadis draining in, and proximal to it (wadi Hanifa), sites W1, W3, and W4 in wadi Wubay, site W8, W9, and W13 in wadi Laban, sites W10 and W12 in wadi Namar on the western side, and site W14 (Al Janadriyah area) on the eastern side (Figure 4). Examples of subsiding areas within lowlands include sites L2, L3, L5, L6, L7, L8, L9 in An Narjis District, L4 in Al Yasmin District, and L10 in An Nakhil District on the northern side, site L1 in Khashm Al An District on the eastern side, and sites L11, L12, and L13 in Tuwaiq District on the western side of Riyadh (Figure 4).

Inspection of temporal (2001–2018) Google Earth images over the above-mentioned subsiding areas indicates that the majority of these areas show evidence for increasing urbanization and associated increase in the soil moisture with time (e.g., sites from L1 to L13 and W5, W6, and W7; Figure 4). The increase of soil moisture can be inferred from the increased vegetation and waterlogging along these valleys and lowlands with time. As soil moisture increases, the soil darkens in tone, and approaches low reflectance values (dark grey to black tones) in waterlogged areas. The observed tonal variations in the wadis are local in nature and should not be confused with the evenly distributed spectral effects caused by seasonal variations in scene illumination or by precipitation-related variations in soil moisture. Both of these seasonal effects will affect the overall brightness of the individual wadis throughout their length and the brightness of the wadi network across the entire scene. In a number of these wadis and lowlands where waterlogging has been problematic there have been attempts to fill the waterlogged areas with loose sediments (e.g., sites W2, W5, W6, W7, W8, W10, W11, and W12; Figure 4).

Figure 7 shows LOS deformation rates (2016–2018) over subsiding areas within Wadi Hanifah (site W5) and displays temporal Google Earth images (March 2001–October 2016) and ArcGIS online base maps that show landscape evolution and progressive urbanization over time (2001: Figure 7a; 2005: 7b; 2014: 7c; 2016: 7d; 2019: 7e). Inspection of Figure 7 shows that in 2001, the site was barren and was dissected by an east–west trending wadi network (Figure 7a). The area witnessed the onset of urbanization by 2005 (Figure 7b), and the appearance of vegetation and burial of sections of the main wadi by 2014 (Figure 7c). Comparison between Figure 7c,d reveals that the progressive urbanization between 2014 and 2016 was associated with waterlogging and increased vegetation. Deformation rates of up to −9 mm/y were observed over the levelled and buried sections on the main wadi in Figure 7e. The LOS deformation rates were plotted over the 2019 ArcGIS online base map. Although not shown, this scenario is typical of sites W1, W3, W4, W9, W13, and W14 (Figure 4).

Similar high LOS subsidence rates and landscape evolution were observed over lower elevation areas, hereafter referred to as lowlands, that witnessed progressive urbanization through time. Figure 8 shows LOS deformation rates over subsiding areas within a lowland in An Narjis District (site L7) and displays temporal Google Earth images (April 2008–November 2018) and ArcGIS online base maps that show evolution of the landscape and progressive urbanization through time. Site L7, a barren area, was divided into parcels by 2008 (Figure 8a) and witnessed progressive urbanization starting in 2012 (Figure 8b). The progressive urbanization caused increased soil wetness and waterlogging, which was observed in the 2016 image (Figure 8c) and apparently required covering by loose soils observed in the 2018 image (Figure 8d). Deformation rates of up to −6 mm/y were observed over the central sections in Figure 8e, where the LOS deformation rates were plotted over the 2019 ArcGIS online base map. Although not shown, this scenario is typical of site L10. In the preceding examples (sites W5 and L7), and in other subsiding lowlands and wadi sections, the urbanization extended over these areas (e.g., sites W1 to W4; W6 to W14; L1 to L6; and L8 to L13; Figure 4). Not all of the subsidence within wadis and lowlands was associated with progressive urbanization. Exceptions include W14, a wadi that is distant from, and upstream of, urban areas.

We attribute the above-mentioned subsidence in lowlands and wadis to release of effluents from septic systems into both alluvial and karstic aquifers. Until recently the overwhelming majority of urban developments had decentralized sewage systems that utilize septic tanks, where the effluents were released in the underlying alluvial or karstic aquifers. It was reported that sewerage system in Riyadh serve only 50% of urban area population, while the remaining 50% uses lined and unlined septic tanks and shallow wells for wastewater disposal [54,75]. As described earlier, these practices led to increase in soil moisture and rise in groundwater levels in Riyadh, discharge of perched groundwater or groundwater in the lowlands and wadis, waterlogging in these areas, and their burial by loose sediments. An average rate of rise in water level of about 0.55 m/y was reported from eastern Riyadh [51]. The wetting of dry, unconsolidated, porous deposits flooring the wadis can cause hydrocompaction [76], a phenomenon that causes subsidence in such areas [76].

Hydrocompaction associated with wetting of loose sediments in the dry wadis or those used to bury the waterlogged areas is here interpreted to be causing the observed high rates of subsidence described above. The wetting of dry and loose sediments could be accomplished by flash floods (e.g., W14) in wadis and by sewage effluents in waterlogged areas (e.g., W5, W6, W7, W9, and W13; Figure 4) and the filling of waterlogged areas within wadis and lowlands with loose sediments will have the effect of broadening the areas affected by the advocated hydrocompaction-related subsidence.

4.3. Subsidence Over Karstic Topography

Figure 4a reveals additional locations that display moderate subsidence rates of up to −9 mm/y (e.g., K1, K2, K3, K4, K5, D1, and D2: Figure 4) to high subsidence rates of up to −21 mm/y (e.g., LF1 to LF4; Figure 4). Inspection of temporal Google Earth images together with our field observations over all these areas reveals that unlike the previously described locations (W1 to W4 and W6 to W14; L1 to L6 and L8 to L13: Figure 4) these sites do not display many of the features typical of a hydrocompaction-related origin. These sites are not found in lowlands or wadis, do not show evidence of progressive increase in urbanization, vegetation, or waterlogging, and do not show an increase in soil moisture with time. A number of these locations (K1, K2, K3, K4, and K5; Figure 4) are found in the Al Qirawan District and lie over and/or are underlain at depth by formations characterized by karstic topography (limestone, dolomite, and anhydrite), namely, the outcrops of the Upper Jurassic Jubalia Formation, the highly fractured limestone and evaporite of the Jurassic Arab Formation, and the dolomite limestone of the Cretaceous Sulaiy Formation. Solution features (voids, open fractures, side wall cavities, caves, scattered collapsed and buried sinkholes, and solution caverns) resulting from chemical leaching of carbonate and evaporite rich formations were reported from these formations using a wide range of techniques and disciplines. These include field, geophysical (electrical resistivity, seismic refraction, and ground penetrating radar) investigations, core drilling, in situ permeability tests (i.e., loss of drilling water and backfill grout), mechanical probes (rate of drilling in relation to depth), geotechnical studies [17,77,78,79], and InSAR analyses [36].

Inspection of temporal Google Earth images (January 2014–December 2018) over one of these areas (site K1) shows that until 2014, this site was a barren land formed largely of low-lying carbonate hills (height: 3 to 5 m above surroundings) that are dissected by wadis (Figure 9a). The removal of the hills, levelling, and subdivision of the area into parcels apparently began around 2016 (Figure 9b), proceeded in 2017 (Figure 9c) and 2018 (Figure 9d), and was completed in 2019 (Figure 9e). By then all the hills were removed and the terrain was levelled and subdivided into parcels in preparation for development of large urban projects. These constructional activities involved the operation of heavy machinery, the addition (and later wetting) of a layer of sediment (thickness: from 1 m up to 15 m), and compaction of the surface and substrate material by heavy machinery [80]. One plausible explanation for the observed deformation is that the compaction and wetting of sediments and/or the collapse of the voids and cavities at depth through the use of heavy machinery [80,81] could have produced collapse of paleo-sinkholes and voids, sagging in the overlying layers (e.g., [82,83,84]), and the LOS subsidence (up to −5 mm/y) observed in Figure 9e. However, we cannot ensure that the observed subsidence is characteristic of the (undisturbed) surface because the period during which the above-mentioned surface-shaping activities occurred coincided with the period of SBAS study (2016–2018).

4.4. Subsidence Over Dumpsites and Landfills

Unlike all of the previously described sites that are located over limestone formations, D1 (Ar Rimayah District) and D2 (Qurtubah District) lie over the Quaternary alluvial gravel, sand, and silt deposits of Wadi AlKharj (Figure 4b,c). Our field observations and inspection of temporal Google Earth images (May 2009–November 2014; Figure 10) reveal that these two sites have been used as unorganized dump sites. Up to 2009, site D1 was barren land formed largely of low-lying alluvium sediments (Figure 10a). By 2011, the area was already being used as a dumping site for building and municipal solid waste materials (Figure 10b) and continued as such into 2013, where a significant increase in the spatial extent of the dump site was observed (Figure 10c). By 2014 the use of the area as a dump site ceased. Instead, the southern sections of the dumpsite were levelled, and a submeter-thick soil layer was added and compacted by heavy machinery; construction proceeded of urban settlements atop the buried sections of the dump site (Figure 10d). The process of erecting structures atop the buried dumpsite continued up to 2019, as shown in Figure 10e, which also shows LOS subsidence of up to −6 mm/y over the buried sections of the dumpsite. The figure also shows varying subsidence rates (−1 mm/y to −6 mm/y) over the buried areas within the dump site (Figure 10e), an observation that suggests that the compacted materials could be subsiding at different rates over relatively small distances. If true, this differential settling could jeopardize the integrity of the overlying constructions.

Sites LF1 (area: 0.80 km²; Al Marwah District; Figure 4), LF2 (area: 0.38 km²; Al Marwah District; Figure 4), LF3 (area: 0.55 km²; As Sulay District; Figure 4), and LF4 (area 0.1 km²; Al Manakh District; Figure 4) were described as locations of organized landfills that were in place prior to year 2000 [56,85,86,87]. Landfills LF1, LF2, and LF3 received waste material from farms, houses, companies, streets, hospitals, workshops, factories, and construction sites, whereas landfill LF4 received industrial waste from a cement factory [56,85,86]. A common practice in these landfills is to spread (with heavy machinery) and compress the waste products to a thickness of 2.5 m, then cover the compressed waste by a submeter-thick layer of loose soil [85].

Figure 11 shows high LOS subsidence rates over sites LF1 (up to −21 mm/y) and LF2 (−11 mm/y). Although not shown, the subsidence rates over LF3 and LF4 are also as high (LF3: up to −21 mm/y; LF4: −9 mm/y, respectively). We attribute the observed high subsidence rates over L1, L2, L3, and L4 to mechanical compression and biochemical processes. Mechanical compression refers to compression of air-filled pores within the waste products, a continuous and long-term process that is caused by compression of low-density waste products. The biochemical processes refer to the continuous decomposition of organic matter [88,89,90]. Similar observations (subsidence over landfills) were reported over landfills (e.g., Noeul and Haneul parks in Seoul, Korea; [90]).

5. Discussion

The majority of the rapidly growing cities in arid lands were not intended to accommodate as many inhabitants as they do today, nor were they expected to expand as much in area. In this respect, the environmental problems associated with this widespread and rapid urban expansion described in this work were not expected, nor were their remedies adequately planned for ahead of time. For example, while the effluents discharged from septic systems in Riyadh did not pose significant environmental problems in the early 1990′s, they do now, as they were found to increase soil moisture and cause waterlogging and subsidence in wadis and lowlands. Similarly, the earlier unregulated dump sites outside the boundaries of the old city did not pose environmental problems at the time, but decades later, these dumpsites are located at the outskirts of, or incorporated within, the new borders of the expanded city, making them potential sites for future construction that could be subject to differential settling.

Integrated studies similar to the one carried out here can provide insights into the ongoing and earlier ground deformation in rapidly growing urban centers in arid lands, investigate the driving forces for the observed deformation, and identify areas most impacted for the purpose of directing and prioritizing remedial efforts. For example, the areas that should be among the first to be switched from septic to centralized sewage systems are those showing high subsidence rates, increase in soil moisture, and are subject to waterlogging. Furthermore, our approach could be used to identify buried dump sites and/or karst undergoing subsidence for the purpose of cautioning against, or adhering to appropriate building codes for, new construction over these sites.

Our approach has its limitations. The LULC observations cover a much longer time span (2001–2019) than that covered by the stack of Sentinel-1 scenes that were used in this study (2016–2018), and thus we could not assess the variations in deformation rates associated with LULC throughout the entire period covered by Google Earth images.

6. Conclusions

Rapid, and in many cases, disorganized urban expansion of the older major cities within the countries of the Saharan-Arabian Desert often causes adverse environmental impacts, one of which is land subsidence. Using Riyadh city, the capital of Saudi Arabia as a test site we assessed the magnitude and nature of land deformation associated with rapid urban growth over the study area and which could be applicable to many of the megacities within the Saharan-Arabian Desert. We adopted an integrated approach that relies on the extraction of deformation rates (period: 2016 to 2018) from Sentinel-1 images using SBAS interferometric analyses, verified the SBAS-deformation rates against GPS measurements, examined LULC variations over deformed areas using multi-temporal high-resolution images viewable on Google Earth in search of causal effects, and collected field observations to verify the inferred factors (anthropogenic versus natural forcings) causing subsidence.

Thirty-eight sites experienced subsidence rates exceeding −4 mm/y and were found over stretches of wadis and lowlands, areas characterized by karstic topography, and over dumpsites and landfills. Subsidence (up to −20 mm/y) over wadis and lowlands was here attributed to discharge of wastewater effluents from septic systems in newly urbanized areas that lead to an increase in soil moisture, rise in groundwater levels, waterlogging, and wetting and hydrocompaction of dry and loose alluvium sediments. Subsidence (up to −5 mm/y) over areas characterized by karstic topography was related to dissolution of karst-bearing formations by wastewater effluents and possibly the collapse of voids and cavities at depth under stresses introduced by heavy construction machinery. Finally, high subsidence rates (up to −21 mm/y) and differential settling were observed over landfills and dump sites due to leveling, compaction, and degradation of municipal and building waste materials. The latter (differential settling) could jeopardize the stability of structures, if erected over these sites.

Our findings illustrate a wide range of adverse environmental impacts due to rapid and widespread urbanization of urban centers in arid environments and highlight the potential applications of the advocated integrated methodology (SBAS interferometric analyses, GPS measurements, LULC variations from multi-temporal high-resolution images viewable on Google Earth, field observations) in addressing these environmental issues. These include assessing the distribution and magnitude of land subsidence associated with rapid urban growth in arid lands, investigating the cause of the observed subsidence, and identifying areas most impacted for the purpose of directing and prioritizing remedial efforts.