Current land subsidence in Jakarta: a multi-track SBAS InSAR analysis during 2017–2022 using C-band SAR data

Abstract

In Jakarta, Indonesia, rapid urbanization and intensive groundwater extraction have led to significant land subsidence, posing challenges for sustainable urban management. This study utilized Interferometry Synthetic Aperture Radar (InSAR) integrated with multi-track data analysis to generate refined time-series and velocities of 2D displacement across Greater Jakarta. This study reveals notable variation in subsidence rates across different areas, with the most significant subsidence observed in Cikarang and the coastal regions of North Jakarta, which remains linear to this day. Notably, while the 2D approach improved accuracy by up to 53% at some locations, discrepancies at others indicated that simpler descending projections might sometimes yield better results. This underscores the necessity for a nuanced application of geodetic methods based on specific site conditions to effectively monitor and manage land subsidence in complex urban environments. Our findings highlight the critical importance of integrating multiple monitoring approaches to comprehensively address both vertical and horizontal displacement.

Keywords:

• InSAR
• SBAS
• 2D displacement
• coastal
• subsidence

1. Introduction

Jakarta, a bustling metropolitan city situated on the north coast of Java Island, is characterized by its location above Quaternary alluvium soil (Cipta et al. 2018a). This geological context renders the city vulnerable to seismic amplification and rapid subsidence, exacerbated by excessive water extraction (Cipta et al. 2018b). Flood inundation in Jakarta and its satellite cities is progressively expanding annually, largely because of land subsidence and natural soil consolidation occurring not only along the coastal areas but also inland (Bott et al. 2021). The primary cause of coastal flooding in Jakarta is land subsidence, evidenced by an annual subsidence rate exceeding 10 cm per year, far surpassing the sea level rise of just a few millimeters per year (Takagi et al. 2016; Tay et al. 2022). Linear land subsidence dominates in Jakarta and Semarang, resulting in amplified consequences, particularly on infrastructure (Andreas et al. 2018). Urgent action to mitigate land subsidence is imperative to prevent escalating damage exposure and the need for more costly and challenging adaptation strategies in the future (Ward et al. 2011).

The continuous expansion of the Jakarta region, encompassing Jakarta, Bogor, Depok, Tangerang, and Bekasi Regency (Greater Jakarta), underscores the pressing necessity for comprehensive and extensive monitoring of land subsidence across broader areas. This necessity aligns with the increasing demand for groundwater among the population, coupled with changes in land cover and reports of excessive groundwater consumption in industrial regions throughout Greater Jakarta (Taftazani et al. 2022). Excessive groundwater withdrawal stemming from residential and industrial activities, alongside alterations in land cover, serves as the principal driver of land subsidence in this region (Chaussard et al. 2013; Bott et al. 2021).

Given the widespread impact of land subsidence on urban infrastructure and water management, robust and accurate methodologies for monitoring subsidence is paramount. Since around 1980, various geodetic measurements, each with unique strengths and limitations, have been utilized to discern land subsidence, particularly in Jakarta. These methodologies include leveling (Abidin et al. 2005), Global Navigation Satellite System (GNSS) (Abdullah et al. 2021), Light Detection and Ranging (LiDAR) (Park et al. 2016), microgravity (Rosid et al. 2021), and widespread usage of Interferometry Synthetic Aperture Radar (InSAR) (Ardha et al. 2021; Rafiei et al. 2022; Suhadha et al. 2023b). Although geodetic leveling boasts high accuracy in detecting subsidence, its labor-intensive and repetitive nature makes it unsuitable for ongoing, real-time monitoring of subsidence (Raspini et al. 2014). The contemporary necessity for subsidence observations across urban growth and metropolitan area expansion demands accuracy, precision over broad regions, and continuous observation. Over the past three decades, InSAR observations have proven their reliability in meeting the current requirements of subsidence observations (Crosetto et al. 2016; Pepe and Calò 2017; Suhadha et al. 2023).

The progression of InSAR capabilities is intricately linked to the development of several validated methodologies, beginning with single-pair InSAR (Hanssen 2001), Persistent Scatterer Interferometry (PSI) (Ferretti et al. 2001), and evolving to Small Baseline Subset (SBAS) (Berardino et al. 2002). Each approach computes Line-of-Sight (LOS) displacements by analyzing the phase change of each pixel (for SBAS) and point (for PSI) in SAR images across different observation epochs. LOS displacement quantifies elevation shifts along the satellite’s LOS, indicating that the observed displacements are projections within the satellite’s viewing geometry. Consequently, the comprehensive displacement, incorporating both vertical and horizontal movements, may not be fully captured (Ferretti 2014; Fuhrmann and Garthwaite 2019; Suhadha and Julzarika 2022).

Although the utilization of InSAR in Jakarta has become increasingly prevalent, the majority, if not all, of InSAR observations in the area still concentrate on LOS measurements. Table 1 outlines the evolution of the InSAR technique for monitoring land subsidence in Jakarta and its environs. Primarily, InSAR monitoring in Jakarta continues to rely on data from a single orbit, with vertical displacement projected from LOS displacement. SAR satellites observe the Earth from oblique angles, capturing images from the northeast (ascending) and northwest (descending) directions. Analyzing vertical displacement in InSAR data becomes problematic when only one orbit is utilized, as it cannot discern between vertical and horizontal deformation (Ng et al. 2012).

Table 1. Analysis of land subsidence in Jakarta and surrounding areas categorized based on input data, employing the InSAR technique, and the resulting displacement geometry.

