Assessing the Impacts of Groundwater Depletion and Aquifer Degradation on Land Subsidence in Lahore, Pakistan: A PS-InSAR Approach for Sustainable Urban Development

Abstract

In various regions worldwide, people rely heavily on groundwater as a significant water source for daily usage. The resulting large-scale depletion of groundwater has triggered surface deformation in densely populated urban areas. This paper aims to employ Persistent Scattered Interferometry Synthetic Aperture Radar (PS-InSAR) techniques to monitor and quantify the land surface deformation (LSD), assess the relationships between LSD and groundwater levels (GWL), and provide insights for urban planning in Lahore, Pakistan, as the research area. A series of Sentinel-1 images from the ascending track between 2017 and 2020 were analyzed. Moreover, the Mann–Kendall (MK) test and coefficient of determination were computed to analyze the long-term trends and spatial relationships between GWL depletion and line of sight (LOS) displacement. Our findings reveal significant increases in land subsidence (LS) and GWL from 2017 to 2020, particularly in the city center of Lahore. Notably, the annual mean subsidence during this period rose from −27 mm/year to −106 mm/year, indicating an accelerating trend with an average subsidence of −20 mm/year. Furthermore, the MK test indicated a declining trend in GWL, averaging 0.49 m/year from 2003 to 2020, exacerbating LS. Regions with significant groundwater discharge are particularly susceptible to subsidence rates up to −100 mm. The LS variation was positively correlated with the GWL at a significant level (p < 0.05) and accounted for a high positive correlation at the center of the city, where the urban load was high. Overall, the adopted methodology effectively detects, maps, and monitors land surfaces vulnerable to subsidence, offering valuable insights into efficient sustainable urban planning, surface infrastructure design, and subsidence-induced hazard mitigation in large urban areas.

1. Introduction

Water resources are a significant source for living organisms [1,2]. Environmental factors, population growth, rapid socioeconomic progress, and intensified agricultural and industrial practices have resulted in growing water demand under massive stress in many cities worldwide [3,4]. Moreover, as a result of the existing circumstances, due to substantial growth in population and urban development, the demand for groundwater has dramatically increased in recent decades [5]. It is estimated that approximately 982 billion cubic meters (BCM) of groundwater are withdrawn globally each year, 70% of which is used for agriculture [6]. In addition, climate change is causing reduced surface water supplies in certain areas, which leads to the increased pumping of groundwater to meet demand [7]. Groundwater and other underlying natural fluids play significant roles in maintaining the Earth’s aquifer system [5]. The large-scale excessive extraction of groundwater reduces pore pressure and causes compaction, resulting in land surface deformation (LSD) and sinkhole formation [5].

This was coupled with a lower recharge of underground aquifers, which resulted in groundwater depletion. When the rate of groundwater extraction exceeds the rate of replenishment, land subsidence (LS) and sinkholes may occur, posing a significant threat to infrastructure and the environment [8,9]. The increasing global population and uncontrolled groundwater extraction exacerbate these hazards worldwide [10,11,12,13]. LS is frequently reported worldwide, particularly in urban areas with high population densities and coastal regions. The overexploitation of groundwater has significant consequences in all major urban areas, especially in developing countries, such as Pakistan, where environmental legislation has not been adequately implemented. Groundwater is an essential resource for Pakistan, and it is estimated that approximately 60–70% of the total population in Pakistan relies directly or indirectly on groundwater for their livelihood [14].

Pakistan has exceeded the sustainable limit of safe yield by extracting 61 billion m³ of water (50 million acre-feet; MAF) from aquifers [15], while the demand for domestic groundwater increasing from 5.20 to 9.70 million acres [16]. The Indus Basin aquifer, which is vital to the country’s water supply, is ranked as the second most overstressed groundwater reserve globally. Over 60% of irrigation, 90% of urban water supply, and 100% of industries in the country depend on groundwater [14]. The over-exploitation of groundwater in densely populated areas with extensive confined aquifers leads to LS [13,17,18]. Pakistan, with a population of over 175 million and a population growth rate of 2.1%, faces this problem, particularly in urban areas such as Lahore. During the 1950s, the groundwater level (GWL) in Lahore, Pakistan, was shallow, at approximately 4.58–4.88 m. However, because of rapid population growth and urbanization, the GWL has now dropped to 21.34 m, declining by an average of 15.24 m from 1959 to 1989 [19]. To access groundwater for domestic, agricultural, and industrial purposes, the water table is depleting at an alarming rate. There is an urgent need to assess the impact of groundwater depletion and aquifer degradation on LS in these urban regions of developing countries. Recent advancements in radio detection and ranging (RADAR) interferometry, specifically Interferometric Synthetic Aperture Radar (InSAR), have been demonstrated to be both efficient and cost-effective for monitoring LSD and LS [20]. InSAR is a powerful technique that enables the measurement of surface topography and its temporal changes [21], particularly in relation to the geological and hydrological factors that characterize aquifer systems [22].

However, InSAR has limitations due to its spatial and temporal decorrelation, atmospheric delays, and acquisition intervals. Persistent Scatterer Interferometric Synthetic Aperture Radar (PS-InSAR) was developed to address these issues and improve subsidence assessment by reducing its ranges from 10–20 mm to 2–3 mm [23]. Advanced methodologies, such as PS-InSAR, are multi-temporal InSAR methods that allow the highly accurate recovery of ground displacement at millimeter precision in urban areas [24,25]. The presence of artificial structures and buildings increases the density of Persistent Scatterer (PS) locations and significantly improves the signal-to-noise ratio of interferograms, making PS-InSAR particularly advantageous for LS detections. In general, the previous studies [26], focused on land subsidence by groundwater extraction in Lahore.