Input data	Periods	Methods	Displacement result	Main findings	References
JERS-1 Ascending/L band	1997–2002	Three-pass DInSAR	Vertical projected from LOS displacement	The subsidence rate was about 5–10 cm from 1993 to 1995. Maximum subsidence was found in the northwestern and central eastern parts of Jakarta.	Abidin et al. (2005)
ALOS PALSAR Ascending/L band, and COSMO-SkyMed Ascending/X band	2007–2011 and 2010–2011	PSI	Vertical projected from LOS displacement	North Jakarta subsidence in the range of 15–25 cm a year.	Koudogbo et al. (2012)
ALOS PALSAR Ascending/L band	2007–2010	PSI	Vertical projected from LOS displacement	PSI measurements vary from −29 to 38 mm/year. Subsidence occurs because of groundwater extraction for industrial use, with subsidence rates ranging from 9.5 cm/year to 21.5 cm/year in coastal areas.	Ng et al. (2012), Chaussard et al. (2013)
Sentinel-1 Descending Sentinel-1 Ascending	2017 2017–2019 2017–2018	Single DInSAR Multiple DInSAR Single DInSAR	LOS displacement and vertical projected from LOS displacement	The average subsidence is −3.685 cm, whereas the highest is about −3.685 cm. The annual average of subsidence is −3.4 cm. The accuracy with DNSS was ± 6.5 cm. The north section of Jakarta has the highest ground subsidence, measuring 6 cm/year. Jakarta shows complex and varied patterns with rates of about 6 cm per year.	Cyntia and Pudja (2018), Widodo et al. (2019), Ardha et al. (2021)
Sentinel-1 Ascending and Descending	2017–2020	PSI	Vertical projected from LOS displacement	The vertical deformation between the two tracks shows a good correlation. The ascending result is more affected by the atmospheric effect.	Hakim et al. (2020)

The fusion of data from both orbits has become imperative to attain more precise vertical and high-precision horizontal measurements, particularly in the east-west (EW) direction (Samsonov and d’Oreye 2012; Fuhrmann and Garthwaite 2019; Morishita and Kobayashi 2022; Suhadha et al. 2023). However, the limitations of InSAR in accurately measuring displacements in the north-south (NS) direction owing to the satellite’s orbital proximity to the pole must be acknowledged. Previous research has primarily focused on analyzing displacement fields caused by sudden events such as earthquakes or volcanic eruptions, or averaging 3D displacement velocities over extended periods. This approach simplifies the analysis by employing only one pair of data rather than conducting a time-series analysis.

The evolution of InSAR techniques for time-series data analysis from multi-track SAR data has witnessed significant contributions from researchers striving to enhance the precision and flexibility of deformation monitoring. Ozawa and Ueda (2011) pioneered time-series data analysis by amalgamating interferograms from multiple orbit tracks, grappling with challenges associated with large datasets. Pepe et al. (2016) introduced a novel approach utilizing minimum acceleration regularization to refine LOS displacement analysis, addressing the integration of datasets with varying acquisition times. Then, Wang et al. (2012) applied Tikhonov regularization to integrate multi-track data for determining vertical velocity by ignoring the horizontal which can raise bias when observed at locations where horizontal displacement is dominant. Previous multi-track combinations of time-series data rely on a robust calculation that needs extensive computation whereas there is a challenge for computational resources and extensive SAR nowadays. Moreover, current open SAR data often correspond to different time acquisitions (Hu et al. 2014).

In this study, we integrate multi-temporal (MT)-InSAR observations from multi-track observations, using a straightforward approach, addressing the challenge of detecting vertical displacement in the extensive dataset, using a study case in Greater Jakarta. By leveraging the nearest acquisition periods of non-redundant observations, as in Wang et al. (2012), we aim to overcome the hurdles of time-series data integration and the applicability of MT-InSAR techniques through this approach. Our results reveal significant subsidence rates, particularly in the Cikarang, North Jakarta, and West Jakarta-Tangerang border areas, with rates reaching up to −10.96 cm/year. The 2D displacement approach has enhanced the measurement precision, reducing standard deviation and RMSE values compared to GNSS data. Consequently, our study underscores the importance of continuous monitoring and a multidisciplinary approach to effectively manage and mitigate land subsidence in urban environments like Greater Jakarta, ensuring sustainable urban development and infrastructure resilience.

2. Regional setting

Jakarta, the capital city of Indonesia and home to nearly 10 million permanent residents has undergone dramatic growth over the past four decades, establishing itself as one of the 28 major cities. Statistics indicate that in 2022, Jakarta was inhabited by 10,679,951 individuals, resulting in a population density of 16,084 people per square kilometer, with an annual growth rate of 0.66% (BPS DKI Jakarta 2023). Looking ahead 30 years, Greater Jakarta will firmly be categorized as a Metacity (Younger Greater Jakarta). This trajectory will further stimulate urban development, including the conversion of land use and cover from agricultural to residential and industrial, impacting the demand for groundwater and building density.

According to the Ministry of Energy and Mineral Resources, the five primary landforms dominating Greater Jakarta are alluvium, andesite, beach ridge deposit, coastal deposit, and floodplain deposits (Figure 1) (Gafoer and Samodra 1993). The formation of these landforms can be attributed to the accumulation of sedimentation from mountains located south of Jakarta. This sedimentation accumulates in the north, undergoes accretion, and forms fluviovolcanic fans at the mouths of rivers. Geologically, the north coast of Java lies on the Sunda shelf, consisting of pre-Tertiary sedimentary rocks, crystalline igneous rocks, and metamorphic rocks (Verstappen 2019). Another important aspect to note is that Jakarta is intersected by approximately 13 natural and artificial rivers, some of which serve as essential sources of public water supply. These geological features, coupled with the intensive extraction of groundwater, are believed to exacerbate land subsidence in Jakarta, affecting both shallow (less than 40 m) and deep (more than 40 m) areas (Abidin et al. 2005).

Figure 1. (a) Geological overview of Greater Jakarta. (b) Elevation profiles of each city within Greater Jakarta overlaid with available GNSS continuous operating Reference system (CORS) stations. (c) Sentinel-1 track footprint on ascending and descending orbits. (Data source: DEM – BIG; Geology – Geological Agency, Ministry of Energy and Mineral Resources, Republic of Indonesia).

Although the utilization of InSAR for monitoring land subsidence in Jakarta is well-established, its methodologies can still be improved. For instance, Hakim et al. (2020) analyzed data from two distinct orbits (ascending and descending) from 2017 to 2020 without integrating them. Despite Jakarta’s relatively gentle topography, notable variations were observed in the InSAR measurements between the two Sentinel-1 orbits. Consequently, this study aims to bridge this gap by combining InSAR images from the Sentinel-1 satellite in ascending and descending orbits to derive both vertical and horizontal displacements in a 2D displacement analysis using long time-series data. As highlighted by Ng et al. (2012), whereas Jakarta is primarily subjected to vertical deformation, the failure to segregate vertical and horizontal data may lead to reduced accuracy.

Therefore, this research endeavors to comprehensively map land subsidence in Jakarta, commencing with the border region encompassing Greater Jakarta, with precision for broader areas and high accuracy for each time series. This approach is anticipated to fulfill the requirement for complete land subsidence mapping, serving as a crucial input to mitigate the risk of tidal floods transitioning to river floods when land subsidence is the primary contributing factor.

3. Material and methods

Before detailing our methodology, contextualizing previous studies on land subsidence in Jakarta is necessary. Monitoring land subsidence using InSAR in Jakarta is not a new concept, as several previous studies have demonstrated its reliability through various methodologies and sensors. The literature review presented in Table 1 outlines the evolution of the InSAR approach, which has been utilized to monitor land subsidence in Jakarta and its surroundings. Overall, InSAR monitoring predominantly relies on single-orbit InSAR, with vertical displacement estimated from projected LOS displacement based on incidence angle. Building upon this previous work, our study leverages Sentinel-1A (S1A) data to further investigate land subsidence patterns in Jakarta using multi-track data to decompose 2D displacement. The subsequent sections provide detailed insights into the datasets and methodologies employed.

3.1 Datasets

This study utilized two tracks of S1A to investigate land subsidence. The tracks comprise the ascending and descending flight directions of the Sentinel-1 satellite. Sentinel-1 employs a C-band with a wavelength of 5.547 cm to acquire images at an oblique angle with a right-looking configuration. Two S1A datasets from each orbit, spanning from 2017 to 2022, were utilized to observe long-term subsidence.

Figure 2 illustrates the time frame of each dataset, encompassing both ascending and descending datasets. Additionally, Table 2 provides details on the acquisition periods, data volume, essential image specifications, and the mean standard error of decorrelation.

Figure 2. Temporal distribution of S1A data used in this study spanning from 2017 to 2022.

Table 2. Basic specifications of the Sentinel-1 dataset and associated uncertainties. (Data source: Sentinel-1 ascending (asc) and descending (des) data acquired from ESA; processed by the authors).

Data	Path	Frame	Orbit	Temporal span	N of data	Incidence angle (^o)	Heading angle (^o)	Mean of standard error (mm)
1	98	1160	ASC	01/01/2017–17/12/2022	170	32.205– 43.483	−11.992	1.024
2	47	609	DES	01/01/2017–13/12/2022	150	30.751–43.082	−167.999	0.094

3.2 Methodology

The methodology comprises four main steps to generate for generating and validating displacement from the SBAS multi-track (refer to Figure 3). First, the SAR data from both ascending and descending orbits were processed using Parallel SBAS (P-SBAS). Second, the resulting Line-of-Sight (LOS) measurements from the P-SBAS processing were decomposed into 2D displacement components. Next, statistical metrics such as standard deviation and Root Mean Square Error (RMSE) were used to evaluate the precision and reliability of the model. Finally, an extensive analysis was conducted using external data sources to contextualize the displacement measurements. This included incorporating land cover data to understand surface characteristics, satellite imagery to visualize changes over time, and historical flood data to correlate observed displacements with known flooding events.

Figure 3. Diagram outlining the four primary stages for analyzing vertical displacement using 2D decomposition in Greater Jakarta.

3.2.1. InSAR processing

InSAR technology leverages the phase differences between consecutive passes over the same location to calculate displacement. This phase difference is attributed to changes in proximity to the Earth’s surface (Hanssen 2001). A multi-temporal InSAR approach was adopted to address conventional D-InSAR’s limitations and acquire time-series displacement over extensive regions while mitigating decorrelation effects (Ferretti et al. 2001; Berardino et al. 2002).

The P-SBAS method, developed by the European Space Agency (ESA) and implemented within the ESA Geohazard Exploitation Platform (GEP), represents an advancement of the standard SBAS technique, specifically tailored for parallel computing environments (Casu et al. 2014). Despite its successful applications in many global geohazard studies, its implementation remains limited in Indonesia (Cigna and Tapete 2021; Suhadha and Harintaka 2024).

In this study, the P-SBAS processing method initiates the collection of data from Sentinel-1 SLC data in both ascending and descending directions, as well as SRTM digital elevation model (DEM) data and orbital information. Subsequently, the images undergo geometric registration, and pairs with high coherence, defined by a coherence threshold of 0.7, are selected to minimize decorrelation effects. The subsequent steps involve refining coregistration and generating interferograms through interferometric parameter evaluation, which is crucial for assessing phase differences. The method proceeds with temporal and spatial phase unwrapping to resolve phase ambiguities, followed by estimating any residual topography not captured by the DEM. To improve the accuracy of the results, a Ground Control Point from the runway of Soekarno–Hatta International Airport is integrated into the process. Finally, the method compiles the time series of displacement for both ascending and descending satellite data, providing comprehensive monitoring of surface movements over time.

The P-SBAS algorithm treats atmospheric correction as a subsequent step, distinct from the primary process of estimating displacement (Casu et al. 2014). Therefore, in line with the direct application of the P-SBAS algorithm, our focus was on evaluating displacement estimation within the Greater Jakarta setting.

3.2.2. 2D Displacement decomposition

In terms of procedure, InSAR observations yield displacements in a 1D geometry, specifically the LOS. The geometry is defined by the incidence angle (θ) and satellite heading (α). LOS measurements are taken in a direction orthogonal to the satellite’s heading, with most current SAR satellites utilizing the right look. The LOS displacement comprises Easting (E), Northing (N), and Upward (U) projection components (Equation (1)(1) $$v_{\mathit{\text{los}}}=\left(-\mathrm{\text{sin}}\theta \mathrm{\text{cos}}\alpha \mathrm{ }\mathrm{\text{sin}}\theta \mathrm{\text{sin}}\alpha \mathrm{ }\mathrm{\text{cos}}\theta \right)\left(\begin{array}{c}{v}_E\\ v_N\\ v_U\end{array}\right),$$(1) ) that can be utilized to convert it into 1D vertical, 2D (vertical UD and EW), and even complete 3D displacement (UD, EW, and NS) (Fialko et al. 2001): (1) $$v_{\mathit{\text{los}}}=\left(-\mathrm{\text{sin}}\theta \mathrm{\text{cos}}\alpha \mathrm{ }\mathrm{\text{sin}}\theta \mathrm{\text{sin}}\alpha \mathrm{ }\mathrm{\text{cos}}\theta \right)\left(\begin{array}{c}{v}_E\\ v_N\\ v_U\end{array}\right),$$(1)