This study aims to investigate the potential environmental threats caused by LS resulting from large-scale groundwater extraction and recommend mitigation strategies for groundwater management challenges, and it represents a substantial advancement in the field of land subsidence analysis in Lahore, Pakistan. It utilizes diverse datasets of long-term InSAR measurements to capture detailed line of sight (LOS) displacement over time from 2017 to 2020. The results of this study would have important implications for effective planning, risk management, subsidence hazard management, and water management strategies in the research area. Additionally, this study sheds light on the extent and nature of subsidence in similar urban areas worldwide, which can help policymakers and urban planners make informed decisions.

2. Study Area

Lahore is located in the eastern part of the country in the Punjab province in the upper Indus plain on the Ravi River, a tributary of the Indus. It is the second largest city in Pakistan and is located in the capital city of Punjab the province [27]. It is situated within the geographical coordinates of 31°13′–31°43′N latitude and 74°0′–74°39.5′E longitude, at an elevation of 150 to 200 m above sea level. The city covers a total area of 1842 sq. km, and is located on the left bank of the Ravi River (Figure 1). Lahore has a warm semi-arid steppe climate characterized by significant fluctuations in both temperature and rainfall. The city experiences winters from November to March, with the coldest months being December and January, with temperatures ranging from 0 °C to 3 °C. The summer season continues from April to September, with an annual mean temperature varying from 18 °C to 38.8 °C. The city is entirely groundwater-dependent and has experienced rapid urbanization and migration of people to Lahore [28]. According to Pakistan bureau of statistics from the 2017 census, the total population of Lahore reached 11.2 million, with an annual growth rate of 3% per year (https://www.pbs.gov.pk/ last accessed: 13 February 2023). This growth led to the division of the city into nine administration zones managed by the Town Municipal Administration (TMA), along with the Lahore cantonment area. Lahore’s population is primarily concentrated in urban areas, with more than 98% of the population living there.

The city’s average annual rainfall is approximately 628.7 mm, with July and August being the wettest months, contributing approximately 40 mm of rainfall, which aids in groundwater recharge in a normal year [29]. The city of Lahore in Pakistan extracts 6.31 million cubic meters of groundwater per day to fulfill its urban, industrial, and agricultural water supply needs. However, the excessive extraction of groundwater combined with land development and decreased aquifer recharge have caused a significant decline in the groundwater table [30]. Currently, the GWL has declined to more than 30.48 m at many locations owing to overexploitation, which has resulted in LS and socio-environmental issues.

Geological Conditions

The Lahore region is mainly characterized by sedimentary rocks, which are overlaid by a significant thickness of alluvial deposits, reaching up to 300 m in depth [31]. They were formed millions of years ago when the area was covered by shallow seas and rivers [32]. These rocks are of different ages and formed during different geological periods. In contrast, alluvial deposits were formed by the Indus and its tributaries, which carried sediment from the Himalayas and deposited it in the Punjab plain (Figure 2). The deposits consist mainly of sand, silt, and clay, and are fertile, making Lahore a major agricultural region [33]. Additionally, Lahore lies in a seismically active zone with the potential for earthquakes because of its location on the Indian plate. However, as it is primarily composed of soft alluvial soil, the region is susceptible to geohazards [34]. The presence of old channels in the Ravi River indicates that there have been periods of excessive flooding during which earlier channels were filled with sediment and the stream was forced to create new channels. This abrupt migration indicates that the stream had to adapt to changes in its environment, such as increased water flow or sediment load.

3. Materials and Methods

3.1. Remote Sensing SAR Dataset

In this study, we utilized Synthetic Aperture Radar (SAR) data obtained from Sentinel-1A satellites operating in Interferometric Wide Swath (IW) mode to acquire SAR images with vertical co-polarization (VV) and C-band (5.6 cm wavelength). The ground resolution of each image was approximately 4 m × 14 m in the range and azimuthal directions, covering an area of approximately 250 km × 200 km. The satellite has a revisit time of 12 days. A total of 128 images were used from the ascending orbit tracks of the Alaska Satellite Facility (https://asf.alaska.edu/about-asf, accessed on 17 September 2022), covering the period from 6 March 2017 to 6 June 2020. The SAR images were processed using the PS-InSAR approach in the SARPROZ 6.1 (https://www.sarproz.com/sarproz-faq) software package, which has been demonstrated to be highly beneficial for InSAR data studies and has been successfully employed in previous studies on LSD [35,36,37,38,39,40,41,42,43,44], including data preparation, preliminary analysis, estimation of atmospheric phase screen (APS), and PS processing of multiple images [45]. To analyze the C-band data, at least 20 SAR images were required using the PS-InSAR method, which can monitor LS over months or years while accounting for factors such as signal noise, atmospheric conditions, and topographic influences. The sensor exhibited a ground resolution of approximately 5 m in the range direction and 20 m in the azimuth direction. It offers diverse acquisition modes, including IW, extra wide swath (EW), strip map (SM), and wave modes. Comparatively, the IW mode requires more data processing to achieve a precise co-registration of images with a high accuracy of up to 0.001 pixels. All images used in this study were obtained in the IW acquisition mode. This mode covers a single scene encompassing an area of 250 km² and is further divided into three sub-swaths using the Terrain observation by progressive scans (TOPS) mode. With a combination of high spatial and temporal resolutions and a short revisiting time, SAR images allow for the investigation of subsidence phenomena from space.

3.2. Research Methodology

The research methodology workflow depicted in Figure 3 consists of four distinct phases: data acquisition, preprocessing, analysis, and classification.

Data Preparation

The research methodology primarily focused on the data preparation steps, importing all the single look complex (SLC) data with accurate orbits, and selecting images with the same rotation of the ascending track. Initially, the images were polarized based on path information, and both slave and master images were chosen. This is a very important step, where orbit and track information were geolocated on Google Earth to remove SAR images of different tracks or different area coverages to create interferograms [46]. To conduct the analysis, we retrieved master images that covered the study region of Lahore, and subsequently acquired slave images covering the same common area as the master image (Figure 4). Subsequently, a star graph was generated, connecting the slave images to the master image for in-depth analysis.

3.3. Groundwater and PS-InSAR

3.3.1. PS-InSAR Data Analysis

The multi-temporal InSAR of PS-InSAR techniques utilizes more than 20 SAR images captured at different times to estimate the surface velocity over a year with a single master image. This technique involves selecting stable pixels as (PS points) points with stable amplitudes and phases. These points are typically dense in rock outcrops and urban areas, but are dispersed in forest and agricultural areas. In total, 128 images of the ascending-order track were used for this process. The processing of this technique was performed using SARPROZ 6.1 software, with the workflow detailed in Figure 3. To better understand the factors influencing surface displacement, this study utilized observational data from 36 groundwater wells and geological features of the study area. Additionally, population and rainfall data from the study region were used to examine the impact of groundwater depletion on the LSD. The co-registration phase, which involves aligning different images to ensure they match accurately a specific region within that specific area, is evaluated and co-registered to achieve precise geometric accuracy [47,48].

3.3.2. Preliminary Analysis

To estimate and remove the Atmosphere Phase Screen (APS) and orbit errors, the phase stability was assessed. Amplitude values are largely insensitive and cause disturbances during processing [49]. Therefore, it was anticipated that the amplitude of the pixels would be consistent across all acquisitions, and that the phase dispersion would be reduced. In SARPROZ data processing, the selection of PS points is based on the amplitude stability index (ASI).

Here, D_A is the amplitude dispersion, m_A represents amplitude the average deviation over time, and the standard deviation of amplitudes in time is given by ${\sigma }_A$.

$$\mathrm{A}\mathrm{S}\mathrm{I}=1-{\mathrm{D}}_{\mathrm{A}}=1-{\mathsf{\sigma }}_{\mathrm{A}}/{\mathrm{m}}_{\mathrm{A}}$$

(1)

3.3.3. Atmosphere Phase Screen (APS) Estimation

The estimation of APS using SAR data is a crucial step in radar remote sensing. SAR images can be affected by various atmospheric phases during data acquisition, resulting in signal interruptions. To mitigate these various disruptions and ensure data integrity, a multi-spatiotemporal filtering technique in SARPROZ is employed to calculate the APS [50]. In all InSAR methods, the fundamental observable is referred to as an interferogram, which represents the phase difference per pixel between two radar acquisitions. It is important to note that in InSAR, the displacement of each PS is measured relative to a reference PS. Once the PS points are selected using the amplitude dispersion index (ADI), the observed phase change in a given interferogram can be expressed by taking into account the influence of topography, atmospheric path delay, orbital error, and other thermal noise errors given in Equation (2).

$${\varnothing }_{int}={\varnothing }_{topography}+{\varnothing }_{movement}+{\varnothing }_{orbit}+{\varnothing }_{atmospheric}+{\varnothing }_{noise}$$

(2)

In Equation (2), ${\varnothing }_{topography}$ represents the phase variation caused by height errors, ${\varnothing }_{orbit}$ represents the error induced by phase discrepancies due to orbit estimation errors, ${\varnothing }_{atmospheric}$ represents the component generated by terrain displacement along the LOS path between two SAR acquisitions, and ${\varnothing }_{topography}$ represents the phase component resulting from variations in atmospheric phase delay. Finally, ${\varnothing }_{noise}$ refers to phase noise, encompassing thermal noise and other error components.

APS estimation is essential in PS-InSAR applications as it enables the calculation of linear deformation velocities and compensates for topographic height effects [51]. To estimate accurate APS, it is recommended to use an appropriate threshold of ASI > 0.75 as a reference for selecting the initial PS points. ASI is a derived metric that combines the amplitude coherence and amplitude information of each pixel that is often used as a thresholding criterion for selecting PS candidates, and an acceptable threshold of ASI > 0.6 was used in our case. Although this stringent parameter selection results in a limited number of PS points, it is necessary to accurately estimate the APS. After selecting the initial PS, a reference network was built using Delaunay triangulation, and the linear model was removed to derive the APS estimates from the phase residual using an inverse network approach. A reference point was established and the velocity was determined.

Temporal coherence analysis is a crucial step in ensuring the accuracy and reliability of APS-corrected SAR data for applications like land subsidence monitoring, infrastructure stability assessment, and geological hazard detection. The temporal coherence analysis of PS was performed to evaluate APS integrity, resulting in acceptable outcomes when the coherence exceeds the 0.6 threshold, as illustrated in Figure 5. During the Multi-Image Sparse Point Processing stage, the identification and selection of second-order PS points is performed using the same criteria of ASI > 0.6 to obtain denser PS points. Similar parameters and reference points were employed to remove APS during the APS estimation process. Finally, the PS points were geocoded and overlaid onto Google Earth, and only points with a coherence threshold of 0.7 or above were utilized for generating the subsidence map.

The deformation regions were then converted to geographical coordinates and combined with geological and groundwater extraction data in a geographic information system (GIS) study to validate the subsidence regions observed using PS-InSAR techniques. The results have been assessed within the geological context of Lahore and integrated with additional information layers in ArcGIS to evaluate the relationships between geological formations, groundwater extraction, and PS-InSAR-estimated subsidence. Furthermore, to identify and quantify the movement in the study area, we selected a stable point as a reference point and compared it with the motion of other points in the region. This was because the reference point served as a monitoring point for movement using PS-InSAR in the study region of Lahore.

3.3.4. Groundwater Level

To evaluate variations in GWL within the region, seasonal GWL data were gathered from monitoring wells by the Punjab Irrigation Department in Pakistan (https://irrigation.punjab.gov.pk, accessed on 23 September 2022) for the years 2003–2020, spanning both the pre- and post-monsoon seasons. The Punjab Irrigation Department of Pakistan monitors the seasonal GWL of groundwater wells in different locations in the province of Panjab, Pakistan. The data included GWL measurements taken at regular intervals over the 17 years, as well as information on the locations and characteristics of each well. However, there are many sporadic gaps in the data record; therefore, only the continuous measurements from 2003 to 2020 of both wet and dry seasons of 36 wells were carefully selected based on their appropriate records (Figure 1).

Subsequently, GWL data were employed to examine the long- and short-term consequences of groundwater variations on LSD measured by PS-InSAR techniques between 2017 and 2020, which coincides with the availability of sentinel-1 data, to examine the influence of GWL changes over time on surface displacement. We evaluate how the land surface responds to long-term variations and analyze the correlation between GWL depletion and surface displacement. Statistical methods, such as regression analysis, may be used for this purpose.

3.3.5. Spatial Interpolation

The data for the unmonitored locations were estimated using geostatistical techniques. The temporal interpolation step yielded a yearly dataset for each groundwater observation well. We spatially interpolated these annual data to generate raster grids with groundwater table values for each grid cell using the Kriging method. This resulted in a raster map for the years 2003, 2008, 2012, 2016, and 2020. We used the spatial interpolation method of ordinary Kriging in ArcGIS 10.8. Kriging is a local estimation technique offering the best linear unbiased estimator of unknown values of spatial and temporal variables [52,53,54,55,56].

Kriging is expressed as:

$${\mathit{Z}}_{\mathit{k}}^{\mathit{*}}=\sum _{\mathit{i}=\mathbf{1}}^{\mathit{n}}{\mathit{\lambda }}_{\mathit{i}}{\mathit{z}}_{\mathit{i}}$$

(3)

where $Z_k^{*}$ is an estimate by Kriging, ${\lambda }_i$ is a weight for $z_i$, and $z_i$ is a variable. The weight is determined to ensure that the estimator is unbiased and the estimation variance is minimal [57].

The unbiased condition of Kriging is:

$$\mathit{E}\left\{{\mathit{z}}_{\mathit{v}}-{\mathit{z}}_{\mathit{k}}^{\mathit{*}}\right\}=\mathbf{0}$$

(4)

where $z_v$ is an actual value and $z_k^{*}$ is an estimated value.

The sum of weights is:

$$\sum _{\mathit{i}=\mathbf{1}}^{\mathit{n}}{\mathit{\lambda }}_{\mathit{i}}=\mathbf{1.0}$$

(5)

The estimation variance of Kriging variance is:

$${\mathit{\sigma }}_{\mathit{k}}^2=\mathit{E}\left\{{\left[{\mathit{z}}_{\mathit{v}}-{\mathit{z}}_{\mathit{k}}^{\mathit{*}}\right]}^{\mathbf{2}}\right\}=\overline{\mathit{C}}\left(\mathit{V},\mathit{V}\right)+\mathit{\mu }-\sum _{\mathit{i}=\mathbf{1}}^{\mathit{n}}{\mathit{\lambda }}_{\mathit{i}}\overline{\mathit{C}}\left({\mathit{v}}_{\mathit{i}},\mathit{V}\right)$$

(6)

where $\overline{C}\left(V,V\right)$ signifies the covariance between sample variables, μ indicates the Lagrange parameter, and $\overline{C}\left(v_i,V\right)$ represents the covariance between the sample variable and the estimates. Several types of Kriging techniques have been developed to align suitable characteristics of user data, that is, ordinary Kriging is applied to stationary data, universal Kriging to nonstationary data, and co-Kriging is suitable for a group of correlated data. Ordinary Kriging was used to generate a graphical representation of GWL data.

3.4. Trend Analysis

Groundwater data from Lahore were analyzed using the Mann–Kendall (MK) test to identify spatiotemporal variations. The test determined both negative and positive trends, allowing the detection of slope values and changes in trends related to GWL depletion in the study area. The MK test was applied to account for the inhomogeneity of the time-series data, revealing spatial and temporal patterns. The calculations involved the use of mathematical equations to compute the MK statistics, V(S), and S, and the standardized test statistic Z.

$$S=\sum _{i=1}^{n-1}· \sum _{j=i+1}^{n}sign\left(x_j-x_i\right),$$

(7)

where n is the number of data points, and x_j and x_i refer to the data points at times j and i, respectively:

$$\mathit{s}\mathit{i}\mathit{g}\mathit{n}\left(\mathit{x}\mathit{j}-\mathit{x}\mathit{i}\right)=\left\{\begin{array}{c}+\mathbf{1}, \mathit{i}\mathit{f} ({\mathit{x}}_{\mathit{j}}-{\mathit{x}}_{\mathit{i}})>\mathbf{0}\\ \mathbf{0}, \mathit{i}\mathit{f} ({\mathit{x}}_{\mathit{j}}-{\mathit{x}}_{\mathit{i}})=\mathbf{0}\\ -\mathbf{1}, \mathit{i}\mathit{f} ({\mathit{x}}_{\mathit{j}}-{\mathit{x}}_{\mathit{i}})<\mathbf{10}\end{array}\right.$$

(8)

$$\mathit{V}\mathit{A}\mathit{R}\left(\mathit{S}\right)=\left[\mathit{n}\right.\left(\mathit{n}-\mathbf{1}\right)\left(\mathbf{2}\mathit{n}+\mathbf{5}\right)-\sum _{\mathit{p}=\mathbf{1}}^{\mathit{q}}{\mathit{t}}_{\mathit{p}}({\mathit{t}}_{\mathit{p}}-\mathbf{1})(\mathbf{2}{\mathit{t}}_{\mathit{p}}+\mathbf{5}\left.)\right],$$

(9)

$$\mathit{Z}=\left\{\begin{array}{c}\frac{\mathit{s}-\mathbf{1}}{\sqrt{\mathit{V}\mathit{A}\mathit{R}\left(\mathit{R}\right)}} \mathit{i}\mathit{f} \mathit{S}>\mathbf{0}\\ \mathbf{0} \mathit{i}\mathit{f} \mathit{S}=\mathbf{0}\\ \frac{\mathit{S}+\mathbf{1}}{\sqrt{\mathit{V}\mathit{A}\mathit{R}\left(\mathit{R}\right)}} \mathit{i}\mathit{f} \mathit{S}<\mathbf{0},\end{array}\right.$$

(10)

where n represents the length of time, t_p signifies the tied values for the pth value, and q denotes the number of tied values. Meanwhile x_j and x_i denote the data values in chronological order. In the present study, the statistical significance of the trend was tested using a Z-critical value > 1.96 with a significance level of 0.05.

To support the null hypothesis of no trend, the condition must be satisfied if −1.96 < $Z_{MK}$ < 1.96.

3.5. Correlation Analysis

To better understand the main factors affecting LS in Lahore, Pakistan, between 2017 and 2020, a correlation analysis between LS and GWL was performed using the PS-InSAR time series and groundwater wells data. To determine the accuracy of the PS-InSAR measurements, the subsidence values obtained from InSAR and the corresponding groundwater observations were compared. This is due to the possibility of a misalignment between the centers of the PS, point pixels, and groundwater observation benchmarks. We employed a 100 m radius of each groundwater observation, using buffer analysis in ArcGIS 10.8, to identify the average subsidence rate in the particular region. Subsequently, we calculated the mean LOS displacement values of the selected InSAR pixels within a 100 m radius, which were regarded as the InSAR measurements at the precise locations of the groundwater observations. To conduct the correlation analysis, we computed the coefficient of determination using an equation following a systematic analytical approach.

$${\mathit{R}}^{\mathbf{2}}=\frac{\mathit{n}\left(\mathit{x}\mathit{y}\right)-\left(\mathit{\Sigma }\mathit{x}\right)\left(\mathit{\Sigma }\mathit{y}\right)}{\sqrt{[\mathit{n}\ast ({\mathit{\Sigma }{\mathit{x}}^{\mathbf{2}}-\mathit{\Sigma }\mathit{x})}^{\mathbf{2}}]\ast [\mathit{n}\ast \left({\mathit{\Sigma }{\mathit{y}}^{\mathbf{2}}-\mathit{\Sigma }\mathit{y})}^{\mathbf{2}}\right]}}$$

(11)

4. Results

4.1. Spatio-Temporal Variations in Land Subsidence

4.1.1. Land Surface Subsidence Scenario

Figure 6 displays the subsidence map obtained from the ascending track orbit (track number 100) during the analysis period, and shows a sufficient number of PS points in the study area. The subsidence map was overlaid onto a base map using ArcGIS 10.5. Figure 6 shows a cluster of points representing ground movement in Lahore.

The blue color indicates comparatively stable regions in the study area of the LS scenario, with a minimal LOS displacement of between +0 and +43 mm/year. This suggests that this particular area has no significant LS, indicating a minor change in land surface without any significant downward subsidence. These areas are primarily rural settlements and farmlands on the outskirts of Lahore. The red color represents urban areas within Lahore, experiencing high subsidence ranging from −217 mm to −325 mm. The color ramp helps visualize the ground movement, with red indicating high subsidence, blue indicating stability, and yellow/light green indicating low subsidence. Overall, urban areas exhibited significant ground sinking compared to more stable rural regions.

4.1.2. Spatial Pattern of Land Surface LOS Displacement

The subsidence area depicted in Figure 7, represented by red dots, reveals notable variations in LS from place to place during the analysis period between 2017 and 2020. To comprehensively evaluate the cumulative displacement of LS in Lahore, we conducted a detailed examination of various subsidence profiles within the study area. To achieve this, we selected 12 PS points in a region highly vulnerable to LS of the study region where the values are in between −0 mm and −325 mm/year. By overlaying on a Google Earth image, we carefully placed the selected points covering the center part of urban areas and close to the existing groundwater well locations. The subsidence profiles were plotted to evaluate the cumulative displacement of the LS. We deliberately selected PS points indicated as P1 to P4, situated in the lower part of the study area from east to west, as shown in Figure 7a. Here, the recorded subsidence values were as follows: P1 (−50.26 mm), P2 (6.45 mm), P3 (−74.79 mm), and P4 (−222.34 mm).

Moreover, the subsidence trend observed in the north–southwest direction along PS points P5 (0.54 mm), P6 (−55.47) mm), P7 (−77.59 mm), and P8 (−111.55 mm) are highlighted in Figure 7b. Similarly, we plotted subsidence profiles in the northernmost part of the study area, moving from west to east for PS points P9, P10, P11, and P12. The corresponding observed subsidence trend values recorded at these points were −168.81 mm, −176.7 mm, −04.5 mm, and 5.26 mm, respectively, as shown in Figure 7c.

However, Figure 7d exhibits a clear description of the subsidence trend at each point from 2017 to 2020. These comprehensive subsidence profiles offer valuable insights into the spatial and temporal dynamics of land subsidence in the study area.

4.1.3. Mean Annual Spatial Pattern in Land Subsidence

Figure 6 (Section 4.1.1) above illustrates the cumulative subsidence within the study area from 2017 to 2020, ranging approximately from +43 mm to −325 mm. Notably the red color indicates within the figure the significant subsidence rate in the high-risk area—about −50 mm to −325 mm of subsidence was recorded during the study year. For a better understanding of the LOS displacement in the study area, the final subsidence maps for each year of the time series (2017, 2018, 2019, and 2020) have been visualized.

The subsidence graph of the time series analysis of PS-InSAR-based LS in Figure 8a represents the observed LS for the year 2017, indicating an annual mean subsidence rate of approximately −39.14 mm/year, with the recorded subsidence measurements rate of about −88 mm/year. As shown in Figure 8b, the LS rates of the same monitoring PS points from 2018 reveal notable fluctuations in readings, with a mean subsidence rate of nearly −42.705 mm/year and an overall subsidence rate of approximately −98 mm/year. Figure 8c illustrates the LS for the year 2019, with an overall subsidence observed at −95 mm/year and a mean subsidence rate of −39.76 mm/year. Meanwhile, in Figure 8d, we see an LS rate of approximately −141 mm/year, and a mean annual subsidence trend of −61.95 mm/year. These graphical representations collectively demonstrate an increase in the LS over the years in Lahore, with an overall subsidence rate of approximately −45.86 mm/year from 2017 to 2020.

However, in the map, the highest deformation rates are indicated in red in the study area. These areas of high deformation were primarily concentrated in the central region of Lahore. Furthermore, the map illustrates that, over time, the extent and intensity of deformation in the central region of Lahore gradually increased and expanded, which is an alarming situation for the sustainable development of urban cities.

4.2. Land Subsidence in Relationship with Subsurface Geology

The subsidence observed in various regions around the world is often attributed to the consolidation properties of soils. Within the study area, it has been observed that urban development has taken place over alluvium deposits, particularly in the region characterized by the Qcu-Mender deposit belt, which was once an ancient stream channel. A literature review revealed that the city’s structures, including residential and commercial buildings, have been erected on alluvium deposits and dolomites, which occupy the lowland region between the ridges. Owing to the growing numbers of settlements and unplanned construction, water channels have been blocked, leading to water stagnation on roads and streets during the rainy season. Furthermore, the blockage of drainage in a relatively broad basin, where different streams originate from the surrounding areas, causes water to percolate into the subsurface layers of the basin. The study results indicate that there is a higher occurrence of LS in the geographical regions of Lahore City, specifically the Qcu-Mender deposit belt and Qcl-Chung areas, which are characterized by soft soil layers consisting of a mixture of sand, silt, and clay. Notably, the observed subsidence rates ranged from −110 to −310 mm during the study period from 2017 to 2020 (Figure 9). This subsidence can be attributed to the saturation of the subsurface layers caused by groundwater fluctuations, which is exacerbated by the additional load exerted by infrastructure elements within the confines of the study area.

4.3. Spatio-Temporal Variations in Groundwater Level

To analyze the potential impacts of groundwater depletion on surface land and the occurrence of LS in Lahore, Pakistan, we selected a sample of 36 monitoring wells whereat we recorded the annual GWL during 2003, 2008, 2012, 2016, and 2020, and derived the temporal change of GWL, as shown in Figure 10. Our results indicate a clear distinction in the Kriging analysis of the GWL that shows variations between the central and southwest regions of the district, with greater groundwater depletion and lower GWL observed in the northwest and southwest areas. Additionally, the Mk test results exacerbate the declining trend in GWL depth in each well across the region, with approximately 0.4 m/year depletion rate during the study period. In 2003, the maximum GWL depth was approximately 9.19 to 35.22 m. In 2020, the GWL depth increased up to 12.56 to 48.15 m in the study region.

These findings suggest that the densely populated central region may require more water consumption, and that groundwater depletion has worsened progressively between 2003 and 2020, with an average GWL decrease of approximately 0.49 m/year. The annual variation in groundwater also revealed fluctuations in the GWL over this time, indicating the occurrence of groundwater depletion.

4.4. Exploring Groundwater Depletion and Its Influential Factors

In Lahore, rainfall, the Ravi River, and irrigation branch canals are potential sources of aquifer recharge. However, the river has recently remained dry, except during the monsoon season. Therefore, the major population of Lahore relies on groundwater to meet its basic needs. This has significantly impacted groundwater resources.

The climate research unit gridded observation time series (CRU TSv4.05) is a commonly used climate dataset derived from a 0.5° latitude × 0.5° longitude grid over all land areas worldwide, except for Antarctica [58]. In this study, we used meteorological datasets of annually observed gridded precipitation from the climate research unit (CRU) with a resolution of 0.5° × 0.5° to compare the annual precipitation at each GWL observation location of the 36 monitoring wells in the study area (https://crudata.uea.ac.uk/cru/data/hrg last access: 25 March 2023). Our study indicates a significant GWL depletion at an average rate of 0.49 m/year from 2003 to 2020 (Figure 11a); however, the study also found that precipitation has been increasing significantly over the last two decades. Nevertheless, the recharge rate in the study area is much lower than the rate of groundwater extraction.

The city of Lahore has undergone significant changes in its land use patterns, owing to extensive industrialization and fast-growing trends in heavy construction. Consequently, the city has experienced a reduction in the number of irrigation fields, which has led to a decrease in infiltration. This reduction played a significant, albeit indirect, role in aquifer depletion. The population of Lahore grew from approximately 3.5 million in 1980 to approximately 9 million in 2010, and population growth has been a significant factor contributing to the increased groundwater consumption (Figure 11b). The growth of urban areas in Lahore has led to an increase in water demand for domestic and commercial use, which has put a strain on the already limited groundwater resources. Unplanned urbanization and the reduction in green spaces have created impermeable surfaces that impede aquifer recharge, which is essential for replenishing groundwater. According to the World Population Review, Lahore, with an estimated population of 13.1 million as of 2021 and a high population density of 6815 people per square kilometer, has been experiencing rapid population growth (https://data.worldbank.org last access: 5 April 2023).

Lahore’s rapidly growing population has likely led to increased demand for groundwater extraction from aquifers. The challenges of overpopulation and the large extraction of groundwater can lead to a lowering of the water level, which has led to significant subsurface arable land. With the notable decline in the groundwater level of aquifers, the land above them can sink, leading to land subsidence. This can result in damage to infrastructure and risk the lives of the urban population. Simultaneously, unplanned urbanization and the reduction in green spaces have created impermeable surfaces that impede aquifer recharge, which is essential for replenishing groundwater. These issues collectively result in a range of negative impacts, including aquifer depletion and land surface deformation. To address these challenges, it is imperative to implement sustainable urban land use planning and groundwater management strategies.

4.5. Correlation Analysis between Land Subsidence and GWL

To understand the spatial relationship, the coefficient of determination was computed based on the cumulative displacement of PS-InSAR and annual average GWL values. Figure 12 shows the coefficients of determination for mean annual GWL and LS. The central part of the city, which showed a strong positive correlation between LS and GWL, was statistically significant at p < 0.05. Based on the spatial pattern of Kriging interpolation for better visualization of the correlation analysis, the northeastern part of the area showed a strong correlation with the GWL. The weakest positive correlation was observed in south and southwest Lahore.

The impact of groundwater depletion and environmental degradation is strongly felt in Lahore’s central area, where high population density and aging infrastructure exacerbate LSD. The risk to the safety of residents is increasing because of concerns about the collapse of buildings and infrastructure. The rapid growth of unregulated urban settlements has exacerbated the strain on groundwater resources and has become a major cause of concern.

5. Discussion

5.1. The Relationship between Groundwater and Land Subsidence

The city of Lahore in Pakistan is facing severe groundwater depletion [59,60] primarily because of the excessive and unwise use of water for domestic and industrial purposes [61]. In response to this issue, the Water and Sanitation Agency (WASA) authorities have started installing more tube wells to address growing public water demands [19]. However, these tube wells were installed at depths ranging from 182.88 to 213.36 m in Lahore. This excessive groundwater extraction resulted in a continued decline in the groundwater table, and threatened the underlying aquifer of Lahore.

The rate of groundwater extraction in Lahore is estimated to be approximately 1.45 million cubic meters per day [62]. According to a WASA report, the water table in Lahore has decreased by an average of 18.59 m since 1961 [19]. In a study conducted by Basharat, Muhammad, and Sultan Ahmad Rizvi (2011), water table depletion has been reported to vary across the different areas of Lahore. This study also found that the average GWL depletion in Lahore was approximately 0.62 m/year [63]. At present, the GWL in Lahore is receding at a faster rate of approximately 0.92 m/year, which is a matter of great concern for the city’s water security [8].

Before 1876, the primary water source in the region was groundwater extracted from open wells. However, currently, the major source of water supply to the public is the groundwater supply system of the WASA tube wells installed in different parts of the city. It is estimated that the total volume of groundwater extracted by the WASA is 280–290 million gallons per day, whereas a significant amount of water (approximately 150 million gallons per day) is extracted by the private sector. The amplified pumping rates resulted in a more substantial decline in water. As a result, the excessive extraction of groundwater has caused significant LS, characterized by an anomalously shaped depression cone that has formed, particularly in the high urban development areas of the city.

The regions where the LOS displacement is lower are represented by blue and green colors, indicating shallow GWL depth and the lower depletion of groundwater. Recent findings from both PS-InSAR and groundwater observations in Lahore are alarming. They reveal a substantial decline in the GWL and a significant occurrence of subsidence over the past few years. Particularly in the northern region of the city center of Lahore, consisting of densely urbanized areas for commercial, industrial and residential purposes, the groundwater level (GWL) variations are the result of a complex interplay of multiple factors; this area, the GW level depth has declined from 35.22 to 48.15 m between 2013 and 2020. Meanwhile, the subsidence rate observed between 2017 and 2020 indicated a range of −100–−325 mm predominantly observed in areas with severe groundwater depletion (Figure 13). The fact that the subsidence rate increased during this period is alarming, as it suggests that the problem worsened over time. This can have a range of negative consequences for sustainable development in the region. The groundwater depletion rate (0.49 m/year) is also a cause for concern, because it indicates a significant loss of groundwater resources in the area. This could have serious consequences for the local population, as well as for agriculture in regions that rely on groundwater. Authorities must take action to address these issues as soon as possible. This could involve implementing measures to reduce groundwater use, such as encouraging the adoption of water-efficient practices in agriculture and industry or promoting the use of alternative sources of water. It may also be necessary to closely monitor the situation to ensure that further subsidence and groundwater depletion do not occur.

5.2. Main Contribution of This Study

This study mainly investigated groundwater depletion associated with LSD and climatic factors, as well as the geological aspects and population of the Lahore region of Pakistan under different groundwater observation points. The findings of this study indicate that groundwater depletion is primarily linked to the rapid growth of urbanization and rainfall variations over time, leading to groundwater demand scenarios. However, the study highlights that LSD and LS occurrences in urban areas are primarily associated with excessive groundwater extraction and the reduction in rainwater infiltration into the aquifer.

The major influencing factors of the LSD identified in this study align with those in prior studies associated with groundwater extraction and geological factors. For example, Blasco et al. (2019) reported a subsidence rate of around −10 to 2 mm/year and −12 to −3 mm/year in the urban areas of Rome. This subsidence was observed particularly along road networks and coastal regions because of extensive urban development and groundwater exploitation [64]. Similarly, a linear correlation between LS with groundwater abstraction, along with the geological structure of Mexico City, was reported, leading to a significant land deformation rate of −20 cm year-1 between 2015 and 2016 [65]. Notably, this phenomenon has also been observed in the coastal regions of Venice and Bangkok, which recorded subsidence rates of up to 1.4 cm y⁻¹ [17] and 12 cm y⁻¹ [66], respectively, in response to groundwater extraction. Overall, various factors, such as the excessive extraction of groundwater, the rapid growth of urbanization, the amount of local rainfall, geological factors, soil conditions, and climatic conditions, could jointly influence the occurrence of ground surface deformation and further determine the extent of environmental degradation. In addition, the majority of coastal cities, such as Shanghai and Tianjin, Karachi, Los Angeles, and Dangjin, have also reported significant LSD [11,18,67,68].

Our study analyzes the LOS displacement in Lahore from 2017 to 2020, using diverse datasets, including long-term InSAR measurements and their driving factors, to comprehensively study subsidence trends over time as they correlate with groundwater level (GWL) depth, in order to provide a completed perspective of land subsidence factors in Lahore. By combining these datasets, we uncover complex relationships between geology, groundwater extraction, and land subsidence, offering new insights not explored in prior studies [26,69]. The results also indicate that densely populated and largely constructed areas of Lahore experienced a higher rate of LS than the green areas. This indicates that rainwater infiltration in the green areas was comparatively higher, allowing for the recharge of the GWL. The use of PS-InSAR techniques enables the optimization of the time-series analysis of land surface LOS displacement in a specific region, facilitating inspection and monitoring for the sustainable development of cities, and risk mitigation of existing infrastructure in urban areas, as well as the construction of new cities, such as the ongoing RUDA (Ravi Urban Development Authority) (https://ruda.gov.pk, accessed on 13 April 2023), project in Lahore. The project examines feasibility studies related to geotechnical investigations and assessing soil and groundwater conditions in the region, as well as the potential LS concerns, which can help in assessing the associated risks and possible consequences of such issues. It also considers the effects of urbanization and development on the local environment and natural resources of groundwater availability in the region, and identifies strategies to mitigate these challenges to ensure sustainable and resilient new city development.

To minimize the adverse effects of groundwater depletion on surface land, the government must establish and implement policies that mandate strict risk management practices. Additionally, policies should incorporate the effective management of groundwater recharge, such as constructing more check dams along canals and rivers. The RUDA Authority proposed the construction of small water reservoirs on the Ravi River to control monsoon floods during the development of new cities. These projects could also aid in managing the groundwater recharge process in the Lahore region.

5.3. Uncertainties and Limitations

In this study, we evaluated the relationship between GWL and LSD using PS-InSAR techniques in Lahore. Generally, geological aspects and groundwater depletion are considered to be major sources of LSD [70]. Similarly, the LS scenario in Lahore is influenced by various factors, including excessive groundwater pumping, rapid urbanization, subsurface geological conditions, and rapid population growth. Hence, in the absence of other known tectonic and geological processes, and based on the above correlations, we suggest that rising GWL depletion is the major factor resulting in LS in Lahore, which requires more effort to mitigate LSD. To this end, there is a need to analyze the most influential factors for LS, such as the rate of groundwater abstraction, rapid urbanization, surface runoff, soil conditions, and rapid population growth. As a limitation of the current study, it is recommended that multiscale research be conducted in the future to comprehensively examine LS conditions in major cities by incorporating the abovementioned factors to achieve the United Nations sustainable development goals for cities and communities. Moreover, geotechnical and aquifer characteristics can be coupled to probe the geological mechanisms behind the LS phenomena because these factors are the major sources of LS.

6. Conclusions

This study sheds valuable light on the critical issue of land subsidence in Lahore, Pakistan, and its relationship to groundwater depletion and aquifer degradation. A reliable supply of groundwater in a particular area is not only becoming an urgent need, but also a challenge for the sustainable development of cities around the world. After the advancement of remote sensing data availability, there have been many studies conducted on groundwater depletion and land surface deformation worldwide. However, in conventional practices, groundwater management and conservation strategies rarely consider aquifer recharge or land surface development in developing countries. In this study, we monitored LS from 2017 to 2020, and highlighted the ability of PS-InSAR to monitor time series subsidence in Lahore. The results indicate that groundwater and LSD varied significantly depending on many factors, including the geological conditions of the region, and spatial variations in local rainfall, water demand, groundwater pumping, surface water availability, surface runoff and climatic conditions. This paper provides a detailed explanation of the processing parameters used and the steps taken to minimize noise and errors. The study also utilized Sentinel-1 SAR data and interferometric processing to analyze surface deformations, revealing heavy LS, at an average rate of −105 mm/year. The correlation between the detected subsidence region and GWL from the selected groundwater wells suggests that groundwater overexploitation is the main factor contributing to LS. The findings indicate that, on average, GWL in the study area depleted at a rate of 0.49 m/year between 2003 and 2020. Additionally, the results indicate that areas with high land surface movement experienced significant groundwater depletion. Overall, this study emphasizes the importance of integrating GWL measurements with PS-InSAR technology to monitor the environmental impacts of groundwater overexploitation and control groundwater abstraction. This research offers a crucial step in that direction, and calls for collaboration among various stakeholders to ensure a more sustainable and secure urban future.