When only one track of InSAR data is available, LOS displacement can be directly projected into a vertical direction using Equation (2)(2) $$v_U=\frac{v_{\mathit{\text{los}}}}{\mathrm{\text{cos}}\theta },$$(2) (Hanssen 2001; Fuhrmann and Garthwaite 2019). However, this approach may introduce bias as it neglects horizontal displacement, despite the primary displacement occurring vertically. Ng et al. (2012) emphasized that although vertical subsidence predominates in Jakarta, horizontal displacement near subsidence bowl rims can distort vertical displacement projections from LOS data, with a horizontal shift of just 10 mm/year potentially causing an 8 mm/year error in vertical measurement.

Moreover, if at least two passes of SAR data are available at the observation site, most recent research suggests 2D displacement modeling (Bechor and Zebker 2006; Pepe and Calò 2017). This recommendation aligns with the significant error contribution of the NS component since the orbits of most current SAR satellites are near-polar. Equation (3)(3) $$\left(\begin{array}{c}{v}_{\mathit{\text{asc}}}\\ v_{\mathit{\text{des}}}\end{array}\right)=\left(\begin{array}{c}{-\mathrm{\text{sin}}\theta }_{\mathit{\text{asc}}}{\mathrm{\text{cos}}\alpha }_{\mathit{\text{asc}}}\\ {-\mathrm{\text{sin}}\theta }_{\mathit{\text{des}}}{\mathrm{\text{cos}}\alpha }_{\mathit{\text{des}}}\end{array}\mathrm{ }\begin{array}{c}{\mathrm{ }\mathrm{\text{cos}}\theta }_{\mathit{\text{asc}}}\\ {\mathrm{ }\mathrm{\text{cos}}\theta }_{\mathit{\text{des}}1}\end{array}\right)\left(\begin{array}{c}{v}_E\\ v_U\end{array}\right),$$(3) is employed to compute vertical and EW displacements while disregarding the NS component (Fuhrmann and Garthwaite 2019; Suhadha and Harintaka 2024). The UD and EW displacements are termed 2D displacements in the subsequent section. (2) $$v_U=\frac{v_{\mathit{\text{los}}}}{\mathrm{\text{cos}}\theta },$$(2) (3) $$\left(\begin{array}{c}{v}_{\mathit{\text{asc}}}\\ v_{\mathit{\text{des}}}\end{array}\right)=\left(\begin{array}{c}{-\mathrm{\text{sin}}\theta }_{\mathit{\text{asc}}}{\mathrm{\text{cos}}\alpha }_{\mathit{\text{asc}}}\\ {-\mathrm{\text{sin}}\theta }_{\mathit{\text{des}}}{\mathrm{\text{cos}}\alpha }_{\mathit{\text{des}}}\end{array}\mathrm{ }\begin{array}{c}{\mathrm{ }\mathrm{\text{cos}}\theta }_{\mathit{\text{asc}}}\\ {\mathrm{ }\mathrm{\text{cos}}\theta }_{\mathit{\text{des}}1}\end{array}\right)\left(\begin{array}{c}{v}_E\\ v_U\end{array}\right),$$(3) where $v_{\mathit{\text{asc}}}$ and $v_{\mathit{\text{des}}}$ represent the LOS displacement in ascending and descending orbits, respectively. The incidence angle ranges from 32.205– 43.483 and 30.751–43.082 for ascending and descending, respectively. Additionally, the heading angle has a constant value both for ascending and descending (Table 2).

Equation (3)(3) $$\left(\begin{array}{c}{v}_{\mathit{\text{asc}}}\\ v_{\mathit{\text{des}}}\end{array}\right)=\left(\begin{array}{c}{-\mathrm{\text{sin}}\theta }_{\mathit{\text{asc}}}{\mathrm{\text{cos}}\alpha }_{\mathit{\text{asc}}}\\ {-\mathrm{\text{sin}}\theta }_{\mathit{\text{des}}}{\mathrm{\text{cos}}\alpha }_{\mathit{\text{des}}}\end{array}\mathrm{ }\begin{array}{c}{\mathrm{ }\mathrm{\text{cos}}\theta }_{\mathit{\text{asc}}}\\ {\mathrm{ }\mathrm{\text{cos}}\theta }_{\mathit{\text{des}}1}\end{array}\right)\left(\begin{array}{c}{v}_E\\ v_U\end{array}\right),$$(3) is applicable for calculating 2D displacements for velocity and can be solved using least square adjustment (LSA). The vertical displacement is denoted by $v_U,$ whereas the EW displacement is denoted by $v_E.$ However, challenges arise when applying this 2D formation to time-series data because of variations in acquisition time between satellite platforms (Hu et al. 2014).

This study employs a 2D construction method for time-series data, utilizing the closest acquisition period. The detailed workflow is depicted in step 2 of Figure 3. The process begins by searching for data pairs. If a pair with the same acquisition date is found, the time series can be decomposed based on these pairs. If not, the search continues for data with the smallest time gap between acquisitions, considering the revisit time of Sentinel-1 is 12 days. Initially, a 12-day interval was used for the first loop determination. However, other areas with different revisit times were also considered to determine the shortest difference in dates for the second loop. After identifying pairs with the same acquisition date, this condition can be expressed mathematically:

Let A = {a₁​, a₂​, … , a_n​} be the set of acquisition dates for the ascending dataset and D = {d₁​, d₂​,…, d_m​} be the set for the descending dataset. We aim to find the pair (a_i​, d_j) such that the absolute difference in days between a_i​ and d_j​ is minimized. We denote this minimum difference as Δdmin​.

Mathematically, the condition to find the optimal pair (a_i​, d_j​) is (4) $${\Delta d}_{\mathrm{\text{min}}}={\mathrm{\text{min}}}_{i,j}|a_i-d_j|,$$(4)

The constraint is that Δd_min​ should be as close to zero as possible, but not greater than the revisit time T of Sentinel-1 for the first loop. For the second loop, the constraint is relaxed to find the smallest Δd_min​ across the entire dataset without the T constraint.

The standard deviation (σ) of each dataset is calculated using the formula: (5) $$\sigma =\sqrt{\frac{\sum _{i-1}^N{(x_i-\overline{x})}^2}{N}}$$(5) where $x_i$ are the individual measurements, $\overline{x}$ is the mean of these measurements, and N is the number of measurements.

3.2.3. GNSS data processing

InSAR observations procedurally yield displacements in one-dimensional geometry (LOS), with the geometry being defined by the incidence angle (θ) and satellite heading (α). To assess the accuracy of the 2D displacement derived from P-SBAS, we compared it with external deformation data obtained from the GNSS Continuous Operating Reference System (GNSS CORS) monitoring provided by the Indonesian Geospatial Information Agency (BIG), considering velocity and time-series contexts. The GNSS CORS data were accessible for sites in Jakarta (CJKT), Cibitung (CBTU), Tangerang (CTGR), and Bogor (BAK1). However, the 2D displacement resulting from the fusion of ascending and descending data did not encompass the broader Bogor area, and the site BAK1 was excluded from the fused displacement.

Susilo et al. (2023) provided processed GNSS data and time-series displacement for each station, processed using MIT/Global Kalman filtering (GAMIT/GLOBK) software. We utilized the time-series displacement data of the GNSS CORS from 2017 to 2021. The processing of GNSS data described involves two main steps using the GAMIT and GLOBK software packages. Initially, the GAMIT software employs double-differencing methods to process ionosphere-free linear combination GPS phase observations, estimating daily station positions, atmospheric parameters, satellite orbits, and Earth orientation parameters. This step entails fixing satellite orbits to International GNSS Services (IGS) final orbits, applying second-order ionospheric corrections, adjusting for various geophysical effects such as ocean and earth tides, atmospheric pressure loading, and integrating data from IGS stations to align with the International Terrestrial Reference Frame 2014 (ITRF2014) reference frame.

In the second phase, the GLOBK software integrates these daily solutions with global GPS solutions, aligning them to the ITRF2014 frame by minimizing position disparities of selected reference sites through the computation of Helmert transformation parameters. This comprehensive approach enables the precise determination of GNSS station positions and the computation of daily time-series coordinates relative to the IGb14 realization of ITRF2014, ensuring high accuracy and alignment with global positioning standards.

The accuracy of the vertical component of the 2D displacement produced by P-SBAS was then assessed using the RMSE method. Table 3 presents the vertical accuracy under time-series conditions.

Table 3. Standard deviation and root mean square error (RMSE) of the P-SBAS InSAR-derived time-series vertical (UD) displacement compared to vertical measurements from GNSS CORS stations at each station from 2017–2021. (Data source: Sentinel-1 ascending and descending - ESA; GNSS CORS - BIG; processed by the authors).

Location	Standard Deviation (cm)				RMSEInSAR vs GNSS (cm)
	UD 2D	UD ASC	UD DES	UD GNSS	UD 2D	UD ASC	UD DES
CBTU	0.235	0.698	0.556	0.955	1.128	1.187	1.264
CJKT	0.619	0.315	0.942	1.154	1.538	2.065	1.295
CTGR	0.357	0.443	0.331	1.072	1.598	1.897	1.425
OVERALL	0.504	0.858	1.049	1.226	1.398	1.757	1.329

4. Result

4.1. 1D And 2D displacement

The P-SBAS was utilized to calculate LOS displacement. Figures 4(a and b) illustrate the projected vertical displacement velocity for each orbit, ranging from −9.29 cm/year to 4 cm/year. These figures demonstrate varying vertical displacement velocities across different orbits, particularly evident around Cikarang, the coastal areas of North Jakarta, and the border region between West Jakarta and Tangerang. These observations are consistent with the findings of Hakim et al. (2020), which also highlight similar differences between ascending and descending observations in the area.

Figure 4. (a–b) 1D Displacement velocities illustrating projected vertical displacement. (c) Scatter plot of 1D displacement between ascending and descending orbital acquisitions. (d) Vertical displacement and (e) EW displacement derived from 2D displacement decomposition. (Data source: ESA, processed on GEP by the authors).

The Cikarang region exhibited the most significant subsidence rates, reaching up to −10.96 cm/year, followed by West Jakarta and North Jakarta with maximum velocities of −5.58 and −8.47 cm/year, respectively. These observations are supported by the scatterplot in Figure 4(c). Using the LSA method, the 2D displacement revealed both UD and EW displacements. Figures 4(d and e) illustrate these displacements, respectively.

4.2. Main subsidence center

Jakarta’s landscape has been significantly influenced by water, with numerous areas of the city situated on floodplains vulnerable to tidal surges and river flooding. These floodplains, primarily located along the northern coast, regularly receive sediment deposits from rivers, notably originating from the Bogor area to the south of Jakarta, a phenomenon extensively documented by researchers such as Verstappen (2000).

Our analysis utilizing 2D displacement derived from InSAR has identified four key regions within Greater Jakarta experiencing substantial subsidence, as delineated in Figure 5:

Figure 5. (a) Primary subsidence centers in Greater Jakarta derived from 2D displacement, comprising at least four main locations, namely, (b) the border between west Jakarta and Tangerang cities, (c) the Western part of Northern Jakarta, (d) the Eastern part of Northern Jakarta bordering Bekasi, and (e) Cikarang, Bekasi, detailed with satellite imagery and their subsidence time series. (f–g) time-series pattern of each location’s vertical displacement over time. (Satellite imagery source: MAXAR imagery).

• West Jakarta and Tangerang City Border: Data spanning from 2017 to 2022 depict significant subsidence in this area. Figure 5(b) illustrates the annual vertical displacement time series for sites a and b, recording subsidence velocities of −5.58 and −4.95 cm/year, respectively.

• Coastal Areas of North Jakarta City: Figure 5(c) highlights sites c and d, registering displacements of −7.44 and −8.47 cm/year, respectively. A notable surge in subsidence occurred around mid-2018, followed by a deceleration by 2022, suggesting a potential stabilization in subsidence rates.

• Bekasi Regency and North Jakarta Coastal Areas: Sites e and f experienced relatively consistent linear subsidence from 2017 to 2022, with velocities of −5.28 and −7.81 cm/year, respectively. The displacement patterns here remained steady, with no significant fluctuations observed over the observed period.

• Cikarang, City of Bekasi: As depicted in Figure 5, this locality exhibited the highest rates of subsidence within the study area. Points g and h displayed average annual displacements of −10.96 and −10.91 cm/year, respectively. The data indicates a significant increase in subsidence rates beginning around mid-2018, highlighting the vulnerability of these areas to subsidence intensification.

Additionally, we analyzed the time-series patterns of each identified significant subsidence point (a-h), categorizing them into two groups as shown in Figures 5(f and g). In Figure 5(f), critical Points (c, f, g, and h) depict a consistent and substantial linear decline in subsidence up to the present, primarily located in North Jakarta and Bekasi City. By contrast, Figure 5(g) demonstrates that although these points exhibited a high linear decrease in subsidence until the end of 2021, the rate of land subsidence steepened beginning from early 2022, characterized as ‘gentle deformation’.

4.3. Validation using GNSS measurement

To evaluate the accuracy of the 2D displacement, we compared it with data from GNSS CORS to derive a time-series displacement. All GNSS stations are situated on pixels with adequate InSAR data, as evidenced by the presence of InSAR measurements around the station, comprising at least three points within a maximum distance of 45 meters. This ensures that the distance between the InSAR observation locations and GNSS falls within the spatial coverage of the P-SBAS observations, within 90 m × 90 m.

Our focus was specifically on assessing the consistency of our vertical displacement derived from InSAR measurements by comparing them with GNSS data in terms of standard deviation and RMSE—two statistical metrics indicative of precision and accuracy. Our analysis revealed that after applying a technique known as 2D displacement decomposition to the InSAR data, the quality of the vertical displacement from InSAR significantly improved.

For instance, the InSAR measurements at the CJKT station exhibited a larger discrepancy from the GNSS data compared with other locations. Table 3 provides a detailed breakdown of the RMSE values for InSAR-derived vertical displacement compared with GNSS measurements from 2020 to 2022. The results demonstrate that our refined InSAR measurements (following the application of the 2D technique) witnessed a decrease in standard deviation of 41.3% over the original ascending measurements and a 51.9% decrease over the descending ones. However, it is noteworthy that the GNSS data displayed a higher standard deviation than all InSAR measurements, suggesting that although GNSS generally offers accuracy, it can sometimes exhibit greater variability.

Regarding RMSE, which offers insights into the accuracy of our approach, the InSAR data closely aligned with the GNSS measurements, particularly at the CBTU station, where the vertical accuracy from 2D improved by approximately 43% and 52% compared to ascending and descending observations, respectively. Nonetheless, inconsistencies were observed in other stations as well as in the overall calculation.

Figure 6 illustrates a comparison of displacement observed through GNSS CORS with time-series SAR acquisitions. The UD displacement derived from the 2D decomposition fitter falls between the ascending and descending observations, with CTGR demonstrating a more consistent measurement. A significant disparity between the ascending and descending observations suggests the presence of EW displacement, consistently observed across all analyzed locations. However, in 2020, the gap between the ascending and descending observations widened to approximately −1 cm/year, accompanied by an EW displacement of around −1 cm/year. This trend persisted with a nearly −2 cm/year displacement in 2022, coupled with a continued increase in EW displacement. Furthermore, the trend lines of vertical displacement derived from GNSS have generally exhibited good agreement with the InSAR observations, with an overall RMSE of approximately 1.398 cm.

Figure 6. Comparison of displacement derived from InSAR over time with GNSS CORS stations at (a) CJKT, (b) CBTU, and (c) CTGR stations. (Data source: BIG, processed by the authors).

Our analysis shows that while the 2D decomposition method generally reduces the RMSE and standard deviation, this is not uniformly true across all stations. In particular, stations such as CJKT exhibited larger discrepancies in 2D velocities compared to direct descending projections. This variability can be attributed to local obstructions, atmospheric conditions, and the inherent limitations of the 2D decomposition method.

5. Discussion

The observed displacement patterns offer a comprehensive perspective on the complexities of land subsidence in Greater Jakarta. The implementation of a 2D displacement analysis has unveiled both vertical and lateral ground movements, providing a detailed understanding of the city’s intricate deformation patterns over time. Subsequently, several implications of the results are discussed below.

5.1. Reliability of time-series 2D approach to estimate subsidence

The ability of InSAR to observe from both ascending and descending orbits confers an advantage in capturing displacement dynamics. As demonstrated in this study, the fusion of multi-track InSAR observations addresses the challenge of discerning vertical displacement in Greater Jakarta by isolating both vertical and horizontal movements. One of the subtleties of this research lies in the pairing strategy of the ascending and descending time-series measurements.

By aligning data based on the smallest time gap, the measurements are treated as if they were acquired almost simultaneously. This approach aims to minimize temporal discrepancies and leverage the strengths of both orbital perspectives. However, pairing data with minor time gaps can introduce variability because of differences in data availability and atmospheric conditions between the two observations. For example, Sentinel-1 did not acquire images during the June period of 2020 while the ascending observation was conducted.

The validation of the 2D displacement decomposition method against GNSS CORS data has underscored the accuracy and reliability of the InSAR measurements. A notable finding from this validation, as evidenced by the data presented in Table 3, is the enhancement in measurement quality of the 2D time-series decomposition. This improvement was particularly evident in the enhanced standard deviation figures, which, in some cases, exhibited greater performance compared with the GNSS measurements.

The reduction in standard deviation for the vertical displacement of 2D compared to individual LOS measurements is due to the fusion of two independent measurements from different perspectives, which only cover pixels observed in both directions. This condition eliminates pixels with high uncertainty; for instance, pixels obstructed by high buildings in one SAR image direction will not be processed in the 2D decomposition. This data fusion inherently reduces the overall uncertainty, leading to a lower variety of vertical displacement deviations. This principle is well-supported in the literature, emphasizing the benefits of multi-track InSAR data fusion (Bechor and Zebker 2006; Pepe and Calò 2017).

For accuracy assessment, 2D InSAR measurements at the CBTU station demonstrated notable accuracy, as evidenced by lower RMSE values compared with both ascending and descending tracks. However, it is important to note that this accuracy improvement, as well as the standard deviation, was not consistent across all stations. Although CBTU, CTGR, and OVERALL showed enhanced accuracy with the 2D method, CJKT exhibited mixed results, with 2D RMSE improving over the ascending track but not the descending track. This variability can be attributed to several factors, including:

1. Obstructions at specific sites: Locations like CJKT, which had the highest standard deviation among the three stations, likely face obstructions that affect measurement reliability.

2. Atmospheric conditions: Ascending tracks are particularly susceptible to atmospheric delays that can distort measurements. This aligns with recent findings by Hakim et al., (2020).

3. The discrepancy between ASC-DES pairs: Significant differences in measurements obtained from ascending and descending passes indicate that each pass captures different aspects of ground deformation. In such cases, the 2D method offers a more balanced representation by integrating the directional biases inherent in each track.

Our findings indicate that the 2D decomposition method enhances accuracy in certain contexts, particularly in areas with significant horizontal displacement components. However, in some cases, such as in heterogeneous building locations or obstructed sites, a simple projection of ascending and descending data may provide more reliable results. This is evidenced by larger discrepancies observed at obstructed stations such as CJKT. This suggests a hybrid approach, where the choice of method is tailored to the specific conditions of each site, may be most effective.

For consistency measurements, we analyzed the temporal evolution of measurement accuracy by assessing RMSE over different periods, as depicted in Figure 7. The 2D method achieved the lowest RMSE in the overall category for 2017. However, as the temporal extent of the data increases, as observed in the periods for 2021 and the comprehensive span of 2020–2021, a rise in RMSE is evident. This increase indicates error propagation inherent in longer time-series analyses from the input dataset (single measurement). Additionally, the accuracy improvement and variation mentioned in the previous paragraph remain consistent across different observation periods, as illustrated in Figure 7.

Figure 7. Extension of Table 3 covering 2 years and a sample for a single year.

5.2. Reliability of GNSS measurements

Although GNSS measurements are typically considered the gold standard for ground positioning, acknowledging their limitations is important. Multipath errors, occurring when GNSS signals bounce off nearby structures before reaching the receiver, can distort readings (Iliev et al. 2020). Atmospheric conditions can also affect the signal’s path, potentially introducing inaccuracies (Yang et al. 2011).

GNSS data exhibits a higher standard deviation than all InSAR measurements. This condition arises because GNSS processing is limited to 30-s intervals. Moreover, the area around the GNSS pillar in Greater Jakarta is surrounded by many buildings, causing obstruction and signal multipath. For instance, the CJKT pillar contains numerous buildings inducing obstruction and signal multipath, resulting in the highest standard deviation and RMSE of InSAR, indicating less consistency of data compared to GNSS measurements. Similar to the multipath effect in GNSS measurements, InSAR observations in this location are influenced by geometrical factors, leading to a complex mix of direct reflections from the detected source (Figure 8).

Figure 8. SAR signal reflection scheme to the obstructed GNSS pillar.

Our analysis, as reflected in Table 3, underscores GNSS’s reliability in our study’s context. Despite a slightly higher standard deviation of 1.765 cm compared with InSAR measurements, the overall RMSE of GNSS data closely aligns with InSAR data trends. This highlights the general reliability of GNSS as a method for validating subsidence patterns detected by InSAR. However, this comparison also emphasizes the variability within GNSS measurements, underscoring the importance of using accurate methods to understand and interpret subsidence accurately.

5.3. Factor contributing to subsidence

The subsidence patterns in Greater Jakarta, illustrated in Figure 5, reveal a complex interplay of factors. Analyzing each sample point considering its geological condition, land use, and proximity to recent flooding events uncovers several interconnected factors influencing land subsidence.

At Point ‘a’, situated atop the volcanic Banten Tuff formation, urban sprawl exacerbates subsidence. Building pressure on this young geological feature contributes to sinking ground, whereas poroelastic behavior, influenced by groundwater extraction, plays a crucial role (Maryadi et al. 2021). Points ‘b’ and ‘c’, located on alluvial fans and based on alluvial deposits, experience subsidence because of sediment composition. Alluvial sediments, comprising sand, silt, and clay, vary in incompressibility. Fine-grained materials, like clay, compact significantly under urban development and groundwater extraction pressures. This compaction involves effective stress and soil consolidation, where reduced pore water pressure increases stress on soil grains, resulting in subsidence (Harvey et al. 2005).

The floodplain deposits at Points ‘d’, ‘f’, ‘g’ and ‘h’ are prone to subsidence because of their fine-grained composition. These sediments compact significantly when water is extracted. Terzaghi’s consolidation theory explains how sediment layers compress over time under load, a process accelerated by groundwater removal that typically supports the sediment structure (Zoysa et al. 2021).

Point ‘e’ is situated on beach ridge deposits with sandy, permeable layers, faces unique challenges. The coarser sand structure allows for rapid particle rearrangement and densification, leading to subsidence. Granular flow principles and packing efficiency help elucidate this process (Otvos, 2005).

Substantial groundwater withdrawal from detrital aquifers exacerbates subsidence across various sediment types. Reduced pore water pressure increases effective stress on sediment particles, promoting compaction. Notable studies by Galloway et al. (1999) analyze groundwater-induced subsidence in different geological contexts.

Besides geological factors, anthropogenic activities significantly contribute. Urban development dominates all points, implying intense human activity, infrastructure development, and potential groundwater extraction, evident from the predominantly urban land cover, reflecting its classification as a Metacity (Figure 9(a)). The concentration of structures and groundwater extraction in urban settings imposes immense stress on underlying soils, especially those of alluvium and floodplains, often leading to soil compaction and subsequent land subsidence. The demand for groundwater for industrial processes is notable in industrial areas, including Points ‘c’, ‘d’ and ‘e’, leading to accelerated groundwater drawdown and land sinking.

Figure 9. (a) Land cover (LC) in 2021 derived from random Forest-supervised classification and significant subsidence locations analyzed in the subsequent chapter, (b) flood occurrences in 2021 obtained from quantitative surveys and photo documentation of remaining flood inundation at significant locations. (LC data source: Planetscope imagery; photo source: field observation by the authors, taken on July 29, 2022).

Environmental determinants, notably flooding, are also crucial (Figure 9(b)). The significant flooding event in 2021, impacting areas close to Points ‘b’, ‘d’ and others, vividly demonstrates this impact. Floodwaters add weight to the soil, causing compaction, especially during a rapid recession, potentially exacerbating ground settling and sinking.

5.4. Remedial measures and recommendations

Consistently employing the 2D displacement approach in urbanized regions like Greater Jakarta is crucial because of its capability to fully capture both vertical and horizontal displacements. However, a single InSAR measurement could be used when the potential obstruction is high. Regular monitoring utilizing InSAR and GNSS techniques should be integrated with data on groundwater levels, land cover changes, flood occurrences, and other geophysical and geodetic measurements, ensuring a comprehensive assessment of subsidence trends.

The significant subsidence in Jakarta is notably driven by groundwater extraction, particularly from the upper confined aquifer system. Urgent measures are needed to regulate extraction rates and establish sustainable limits. Additionally, urban planning requires reevaluation, particularly in high-subsidence areas. Construction regulations should restrict the height and weight of new buildings to alleviate further land pressure, whereas foundation designs should be adjusted to minimize their impact on subsidence.

Community Focus and Flood Management: Recognizing the strong correlation between subsidence and flooding, enhancing flood defenses, such as constructing large-scale embankments in vulnerable areas, is imperative. Community relocation should be considered for the most vulnerable regions. Flood preparedness education and the establishment of early warning systems are vital for residents who remain. Public awareness campaigns and collaborations among stakeholders can foster a unified approach to subsidence management.

6. Conclusions

In summary, our study demonstrates the value of integrating multi-track InSAR data to derive 2D displacement measurements in Greater Jakarta. While the 2D approach generally enhances measurement precision, the variability observed at certain stations indicates that a simple projection may sometimes yield more accurate results. This highlights the need for a nuanced application of different methodologies based on specific conditions to effectively monitor and manage land subsidence. However, if the monitoring considers not only the vertical displacement but also the horizontal displacement, this approach is advised.

Our analysis highlights the highest subsidence rates in Cikarang, Bekasi, reaching up to −10.96 cm/year, with substantial displacement also observed in North Jakarta and the West Jakarta-Tangerang border areas. By employing a 2D displacement approach, we achieved enhanced measurement precision, significantly reducing standard deviation and RMSE values compared to GNSS data, particularly evident at the CBTU station. This method effectively isolates both vertical and horizontal movements, providing a comprehensive understanding of ground deformation dynamics.

Despite the overall improvement, certain locations, such as the CJKT station, exhibited variability due to factors like obstructions and atmospheric conditions, underscoring the need for further refinement and application to other conditions, like semi-urban or rural areas. The validation of our approach against GNSS measurements has demonstrated the reliability of the InSAR-derived vertical displacement. This comparative analysis underscores the importance of integrating multiple measurement methods to accurately assess subsidence.

Our findings attribute the subsidence in Jakarta to a complex interplay of geological conditions and anthropogenic activities, notably groundwater extraction and urban development. The compressibility of alluvial sediments and the impact of flooding events further exacerbate the issue, highlighting the need for comprehensive mitigation strategies.

To address these challenges, we recommend regulating groundwater extraction, reevaluating urban planning to restrict building specifications in high-subsidence areas, and enhancing flood defenses. Integrating continuous InSAR and GNSS monitoring with additional geophysical and geodetic measurements will provide a robust framework for managing subsidence trends. Our study underscores the critical importance of a multidisciplinary approach and continuous monitoring to effectively manage and mitigate land subsidence in urban environments like Greater Jakarta, ensuring sustainable urban development and infrastructure resilience.

Author contributions

Harintaka contributed to manuscript writing, review, and editing, as well as supervising the work. Argo Galih Suhadha played a pivotal role in conceptualization, data curation, methodology, validation, visualization, and original manuscript writing, review, and editing. Agung Syetiawan participated in field data investigation, formal analysis, project administration, and original draft writing. Mohammad Ardha contributed to field data investigation, data curation, and original draft writing. Anisa Rarasati contributed to data curation, formal analysis, and manuscript review and editing. All authors have reviewed and approved the final manuscript.

Acknowledgements

The authors gratefully acknowledge the European Space Agency (ESA) for providing Sentinel-1 data. Additionally, ESA and Terradue are acknowledged for their provision of P-SBAS processing on the Geohazard Exploitation Platform (GEP) under a cooperative research agreement involving Kemen-PUPR, the Asian Institute of Technology (AIT), the National Research and Innovation Agency (BRIN), and funding support from the Asian Development Bank (ADB). We also thank the Indonesian Geospatial Agency (BIG) for providing GNSS and DEMNAS data.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

Sentinel-1 imageries obtained from ESA, can be accessed via https://scihub.copernicus.eu. GNSS CORS can be accessed via http://srgi.big.go.id/. Processed data will be made available upon request.

Additional information

Funding

This work was partially supported by the National Research and Innovation Agency (BRIN) through the Seed Funding ‘Purwarupa Sistem Pengambil Keputusan Berdasarkan Analisis Citra Satelit Batch 1 2024’ numbered 1/III.6/HK/2024 and Universitas Gadjah Mada.