Land subsidence on Java Island and its contributions to relative sea level change.

Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Sci Adv . 2026 Apr 8;12(15):eaec0172. doi: 10.1126/sciadv.aec0172 Search in PMC Search in PubMed View in NLM Catalog Add to search Land subsidence on Java Island and its contributions to relative sea level change Leonard O Ohenhen Leonard O Ohenhen 1 Department of Earth Systems Science, University of California, Irvine, CA, USA. Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing - original draft, Writing - review & editing Find articles by Leonard O Ohenhen 1, * , Manoochehr Shirzaei Manoochehr Shirzaei 2 Department of Geosciences, Virginia Tech, Blacksburg, VA, USA. 3 Institute of Water, Environment and Health, United Nations University, Ontario, Canada. Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Writing - review & editing Find articles by Manoochehr Shirzaei 2, 3 , Praveen Kumar Praveen Kumar 4 Department of Earth and Planetary Sciences, Rutgers University, Piscataway, NJ, USA. 5 Rutgers Climate and Energy Institute, Rutgers University, New Brunswick, NJ, USA. Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Writing - review & editing Find articles by Praveen Kumar 4, 5 , Arif Aditiya Arif Aditiya 6 Geospatial Information Agency, Cibinong, Indonesia. Investigation, Resources, Writing - review & editing Find articles by Arif Aditiya 6 , Ashutosh Tiwari Ashutosh Tiwari 7 Texas A&M AgriLife Research Center, Corpus Christi, TX, USA. Investigation, Methodology, Writing - review & editing Find articles by Ashutosh Tiwari 7 , James L Davis James L Davis 8 Lamont-Doherty Earth Observatory of Columbia University in the City of New York, New York, NY, USA. Formal analysis, Methodology, Writing - review & editing Find articles by James L Davis 8 , Folarin Kolawole Folarin Kolawole 8 Lamont-Doherty Earth Observatory of Columbia University in the City of New York, New York, NY, USA. Formal analysis, Investigation, Supervision, Writing - review & editing Find articles by Folarin Kolawole 8 , Estelle Chaussard Estelle Chaussard 9 FM—Research Division, Norwoord, MA, USA. Investigation, Resources, Writing - review & editing Find articles by Estelle Chaussard 9 , Nitheshnirmal Sadhasivam Nitheshnirmal Sadhasivam 2 Department of Geosciences, Virginia Tech, Blacksburg, VA, USA. Investigation, Writing - review & editing Find articles by Nitheshnirmal Sadhasivam 2 , Oluwaseyi Dasho Oluwaseyi Dasho 2 Department of Geosciences, Virginia Tech, Blacksburg, VA, USA. Investigation, Writing - review & editing Find articles by Oluwaseyi Dasho 2 , Wen Zhong Wen Zhong 2 Department of Geosciences, Virginia Tech, Blacksburg, VA, USA. Investigation, Writing - review & editing Find articles by Wen Zhong 2 , Roselyn H James Roselyn H James 10 10Alytics, Leicester, UK. Investigation, Writing - review & editing Find articles by Roselyn H James 10 , Samuel Daramola Samuel Daramola 11 Department of Civil and Environmental Engineering, Virginia Tech, Blacksburg, VA, USA. Writing - review & editing Find articles by Samuel Daramola 11 , Robert J Nicholls Robert J Nicholls 12 Tyndall Centre for Climate Change Research, University of East Anglia, Norwich, UK. 13 School of Engineering, University of Southampton, Southampton, UK. Investigation, Writing - review & editing Find articles by Robert J Nicholls 12, 13 , Philip SJ Minderhoud Philip SJ Minderhoud 14 Soil Geography and Landscape Group, Wageningen University and Research, Wageningen, Netherlands. 15 Department of Civil, Environmental and Architectural Engineering, University of Padova, Padova, Italy. 16 Department of Subsurface and Groundwater Systems, Deltares Research Institute, Utrecht, Netherlands. Investigation, Writing - review & editing Find articles by Philip SJ Minderhoud 14, 15, 16 Author information Article notes Copyright and License information 1 Department of Earth Systems Science, University of California, Irvine, CA, USA. 2 Department of Geosciences, Virginia Tech, Blacksburg, VA, USA. 3 Institute of Water, Environment and Health, United Nations University, Ontario, Canada. 4 Department of Earth and Planetary Sciences, Rutgers University, Piscataway, NJ, USA. 5 Rutgers Climate and Energy Institute, Rutgers University, New Brunswick, NJ, USA. 6 Geospatial Information Agency, Cibinong, Indonesia. 7 Texas A&M AgriLife Research Center, Corpus Christi, TX, USA. 8 Lamont-Doherty Earth Observatory of Columbia University in the City of New York, New York, NY, USA. 9 FM—Research Division, Norwoord, MA, USA. 10 10Alytics, Leicester, UK. 11 Department of Civil and Environmental Engineering, Virginia Tech, Blacksburg, VA, USA. 12 Tyndall Centre for Climate Change Research, University of East Anglia, Norwich, UK. 13 School of Engineering, University of Southampton, Southampton, UK. 14 Soil Geography and Landscape Group, Wageningen University and Research, Wageningen, Netherlands. 15 Department of Civil, Environmental and Architectural Engineering, University of Padova, Padova, Italy. 16 Department of Subsurface and Groundwater Systems, Deltares Research Institute, Utrecht, Netherlands. * Corresponding author. Email: [email protected] Roles Leonard O Ohenhen : Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing - original draft, Writing - review & editing Manoochehr Shirzaei : Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Writing - review & editing Praveen Kumar : Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Writing - review & editing Arif Aditiya : Investigation, Resources, Writing - review & editing Ashutosh Tiwari : Investigation, Methodology, Writing - review & editing James L Davis : Formal analysis, Methodology, Writing - review & editing Folarin Kolawole : Formal analysis, Investigation, Supervision, Writing - review & editing Estelle Chaussard : Investigation, Resources, Writing - review & editing Nitheshnirmal Sadhasivam : Investigation, Writing - review & editing Oluwaseyi Dasho : Investigation, Writing - review & editing Wen Zhong : Investigation, Writing - review & editing Roselyn H James : Investigation, Writing - review & editing Samuel Daramola : Writing - review & editing Robert J Nicholls : Investigation, Writing - review & editing Philip SJ Minderhoud : Investigation, Writing - review & editing Received 2025 Sep 3; Accepted 2026 Feb 27; Collection date 2026 Apr 10. Copyright © 2026 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC). This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license , which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited. PMC Copyright notice PMCID: PMC13060580  PMID: 41950319 Abstract Rising sea levels and land subsidence combine to determine relative sea level (RSL) rise, which is intensifying coastal hazards. However, many densely populated regions lack the observational infrastructure to identify and quantify land subsidence contribution to RSL, hindering effective planning of responses. Here, we used satellite radar observations to generate a high-resolution assessment of land subsidence across Java Island, Indonesia, and evaluate its contribution to 21st-century RSL change. We identify widespread and temporally evolving subsidence with rates ranging from 1 to 15 cm/year in multiple coastal cities. Using machine learning spatiotemporal clustering and ancillary datasets, we attribute the dominant subsidence mechanisms to resource extraction across various geographic and geological settings. We further construct virtual tide gauges at 5-km intervals along the northern coastline, revealing that contemporary subsidence will dominate RSL budgets over the next 25 years along >75% of the coast. These findings underscore the urgent need to integrate subsidence into sea level risk and adaptation assessments in vulnerable coastal regions. Land subsidence will drive up to 85% of relative sea level rise along Java’s coast by 2050, outpacing global ocean rise. INTRODUCTION Globally, low-elevation coastal zones face an escalating array of natural hazards that threaten infrastructure, livelihoods, and ecosystems ( 1 – 3 ). These areas are subject to a confluence of acute events such as tropical cyclones, extratropical cyclones and storm surges ( 4 ), extreme rainfall ( 5 ), and wave inundation events ( 6 ); alongside chronic stressors such as shoreline erosion ( 7 , 8 ), saltwater intrusion ( 9 ), and recurrent tidal flooding ( 5 , 6 ) that progressively degrade coastal environments and communities through instantaneous, persistent, and incremental impacts ( 3 , 6 , 10 ). Together, these processes generate compound, spatially variable, and temporally evolving risks that demand both short-term emergency responses and long-term adaptation strategies that balance immediate protection and risk reduction with long-term sustainable development goals ( 11 ). Central to these evolving risks is the interaction between the rising sea surface and vertical land motion (VLM) ( 12 ). Their combined effect determines relative sea level (RSL) change ( 13 ), a key determinant of flood exposure, drainage capacity, and shoreline stability ( 3 , 10 , 14 ). While sea level rise driven primarily by climatic effects of thermal expansion, ocean mass changes associated with ice mass loss, and glacial isostatic adjustment creates spatially variable oceanic changes worldwide ( 3 , 15 ), RSL rise reflects localized conditions experienced at the coast. In many locations, land subsidence (i.e., negative VLM), often anthropogenically driven, can exceed oceanic processes by an order of magnitude ( 12 , 16 , 17 ). This distinction is critical, as land subsidence amplifies RSL rates at local scales, accelerating chronic hazards and exacerbating exposure to acute events ( 6 , 12 , 17 , 18 – 20 ). While the relative contributions of oceanic and terrestrial processes vary widely across spatial and temporal scales ( 15 ), accurately quantifying this heterogeneity requires an integrated observational framework that simultaneously accounts for oceanographic dynamics and land-surface deformation ( 3 , 21 ). However, tide gauge networks essential for measuring long-term RSL trends and capturing the combined effects of ocean processes and VLM remain sparsely distributed across global coastlines, especially in low- and middle-income nations where coastal vulnerability is often greatest due to rapid urbanization and limited adaptation capacity ( 22 , 23 ). For example, only 18 permanent service for mean sea level (PSMSL) stations from the African continent provide observational data for understanding current and projected regional sea levels in the Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report (AR6), while the 130,000 km south-eastern Asia coastline hosts just 49 stations ( 24 ). Similarly, the limited availability of Global Navigation Satellite System (GNSS) stations further constrains our ability to resolve the highly heterogeneous subsidence patterns that dominate RSL variability. These observational infrastructure gaps leave vast stretches of vulnerable coastline without sufficient long-term observational data to assess current and future local RSL dynamics. Furthermore, because VLM rates can vary substantially over short distances, even dense in situ observation networks are insufficient on their own to accurately quantify RSL. Java island, Indonesia—home to ~2% of the global population (150 million inhabitants)—represents a critical hot spot where these challenges converge, creating an archetype of compound coastal vulnerability driven by both climate and anthropogenic factors ( 25 ) ( Fig. 1 ). The Java Sea exhibits absolute sea level rise rates of 5 to 6 mm/year since 2000 ( Fig. 1 ), while urban subsidence in cities such as Jakarta reaches up to 15 cm/year ( 26 , 27 ). The resulting RSL rise rates far surpass global projections, transforming long-term risks into an immediate threat for millions of people. However, this predominant narrative of “Jakarta is sinking” ( 26 , 28 ) has obscured the broader, island-wide impacts. Existing studies have documented land subsidence in major urban centers across the island ( 26 , 27 , 29 – 31 ), particularly along the northern low-lying coastline—a densely populated corridor that is the country’s economic and population hub extending ~1500 km. With projected sea levels along this coastline potentially reaching 1 m by the end of the century ( Fig. 1, B to D ), the vulnerability of this region could be substantially magnified by subsidence. While previous studies documented land deformation patterns across Java ( 27 , 31 ), these studies have relied on single satellite orbit geometries, limiting comprehensive attribution of spatio-temporal deformation trends. This knowledge gap is further enhanced by Indonesia’s lack of representation in the PSMSL stations used for IPCC AR6 sea level projections (although local monitoring stations may exist but are not publicly accessible). Fig. 1. Study area and datasets. Open in a new tab ( A ) Population density (population/km 2 ) and synthetic aperture radar (SAR) frame boundaries across Java Island. The SAR data include 10 ascending (solid gray rectangles) and 11 descending (dashed gray rectangles) frames. Observed (2001–2024) and projected (2020–2050) absolute sea level (ASL) change for ( B ) Jakarta, West Java; ( C ) Semarang, Central Java; and ( D ) Surabaya, East Java. The locations of (B), (C), and (D) are shown in (A). The population density data is based on the 2021 WorldPop Global High Resolution Population dataset (WPGP) available at www.worldpop.org/ . The observed ASL trends (2001–2024) are from satellite altimetry dataset obtained from Copernicus Marine Environment Monitoring Service ( 79 ), while the projected ASL trends (2020–2050) are obtained from ( 24 ). The underlying basemap in (A) is provided by Esri, TomTom, Garmin, Food and Agriculture Organization (FAO), National Oceanic and Atmospheric Administration (NOAA), and U.S. Geological Survey (USGS). In this study, we address these knowledge gaps by developing a comprehensive island-wide assessment of land subsidence patterns across Java Island and evaluating its contribution to RSL rise along the northern coastline. We quantify both the spatial extent and temporal evolution of subsidence hot spots using satellite geodetic measurements, classify the dominant anthropogenic drivers, and construct historical and projected 21st century RSL trends at 5 km resolution along the northern coast. In addition, we introduce a process-based approach to interpreting drivers of land motion, shifting from static rate maps to a dynamic understanding of hazard evolution. Our findings provide critical insights for risk assessment, adaptation planning, and sustainable long-term urban development across the region, supporting more targeted and effective responses in the face of intensifying coastal risks. RESULTS Spatio-temporal variability of land subsidence on Java Island To estimate the spatial and temporal trends of land subsidence across Java, we used a multitemporal interferometric synthetic aperture radar (InSAR) approach to ~4000 SAR images from 10 ascending and 11 descending frames acquired by the Sentinel-1 C-band satellite from 2017 to 2023 ( Fig. 2 and table S1). We combined these SAR datasets with VLM rates inferred from GNSS observations at 28 sites to generate 75-m spatial resolution two-dimensional (2D) displacement fields of horizontal (east-west) (fig. S1) and VLM ( Fig. 2A ) across the island (see Materials and Methods). Fig. 2. Spatial patterns and temporal evolution of land subsidence across Java Island. Open in a new tab ( A ) VLM rates (2017–2023) for Java Island in the IGS14 geodetic reference frame. Negative VLM values indicate subsidence. ( B ) Along-coast transect of VLM sampled at 5-km intervals along the northern coastline, color-coded with acceleration, and overlaid on coastal population within 10 km (gray bars). Topography elevation is plotted in the lower panel. Positive acceleration values indicate accelerating subsidence rates, while negative values indicate deceleration. ( C ) Histogram of subsidence acceleration across Java Island, derived from Sentinel-1 observations (2017–2023). The spatial map of the acceleration is shown in fig. S3A. ( D ) Distribution of differences between Sentinel-1 (2017–2023) and ALOS-1 (2007–2011) VLM rates highlighting percentages of intensifying, stable, or slowing subsidence. The spatial map of the difference is shown in fig. S3B. The underlying basemap in (A) is provided by Esri, TomTom, Garmin, FAO, NOAA, and USGS. The horizontal land motion in a local reference frame exhibits negligible average velocity (<0.1 cm/year) across the island, reflecting overall stability (fig. S1). However, localized zones of enhanced east-west horizontal motion with peak velocities of up to 1 cm/year are observed in urban centers such as Jakarta, Bekasi Regency, Bandung, and Semarang (fig. S1, B to D). These zones of elevated horizontal displacement correlate with urban centers and areas experiencing high subsidence rates (>5 cm/year) ( Fig. 2A ) and coincide spatially with previously documented cases of structural damage to buildings, roads, and other infrastructure ( 29 , 32 ). Regions of high spatially variable horizontal displacement and subsidence rates are associated with localized crustal strain, which can induce hogging and sagging effects in structure-bending deformations that lead to cracking, tilting, and structural failure ( 33 ). The VLM data across Java indicate widespread subsidence at an average rate of −0.5 cm/year and an SD of 1 cm/year, indicating notable spatial variability across the island ( Fig. 2A ). We identify hotspots of high subsidence rates (>1 cm/year) in coastal and inland urban centers, including Jakarta, Bekasi, Bandung, Tegal, Semarang, Cilacap, Surabaya, Sidoarjo, and Jember. Other rapidly subsiding nonurban regions include Pekalongan, northern Subang, Brebes, Madiun, Demak, areas adjacent to the Lusi mud volcano, and large agricultural fields along the Bengawan Solo River in Pasi village, Lamongan Regency. The detailed analysis of VLM distributions for 12 major districts reveals that more than 20% of the urban area in the majority of these cities (9 of 12) is subsiding faster than 1 cm/year, and in Bandung, Pekalongan, and Demak, more than 30% of the city area is experiencing subsidence rates greater than 5 cm per year (fig. S2). In the majority of these 12 cities, mostly located along the northern coastline, we observe median subsidence rates of more than 0.5 cm/year and peak subsidence rates exceeding 10 cm/year (table S2). The analysis of the spatially varying VLM along the northern coastline shows that more than 25% of the 1500 km stretch of the coastline (defined as a 10-km buffer extending inland from the shoreline) is experiencing subsidence at rates greater than 1 cm/year, with the most severe sinking concentrated in low-lying elevation and highly populated regions ( Fig. 2B ). The hot spots of subsidence exposure (high subsidence rate, high population, and low elevation) along the coastline (within ~1-km inland) include Jakarta (maximum subsidence: 3.6 cm/year; population: 100,000 people; average elevation: 0.1 m), Cirebon (maximum subsidence: 3.8 cm/year; population: 12,000 people; average elevation: 0 m), Pekalongan (maximum subsidence: 10 cm/year; population: 5000 people; average elevation: −0.1 m), and Semarang (maximum subsidence: 8 cm/year; population: 10,000 people; average elevation: −0.1 m). Furthermore, we characterize the temporal evolution of VLM at each InSAR pixel as accelerating, decelerating, or monotonic (linear). For this, we applied a time-dependent quadratic regression model to the Sentinel-1 VLM time series to quantify recent nonlinear trends (see Materials and Methods; figs. S3A and S4). We then tested for the existence of a significant nonlinear pattern in the VLM data using a t test ( P < 0.05), ensuring that acceleration or deceleration trends were statistically significant (see Materials and Methods). To examine long-term changes in subsidence rates, we compared historical VLM rates from ALOS-1 L-band satellite data (2007–2010) ( 27 ) with the Sentinel-1 estimates, classifying locations where land subsidence has intensified, remained stable, or slowed relative to past trends (see Materials and Methods; fig. S3B). Our analysis reveals that 17% of currently subsiding areas exhibit statistically significant acceleration in VLM over the short timescales (2017–2023), while 45% show a deceleration in the subsidence rate ( Fig. 2C ). Along the coastline, the accelerating subsidence trend is concentrated in some high-population regions such as Cirebon (mean acceleration: 0.2 to 1.0 mm/year 2 ), Tegal (0.2 to 0.4 mm/year 2 ), and central Semarang (0.5 to 2.0 mm/year 2 ), whereas we observe linear or decelerating trends in Jakarta, Demak, Surabaya, and Sidoarjo ( Fig. 2B ). The analysis of the long-term trends reveals that ~13% of previously identified subsiding zones have transitioned from moderate subsiding conditions (<0.5 cm/year) in the late 2000s to rapid subsidence over the past decade ( Fig. 2D ). These intensifying trends are observed mostly in urban areas such as Bekasi, Tegal, Pekalongan, Semarang, Bandung, and Madiun; other rural areas including Pamanukan; and agricultural areas in Demak and Pasi, with mean subsidence rate increasing by 1 to 7 cm/year (fig. S3B). In contrast, slowed subsidence or reversed trends are observed in 22% of the land area, notably in areas of central Jakarta, Cirebon, Cilacap, Surabaya, and several rural districts, which may reflect the effects of localized subsidence control efforts such as reduced groundwater extraction. The remaining 65% of the land area exhibited no notable change in subsidence rates over the past decade, indicating persistent deformation at rates comparable to those observed in the late 2000s ( Fig. 2D ). The comparison of these long-term trends with annual flood extent (2012–2023) derived from VIIRS (Visible Infrared Imaging Radiometer Suite) observations indicates that regions experiencing a transition from stable to subsiding conditions or where subsidence has intensified such as Madiun, Pekalongan, Demak, and Pasi, also exhibit concurrent increases in the rate of flood frequency and extent (fig. S5). Dynamics and attribution of subsidence hot spots on Java Island To investigate the spatial clustering of subsidence patterns (i.e., contiguous areas with consistent deformation characteristics) and their underlying drivers across Java Island, we implemented a hybrid unsupervised machine learning framework that integrates temporal variability in VLM, spatial proximity, and dynamic deformation behavior, followed by post hoc interpretation using geospatial and geologic datasets and previous studies (see Materials and Methods; fig. S6). This approach leverages the complementary strengths of temporal pattern analysis and spatial contextualization to delineate cohesive subsidence zones with physically distinct deformation patterns. To that end, we first applied k -means temporal clustering on the horizontal and vertical displacement time series after dimensionality reduction using principal components analysis (PCA) to isolate dominant temporal trends between 2017 and 2023. Next, we refined the spatial coherence of the clusters using density-based spatial clustering of applications with noise (DBSCAN) to aggregate adjacent pixels from shared temporal clusters into contiguous spatial units based on a 1-km neighborhood radius and a subsidence threshold (>0.5 cm/year). Last, we used multivariate dynamic time warping (DTW) to compare the mean detrended cluster time series, identifying regions with analogous temporal structure (e.g., delayed onset, curvature, inflection points, or phase shifts). The k -means clustering revealed seven dominant temporal patterns, partitioned by DBSCAN into 30 spatially coherent deformation clusters ( Fig. 3A and fig. S7). These clusters capture the primary deformation hot spots across the Island, each with unique spatiotemporal deformation characteristics ( Fig. 3A and figs. S7 and S8). Fig. 3. Spatio-temporal classification and attribution of subsidence hot spots across Java Island. Open in a new tab ( A ) Spatial distribution of subsidence hotspot derived from k-means clustering of InSAR-derived VLM time series and spatial coherence refinement using DBSCAN. Each color represents a unique cluster (cluster ID). Detailed figure of each cluster is shown in fig. S9. ( B ) Conceptual model showing the dominant land use or geologic mechanism associated with each group of clusters, determined using multivariate DTW of their temporal signatures. Clusters are grouped into six classes: coastal alluvial compaction (clusters 1 to 4, 17, and 21), deltas and river channels (clusters 5, 11, and 23), industrial areas (clusters 7, 10, 12 to 14, and 24), residential areas (clusters 6, 8, 9, 15, 16, and 25 to 27), agricultural areas (clusters 18 to 20, 22, and 28), and the Lusi mud volcano (cluster 29). The underlying basemap in (A) is provided by Esri, TomTom, Garmin, FAO, NOAA, and USGS. Cluster 0, comprising areas with background VLM rates (<−0.5 cm/year) and minimal spatial coherence (<5 km), was excluded from further analysis. Broadly, we find that multiple clusters overlap spatially, likely reflecting shared land deformation driving processes across different areas. In the Western Java region, eight spatial clusters (clusters 5, 6, 7, 8, 9, 12, 13, and 15) co-occur within a 600-km 2 zone, extending from Banten-Jakarta to eastern Bekasi Regency (fig. S7B). Similarly, six clusters (clusters 1, 2, 3, 4, 17, and 21) concentrate along Central Java’s northern coastline in Cirebon, Beres, Tegal, Pekalongan, Semarang, and Demak (fig. S7A). Other multicluster regions include the Sidoarjo-Surabaya metropolitan area and Jember Regency (clusters 23, 24, and 29) and occupy the narrow southern coastal region (fig. S7D). Apart from cluster 29, cluster-level east-west motion shows a median rate of 0 cm per year while the median VLM rates vary significantly (median ranging from −1 to −7 cm/year; SD: 0.3 to 8 cm/year) (fig. S8). To classify shared deformation behavior, we grouped the 29 spatial clusters into six ( 6 ) distinct spatiotemporal groups using multivariate DTW ( Fig. 3B and fig. S9). Lower DTW distances reflect greater similarity in temporal evolution, while higher distances indicate divergent behaviors. A DTW distance threshold of 200 was selected to balance inter-group dissimilarity and intragroup cohesion. Below we discuss the deformation characteristics of the six distinct groups with inferred drivers interpreted from their temporal signatures, spatial distribution, regional geological and anthropogenic contexts, and comparisons with previous site-specific studies (table S3). Group 1 exhibits sharp, high-amplitude VLM oscillatory patterns (−1.5 to +2 cm), likely reflecting a complex mixture of spatially heterogeneous surface and subsurface processes (fig. S10). These clusters correspond to agricultural, urban, industrial, and peri-urban regions located on Central Java’s coastal alluvial deposits ( Fig. 3 and fig. S11) ( 34 ). The irregular peaks and troughs may suggest asynchronous groundwater use and recharge cycles across the group’s spatial footprint. While similar geologic deposits occur in other areas of the Island (fig. S11A), the temporal pattern of group 1 likely reflects the heterogeneous anthropogenic stressors and diverse extraction regimes ( 27 , 35 , 36 ). Group 2 is characterized by multiyear VLM oscillations (−2.5 to +1.5 cm), suggesting deformation patterns modulated by both extraction and natural hydrologic variability (fig. S10). Clusters in group 2 are found predominantly in deltas and along river channels, such as Pamanukan and the Brantas Delta in Sidoarjo, these clusters likely reflect regional climate patterns combined with human extraction responses to those conditions ( Fig. 3 and fig. S12) ( 27 ). Group 3 features moderate-amplitude oscillations (−1 to +1.3 cm) with mild seasonal variability and no extreme amplitudes, characteristic of residential zones experiencing consistent but modest land motion driven by domestic groundwater usage (fig. S13). Although the spatial and temporal characteristics of group 3 strongly suggest groundwater withdrawal as the dominant mechanism, contributions from shallow or deep sediment compaction cannot be ruled out. These clusters are distributed across urban residential zones such as central Jakarta, parts of Bandung, and Jember regency ( 26 , 27 , 29 ). Group 4 displays the least seasonal variation among all the groups (−1 to +1 cm) but exhibits a monotonic downward trend, suggesting continuous accelerating subsidence (fig. S10). These clusters are mostly located in industrial zones or mixed industrial and residential areas like Jakarta, Bekasi, Bandung, and Surabaya (fig. S12, B and D). The temporal trend of this group may reflect deformation associated with sustained, intensive deep aquifer resource (likely groundwater or hydrocarbon) extraction ( Fig. 3 and fig. S13) ( 27 , 37 ). Group 5 is characterized by large-amplitude seasonal oscillations (−3 to +2 cm), repeating temporal fluctuations representative of seasonal groundwater depletion and recharge linked to irrigation practices (fig. S10) ( 37 ). Group 5 is primarily observed in the agricultural zones along the Bengawan-Solo River plains in Pasi and agricultural fields in Madiun ( Fig. 3 and fig. S12C). Last, group 6, corresponding to the Lusi mud volcano in Sidoarjo ( Fig. 3B ), exhibits complex, episodic, nonlinear vertical (−1.5 to +1.8 cm), and horizontal (−1.5 to +3.2 cm) motion (amplitude ratio = 0.7) possibly indicative of anthropogenic signals superimposed on natural (tectonic) deformation events (fig. S10, A and B) ( 38 ). Collectively, these spatiotemporal groups reveals that land subsidence across Java Island is driven by several key mechanisms: intensive groundwater withdrawal in urban coastal centers including Jakarta, Bekasi, Semarang, and Surabaya (groups 3 and 4), seasonal agricultural extraction cycles (Group 5), natural sediment compaction and loading in deltaic environments such as Pamanukan and the Brantas Delta (groups 1 and 2), and mixed anthropogenic-tectonic processes (group 6) ( Fig. 3B ). RSL along the northern coastline of Java Island We quantified 21st century historical, present, and future (projected) RSL change along the northern coastline of Java Island by combining InSAR-derived VLM estimates with absolute sea level trends and probabilistic sea level projections to create time-resolved virtual tide gauge records at 5-km intervals along the coastline (see Materials and Methods). For historical (2001 to 2013) and present-day (2014 to 2024) RSL estimates, we integrated gridded satellite altimetry data with satellite-derived VLM rates to create spatial dense fields of observed RSL trends. For future projections (2024 to 2050), probabilistic projections of RSL were obtained using the Framework for Assessing Changes to Sea-level (FACTS) ( 39 ), presented as medium and low confidence under multiple emission scenarios [Shared Socioeconomic Pathways (SSPs) 1-1.9, 1-2.6, 2-4.5, 3-7.0, and 5-8.5). The VLM component incorporated into FACTS is assumed to be temporally stable (linear projection) over the next two decades, while we also explored potential sensitivity of the projection to nonlinear VLM. Our analysis reveals substantial spatial heterogeneity in RSL rates across the coastline ( Fig. 4 ). During the historical period, observed RSL rise rates ranged from 0.4 to 8.9 (average: 1.5) cm/year ( Fig. 4A ). Current RSL rise rates (measured over the past decade) show a notable increase compared to historical measurements, varying between 0.2 and 9.0 (average: 1.7) cm/year, with the highest rates observed in Jakarta, Cirebon, Tegal, Pekalongan, Semarang, and Demak ( Fig. 4B ). These elevated rates above historical measurements reflect increases in the land subsidence rates (5 to 100%) over the past decade despite a decrease in sea level trends over the same period (1 to 70%) (fig. S14). By 2050, under SSP2-4.5 (current emissions trajectory), projected RSL rates are expected to reach 0.4 to 9.0 (average: 1.8) cm/year. These rates represent a two- to tenfold increase compared to those attributable to sea level rise alone (fig. S14). Fig. 4. 21st century RSL rise rates along the northern coastline of Java Island. Open in a new tab Observations of RSL rise rates during ( A ) the historical (2001–2013) and ( B ) current (2014–2024) periods. ( C ) Projected RSL rise rates under shared socioeconomic pathways 2-4.5 (SSP2-4.5) scenario for the year 2050. Each circle represents a virtual tide gauge constructed at 5-km coastal segments interval. Elevation labels for major coastal cities in (B) and (C) (Jakarta: 0.2 m; Pamanukan: 0.3 m; Cirebon and Tegal: 0.3 m; Pekalongan: −0.5 m, Semarang and Demak: −0.3 m; and Surabaya: 0.4 m) indicate average elevation relative to mean sea level in meters. The underlying basemap in (A) to (C) is provided by Esri, TomTom, Garmin, FAO, NOAA, and USGS. Figure 5 presents the virtual tide gauge records for multiple areas in four representative coastal cities (Jakarta, Cirebon, Semarang, and Surabaya) (see data S1 for all locations). Across these locations, cumulative RSL shows an increase between 0.1 and 0.7 m since 2001 (start of the 21st century). The current short-term trajectories across 70% of the locations align most closely with either SSP1-2.6, representing a low-emissions pathway, or SSP2-4.5, an intermediate-emissions scenario. By 2050, RSL is projected to rise by 0.4 to 1.9 m across all scenarios (17th to 83rd percentile), with the likely range for SSP2-4.5 being 0.5 to 1.7 m, relative to the year 2001. Under low-confidence assumptions, SSP5-8.5 projections by 2050 approach 2.0 m in Jakarta, Cirebon, and Semarang. Incorporating temporally evolving VLM patterns alters projected RSL rates and trajectories under SSP2-4.5, producing differences in the rates exceeding 1 cm/year and RSL change between 0.3 and 0.7 m by 2050 (figs. S15 and S16). Fig. 5. RSL change on Java’s northern coastline. Open in a new tab 21st century measurements (gray points) and projections until 2050 (blue to red colors and shaded areas) of RSL change for ( A and B ) Jakarta, ( C and D ) Cirebon, ( E and F ) Semarang, and ( G and H ) Surabaya. Projected RSL for five Shared Socioeconomic Pathways (SSPs): SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 (medium confidence), and SSP5-8.5 (low confidence). Shaded bands represent the likely range (17th and 83rd quantile interval), and bold lines indicate the median (50th quantile) projection. Panels on the top (A, C, E, and G) represent sites with higher subsidence rates, while those on the bottom (B, D, F, and H) reflect lower-subsiding or less affected nearby sites. The time series of the observed monthly trends are shown in fig. S20. The locations of all sites are shown on Fig. 4A . To assess the relative contributions of various processes to RSL change along the coastline, we quantified the fractional contributions of climate-related (sterodynamic, glaciers, and ice sheet mass) and nonclimatic (VLM and land water storage) driven components. Now, VLM dominates the RSL budget, accounting for 25 to 91% of total RSL change under the high-emissions scenario (SSP5-8.5) (medium confidence) across 77% of the coastline ( Fig. 6A ). Over short timescales (projections by 2050), VLM remains the primary driver in more than 75% of coastal segments, contributing 27 to 85% of total RSL variability, particularly in densely populated urban centers ( Fig. 6, B and C ). Sterodynamic effects—primarily thermal expansion and regional ocean circulation shifts—and glacier loss constitute the second-largest driver but contribute less than 25% to the total RSL budget across the coastline. Ice sheet mass loss contributes minimally (<20% combined), while changes in terrestrial water storage (e.g., groundwater and reservoir fluxes) play a negligible role (<1%) in RSL variance across all time horizons and segments. Fig. 6. Component contributions to projected RSL change along Java’s northern coastline. Open in a new tab Stacked area plots of the fractional contributions of VLM, land water storage (LWS), Antarctic ice sheet (AIS), Greenland ice sheet (GIS), glaciers, and sterodynamic processes to total projected RSL change at each 5-km coastal segment for years ( A ) 2024 and ( B ) 2050 under the SSP5-8.5 scenario (medium confidence). The horizontal dashed line indicates the 0.5 threshold, below which VLM contributes more than half of the total RSL budget. ( C ) Longitudinal distribution of VLM rates, overlaid with coastal population (gray bars). DISCUSSION Java Island’s coastal plains are experiencing RSL rise at rates that exceed regional and global averages, driven by the compounding interplay of climate-driven oceanic processes and natural or anthropogenically accelerated land subsidence. The Java Sea shows current (2001 to 2024) rates of 5 to 7 mm/year, which is 25 to 75% higher than the global mean ( Fig. 1 ). However, our analysis reveals that land subsidence rates of 1 to 15 cm/year dominates present and presumably future RSL budgets along much of the densely populated northern coastline, outpacing present and predicted climate-driven sea level rise rates even under the highest-emission scenarios. These findings align with global observations that local land subsidence currently acts as the primary driver of increasing coastal vulnerability in some low-elevation regions, with South, East, and Southeast Asia being major hot spots ( 12 , 17 ). This hazard is amplifying the threat to Java Island’s coastal and inland areas far beyond Jakarta’s well-documented coastal challenges, compromising critical infrastructure, intensifying present-day chronic flooding, exacerbating 21st century flood exposure, and enhancing socioeconomic and environmental risks across the region. For Java’s coastal and inland urban areas, land subsidence also emerges independently as a severe hazard to infrastructure, driving systemic risks that extend beyond RSL. As the ground progressively sinks, critical infrastructure networks such as water supply lines, transportation routes, energy grids, and building foundations are strained to the point of failure ( 33 , 40 ). While vertical displacement directly undermines structural stability, lateral ground motion (or horizontal motion)–often underappreciated in infrastructure risk assessments–may induce differential shearing, tilting, and stress accumulation in buildings and linear infrastructure ( 33 ). We observed these dual modes of ground deformation in Jakarta, Bandung, and Semarang ( Fig. 2A and fig. S1), where damage manifests as warped roads and rail lines, ruptured pipelines, and multilevel structural failure ( 29 , 32 , 41 ). Land subsidence also subverts flood management systems, compromising the physical integrity of protective infrastructure such as levees, seawalls, and drainage networks, making flood management policies and approaches misaligned with present-day reality ( 20 , 42 ). Coastal regions that stood above high tide just a few decades ago are now routinely inundated during normal high tides, even in the absence of rainfall, while inland cities face their own variants of flood risk exacerbation as subsidence alters drainage gradients and amplifies runoff exposure (fig. S5) ( 43 ). In Semarang, for example, banjir rob (tidal flooding) has evolved from a sporadic nuisance into a chronic crisis, inundating roads and submerging homes ( 44 ). Similarly, in Bandung, subsidence along the Citarum River has reduced channel gradients, slowing drainage and trapping floodwaters in industrial and residential zones already burdened by aging infrastructure ( 29 ). These examples highlight a widespread challenge around the world, where land subsidence reinforces multiple dimensions of vulnerability, emphasizing the need for integrated and sustainable subsidence management, and forward-looking urban planning and adaptive flood control strategies in coastal regions that consider all the climate and nonclimate drivers of risk ( 11 , 42 , 45 , 46 ). In cases where subsidence is anthropogenically driven, local and city-specific groundwater policies, governance capacity, and infrastructure systems can offer effective solutions that reshape subsidence trajectories ( 16 , 45 ). On Java Island, the observed geographic distribution of VLM reveals both encouraging transitions and concerning new patterns. In parts of Jakarta, our analysis shows a current decelerating trend and a shift from rapid subsidence during 2007–2010 to near stability or even uplift in 2017–2023 in 50% of the land areas ( Fig. 2 and fig. S3A). This transition coincides with the phased implementation of groundwater extraction bans in high-risk zones, industrial relocation incentives, and expansion of piped water infrastructure to reduce aquifer dependence (fig. S13A) ( 47 ). This demonstrates that regulatory and economic measures can halt or even reverse human-induced subsidence trajectories as shown in other Asian locations ( 48 ). In addition to this policy-driven recovery, we also observe a natural uplift of the land (~2 mm/year) in the Ciliwung River delta region east of Jakarta ( Fig. 2A ), likely reflecting sediment deposition or possibly elastic unloading processes, highlighting the diverse VLM mechanisms operating across the region. Yet, these localized gains contrast with the emergence of new subsidence hot spots elsewhere on the island. As regulations tightened in core urban centers, groundwater extraction and water-intensive industries (e.g., agricultural areas) appear to have migrated to peri-urban and rural regions as observed previously around Bangkok, Thailand ( 49 ). Our analysis reveals accelerating or emerging subsidence trends in previously stable districts of Madiun, Jember, and Pasi, which previously showed minimal deformation ( Fig. 2 and fig. S3, A and B). This pattern suggests a concerning “balloon effect” and a critical gap in subsidence policy and management, where subsidence is not being eliminated but merely displaced geographically, transforming once-stable areas into zones with high rates of RSL rise, infrastructural stress, and increased socioeconomic and environmental degradation. The compounding effects of land subsidence and sea level rise create a dual hazard for many coastal communities, which poses challenges for observational monitoring and predictive modeling of coastal risk, particularly as the subsidence component is underappreciated. A major limitation lies in the sparse and uneven global distribution of observational networks (tide gauges and GNSS), restricting both the spatial resolution and coverage of RSL observations ( 21 , 22 ) . Even where measurement stations exist, they typically offer only point measurements, providing limited spatial coverage of complex coastal systems and failing to capture the heterogeneous nature of both land and oceanographic processes ( 20 ). This limitation is particularly consequential in urban centers, where human activities drive variable land motion in space and time, producing different relative sea level conditions along short coastal transects. Our analysis, along with other recent studies ( 50 , 51 ), demonstrates that satellite observations can overcome these limitations, providing a framework for creating virtual tide gauges at high spatial resolution across most coastal zones, thereby enhancing the detection of RSL hot spots that are likely to remain undetected by conventional monitoring networks. Beyond recent uncertainty quantification of nonlinearity in VLM trends ( 51 , 52 ), we show that VLM exhibits complex temporal behaviors (linear, accelerating, and decelerating phases) that are not captured by standard linear projections. This variability has direct implications for sea level projections, as extrapolations based on historical linear trends can under- or overestimate future RSL change. Our sensitivity analysis reveals that, in some locations, linear versus nonlinear projections diverge by over 40 cm by 2050, underscoring the importance of considering temporal dynamics in risk assessment and adaptation planning (figs. S15 and S16). Developing scenario-based, driver-specific subsidence projections—methodologically aligned with existing multiscenario sea level frameworks—represents an important direction for future research ( 17 , 45 ), which would equip decision-makers with a more realistic range of land-surface change futures for planning effective adaptive responses. The widespread confluence of accelerating land subsidence and sea level rise is reshaping coastal risks, demanding a paradigm shift in risk assessment and adaptation strategies. The systemic archetypes revealed in our analysis demonstrate that growth engines like resource extraction must be constrained to prevent cascading hazards across water, land, and infrastructure systems. Thus, maintaining and enhancing resilience in subsidence-prone coastal regions requires reallocating limited adaptation resources to prioritize subsidence mitigation measures along with broader coastal adaptation strategies ( 11 , 53 , 54 ). Without this integration, coastal communities risk perpetuating cycles of maladaptation, where short-term infrastructural fixes inadvertently accelerate long-term vulnerabilities by creating false impressions of security while baseline hazards continue to intensify. MATERIALS AND METHODS Land surface deformation on Java Island SAR interferometric analysis We used SAR datasets from the Sentinel-1 C-band satellite to generate surface land deformation dataset for Java Island, Indonesia. The SAR datasets include 4000 images acquired in ascending and descending orbit geometries spanning 2017 to 2023. The datasets contain 10 and 11 frames for the ascending and descending satellite orbits, respectively ( Fig. 1 and table S1). For each SAR frame, we used a multitemporal SAR interferometric approach to create the surface deformation time series in the satellite’s line-of-sight (LOS) using Wavelet-Based InSAR (WabInSAR) algorithm ( 55 , 56 ). We applied a multi-looking factor of 32 in range and 6 in azimuth directions, resulting in a pixel resolution of ~75 m. We first performed coregistration of the SAR images using precise ephemeris orbit data and Shuttle Radar Topography Mission digital elevation model and apply an enhanced spectral diversity (ESD) algorithm to minimize the interferometric phase error ( 57 ). Next, we generated more than 11,000 interferograms from the coregistered SAR images using the GAMMA software ( 58 ), using a pair selection algorithm optimized via dyadic temporal downsampling and Delaunay triangulation ( 59 ). The interferometric pairs were constrained to a maximum threshold of 150 m for perpendicular and 300 days for temporal baselines, respectively, to minimize decorrelation errors. Next, we retained so-called elite pixels with average coherence greater than 0.7 for distributed scatterers and amplitude dispersion less than 0.35 for permanent scatterers using a statistical framework ( 59 ). We then unwrapped the interferogram phases of the elite pixels using a 2D minimum cost flow algorithm ( 60 ) optimized for sparse coherent pixels ( 61 ). We corrected all unwrapped interferograms for the effect of orbital error ( 62 ) and reduced the effects of spatially uncorrelated topography error and the topography-correlated component of atmospheric delay ( 55 , 63 ). To estimate the LOS time series and rates of each pixel, we use a reweighted least squares optimization ( 55 , 56 ) using stable GPS-affiliated or stable zero-velocity pixels as local reference point. The final LOS for each SAR frame is mosaiced to generate two large-scale maps of LOS displacements for the ascending and descending datasets, following the approach detailed in ( 64 ). 2D land motion decomposition Given the satellite’s near-polar orbit geometry, which provides minimal sensitivity to north-south motion, we decompose the LOS time series into two components of land motion: horizontal (east-west) and vertical (VLM), by jointly inverting the ascending and descending datasets ( 40 ). Consequently, we first resampled the LOS displacement data from the descending track onto the location of pixels within the ascending dataset by identifying colocated pixels within a 10-m spatial radius. Let { y A , y D } and { σ A 2 , σ D 2 } be the LOS displacement and associated variances for a given pixel, respectively, where the subscripts A and D denote ascending and descending track geometries. The model for the horizontal ( d e ) and vertical ( d v ) components of displacement is expressed in Eq. 1 and solved using a weighted least squares adjustment, given in Eq. 2 [ y A y D ] = [ C A e C A v C D e C D v ] [ d e d v ] (1) X = [ G T PG ] − 1 G T PL (2) C is the unit vector for projecting the 2D displacements { d e , d v } onto the LOS ( 65 ), X is a matrix of unknowns { d e , d v }, L are the observations { y A , y D } , G is a design matrix, and P is a weight matrix, the elements of which are inversely proportional to the observant variances ( σ 2 ) . After obtaining the 2D displacement for each pixel, we performed an affine transformation to translate the vertical velocities to the IGS14 global reference frame ( 64 , 66 ). We used the vertical velocities of 10 GNSS stations (35% of the available GNSS data) provided by the Nevada Geodetic Laboratory ( 67 ) and previous studies ( 30 ). The final VLM and east-west rates are shown in Fig. 2 and fig. S1, respectively. Note that the east-west velocity is provided in the local reference frame. Error analysis , temporal variability , and validation of VLM measurements We evaluated the precision of the displacement velocities using the associated parameter uncertainties (i.e., SD), quantified the temporal variability in VLM trends, and validated the VLM rate estimates by comparison against GNSS vertical observation. We use the concept of error propagation ( 68 ) to obtain the variance-covariance matrix ( Q XX ) using Eq. 3 Q XX = [ G T PG ] − 1 (3) The SDs (formal uncertainties) for VLM and east-west velocities are shown in fig. S17. In addition to the formal uncertainties, we quantify temporal variability in the d v by estimating the posterior uncertainty associated with d v at each observation epoch. For each time step t i , the SD characterizes the uncertainty in the estimated d v , capturing deviations from a linear trend. The temporal variability for each pixel is defined using Eq. 4 σ temporal = 1 n ∑ i = 1 n σ d v ( t i ) (4) where n is the total number of time steps. This metric offers a quantification of time-dependent deformation behavior, distinguishing temporally variable motion from steady or linear displacement patterns (fig. S18). To validate the VLM rates, we compared the VLM rates in IGS14 reference frame with VLM observations from 26 available GNSS stations (including GNSS observations used for aligning the reference frame). The comparison shows a correlation of 0.83 to 0.96 and a mean and SD of 0.1 and 0.4 cm/year, respectively, for the difference between the two datasets (fig. S19). Analysis of multi-epoch VLM trends To evaluate the temporal evolution of VLM, we implemented a dual-epoch trend analysis that separately characterizes (i) recent short-term accelerations or decelerations in VLM and (ii) long-term changes in VLM trends over decadal timescales. For the short-term (2017–2023) analysis, we modeled the time-dependent displacement at each InSAR pixel using the quadratic polynomial in Eq. 5 to identify accelerating, decelerating, or monotonic (linear) deformation patterns d u ( t ) = 1 2 a t 2 + vt + ε ( t ) (5) where d u ( t ) is the cumulative vertical displacement, a is the acceleration in cm per year 2 , v is the velocity in cm per year, and ε ( t ) represents the observation error. To ensure the robustness of our analysis, we tested for the existence of a statistically significant ( P < 0.05) nonlinear pattern in each time series using a t test. If a statistically significant relationship exists, then we classify the pixels by the sign of a into accelerating (positive) or decelerating (negative) subsidence trends and assign “linear” if the test fails to reject the null hypothesis. For the long-term (multi-epoch) analysis, we compared contemporary VLM rates derived from Sentinel-1 observations (2017–2023) with historical rates derived from ALOS-1 L-band data spanning 2007–2010. The ALOS-1 dataset was obtained from Chaussard et al. ( 27 ), which applied InSAR observations to map VLM trends across Java and parts of Sumatra Islands, Indonesia. After resampling the ALOS-1 measurements to the Sentinel-1 pixels, we computed the velocity difference for each pixel using Eq. 6 ∆ v = v S − v A (6) where, v S is the VLM rate derived from Sentinel-1 observations and v A is the corresponding rate from the ALOS-1 dataset. Based on the magnitude and sign of ∆ v , we adopted the scheme in Eq. 7 to classify the subsidence trends as intensifying/emergent (IE), stable (ST), or slowing (SL) ∆ v = { IE if ∆ v < − 0.5 cm per year ST if ∣ ∆ v ∣ ≤ 0.5 cm per year SL if ∆ v > 0.5 cm per year (7) We selected these thresholds to account for systematic variations in processing methodologies between the two datasets, differences in the orbit geometries, and to accommodate the measurement uncertainties inherent in the ALOS-1 observations. This comparison enabled identification of areas where land subsidence has worsened, stabilized, or improved over the past decade (figs. S3 and S4)—reflecting both policy-driven interventions (e.g., groundwater regulations) and emergent stressors from expanding urban and resource demand (fig. S13). Classification and attribution of deformation hotspots on Java Island To identify and characterize dominant patterns of land deformation across Java Island, we developed a hybrid unsupervised machine learning framework that integrates temporal dynamics of deformation, spatial contiguity of subsiding regions, and multivariate time series comparison of the deformation trends. The objective of this classification was to delineate spatio-temporal coherent subsidence hot spots and to attribute these deformation behaviors to plausible underlying anthropogenic and geophysical/natural drivers. Information mining methods using machine learning have proven valuable for extracting contextual insights from InSAR-derived deformation data, including detection of deformation signals ( 69 ), spatial, context-driven interpolation ( 70 ), and decomposition of complex signals into physically meaningful modes ( 71 – 73 ). We implemented a four-step methodology (fig. S6): (i) PCA for reducing dimensionality in the displacement time series, (ii) k -means temporal clustering refined with DBSCAN to identify spatially contiguous deformation zones, (iii) multivariate DTW to quantify temporal similarities among clusters, and (iv) the resulting deformation trends are interpreted using ancillary datasets including geospatial and geologic information and domain-specific contextual interpretation from previous studies to identify underlying physical drivers and their attribution. First, we applied PCA to the VLM and east-west displacement time series derived from the Sentinel-1 datasets to extract their dominant temporal features. PCA transformed the original time series into an orthogonal feature space dominated by the principal modes of variability ( 74 ). For each pixel, the first two principal components explaining ≥90% of the cumulative variance were retained, preserving both long-term and short-term variability while providing a compact representation of deformation signals for subsequent clustering. Next, we applied k -means clustering to the PCA-reduced data to group pixels with similar temporal deformation profiles. The k -means algorithm partitions the dataset into disjoint K clusters by minimizing intracluster variance J = ∑ k = 1 K ∑ x i ∈ C k ‖ x i − μ k ‖ 2 (8) where x i is the 2D PCA-transformed feature vector for pixel i and μ k is the centroid of cluster C k . We determined seven optimal K based on the elbow method, yielding clusters with consistent intra-group temporal behaviors. While this step generated clusters reflecting temporal similarity, spatial contiguity was not enforced, leading to fragmented groupings. To ensure spatial continuity and remove noisy outliers, we refined temporal clusters using DBSCAN ( 75 ). The DBSCAN algorithm aggregates spatially connected pixels with similar deformation patterns into contiguous regions by connecting clusters to density-reachable neighbors. We used a spatial search radius of 5 km and a minimum subsidence threshold of 0.5 cm/year, yielding 30 spatially coherent clusters, each representing distinct temporal and spatial footprint ( Fig. 3 and figs. S7 and S8). To compare temporal deformation trends across the 30 clusters, each cluster was represented by its mean time series. The VLM component was detrended to accentuate short-term variability and structural features such as curvature, inflection points, and seasonal variations associated with anthropogenic processes (e.g., seasonal groundwater extraction), while the east-west displacement retained its original trend to preserve geophysically meaningful horizontal signals (e.g., tectonic signals). Each cluster was then represented as a multivariate time series vector in Eq. 9 X j = [ v t j , e t j ] t = 1 T (9) where v t j and e t j represents average VLM and horizontal displacement at time t for cluster j , respectively. We then quantified the intercluster similarity in deformation behavior between pairs of clusters using DTW ( 76 ). Unlike univariate DTW, multivariate DTW considers both deformation channels simultaneously, aligning sequences to minimize cumulative Euclidean distance across both vertical and horizontal dimensions. Given two clusters 1 and 2 , represented by time series X 1 , 2 = v t ( 1 , 2 ) , e t ( 1 , 2 ) t = 1 T , the DTW distance between them is defined by Eq. 10 DTW ( 1 , 2 ) = min π ∑ ( i , j ) ϵπ X i ( 1 ) − X j ( 2 ) 2 (10) where π represents the optimal warping path aligns points in time between the two sequences (1) and (2), ‖ · ‖ 2 is the Euclidean norm over the 2D deformation vector, and i and j are the time indices from the two sequences being compared. The resulting pairwise DTW distance matrix between all 29 clusters (excluding cluster 0) was subjected to agglomerative hierarchical clustering, generating a dendrogram spatiotemporal similarity (fig. S9). In the resulting dendrogram, lower DTW distances indicate higher similarity between deformation trajectories, while higher distances reflect greater dissimilarity. This suggests that clusters that merge at low DTW distances exhibit closely matched deformation dynamics, while those connected at higher thresholds show divergent patterns. A DTW distance threshold of 200 was selected empirically to define six meta-clusters or groups refined by contextual knowledge and ancillary datasets, balancing intragroup cohesion with intergroup separation. The resulting groups reflect common deformation trajectories ranging from steady, gradual trends to multi-phased behaviors even among geographically distant clusters. We interpreted the physical attribution and underlying drivers of the resulting spatiotemporal groups through post hoc integration of multiple geospatial and geologic datasets and previous studies ( 26 , 27 , 29 , 35 – 38 ). Specifically, we obtained geologic data of Indonesia to delineate lithologic information across Java Island ( 77 ) (fig. S11A). High-resolution Google Earth imagery and land cover dataset ( 78 ) were used to classify land use patterns across the island and differentiate residential, industrial, and agricultural zones (figs. S11B and S12). Last, previous studies were used to independently validate and corroborate the inferred deformation mechanisms. The framework presented here is designed for applicability in data-limited environments where direct measurements of hydrological, geological, or anthropogenic forcing are sparse. Deformation mechanisms are therefore inferred from the integration of spatial patterns, temporal signatures, geological context, and comparison with previous site-specific studies. While this approach provides valuable insights into the most plausible dominant physical processes, some deformation classes may represent mixtures of processes including groundwater withdrawal, shallow or deep sediment compaction, and tectonic contributions that cannot be fully disentangled with current datasets. This represents a common challenge in regional-scale subsidence studies and does not diminish the value of the spatiotemporal classification approach for identifying deformation hot spots and their likely dominant mechanisms. Estimation of RSL trends across Java’s northern coastline We estimated historical (2001–2013), present-day (2014–2024), and future (until 2050) RSL rise along the northern coastline of Java Island by integrating satellite-derived VLM rates with absolute sea level anomalies from satellite altimetry and probabilistic sea level rise projections. To resolve spatial variability of RSL, we constructed high-spatial resolution (5 km), time-resolved virtual tide gauge records along the coastline. To estimate the observed RSL change (historical and present-day epochs), we combined absolute sea level anomalies (SLA) from satellite altimetry data with InSAR-derived VLM estimates. The altimetry data were obtained from the Copernicus Marine Environment Monitoring Service ( 79 ), which provides gridded monthly SLA at 0.25° spatial resolution. The SLA dataset includes all standard corrections, including adjustments for instrumental drift, tidal and atmospheric effects, and geophysical processes ( 79 , 80 ). Because the SLA resolution is coarser than our 5-km coastal segment grid spacing, we selected the closest SLA grid point to each coastal segment using a minimum distance spatial proximity threshold. SLA values represent instantaneous absolute sea surface height at each time step ( i ) and were converted to RSL using the VLM estimates, applied cumulatively as a deterministic, constant-rate process across the SLA time series ( s i ) . Let v represent the annual VLM rate assigned to each 5-km coastal segment and N the number of time steps per year (i.e., N = 12 for monthly SLA dataset). The RSL at each time step ( r i ) is given by Eq. 11 r i = L v ( s i ) = s i + ( v N ) · i (11) where L v denotes a linear operator applying a cumulative displacement from VLM, distributed uniformly across subannual intervals assuming land motion progresses linearly and independently of the SLA signals. To account for differences in VLM across observational epochs, let v i represent the VLM rate assigned to time step i , applied as a piecewise-constant function such that v i = v A for i ∈ [ 1 , T 1 ] (corresponding to 2001–2013) and v i = v S for i ∈ ( T 1 , T ] (corresponding to 2014–2024). The trend of the reconstructed RSL series for each epoch satisfies Eq. 12 dr dt = ds dt + v (12) where ds dt is the trend of absolute SLA, dr dt , and is the trend of the RSL (fig. S20). Validation of virtual tide gauge records To validate the reconstructed RSL time series, we compared the virtual tide gauge estimates against observations from the PSMSL station at Cilacap. Although this station is located on the southern coast of Java and outside the primary analysis domain, the Cilacap tide gauge provides the longest continuous sea level record in the region and serves as a critical benchmark. We find strong agreement in both temporal structure and magnitude, with a high correlation in trend ( ρ = 0.9) and a rate difference of 0.1 mm/year between the virtual tide gauge and the PSMSL record (fig. S21). For projections of RSL through mid-21st century, we applied the FACTS ( 39 ), consistent with the IPCC AR6 sea level rise projections. This framework combines probabilistic projections of oceanographic processes such as thermal expansion (sterodynamics), glacier and ice sheet mass loss, and changes in land water storage with localized VLM to produce site-specific RSL change under multiple SSPs. We incorporated both the VLM rates and their temporal variability (estimated from Eq. 4 ) to ensure that local VLM trends and confidence bounds reflect nonlinear land motion in the RSL projections. Following the procedures outlined in ( 39 ), we generated full probability distributions of projected RSL for five SSPs (SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5), providing estimates of both medium confidence and low confidence (high-impact outcomes) scenarios. Seven projection workflows, corresponding to different assumptions about ice sheet dynamics and greenhouse gas forcing, were used to account for deep uncertainty in future sea level contributions ( 39 ). To combine these outputs into confidence-bound projections, we constructed probability boxes (p-boxes) in accordance with the AR6 guidance on uncertainty quantification. For medium-confidence RSL projections, we used workflow 1e (based on the emulandice emulator for Antarctic ice sheet response) and workflow 1f for RSL rates; for low-confidence projections, we used workflows 2e and 2f, respectively. Each coastal segment’s projected RSL reflects the cumulative contributions from both oceanographic processes and contemporary (Sentinel-1 derived) land subsidence. We emphasize that in these projections, VLM is assumed to be temporally stable (linear rates), consistent with the current AR6 framework. To investigate the sensitivity of RSL projections to nonlinear land motion, we performed an auxiliary analysis to incorporate VLM accelerations (observed in some parts of the coastline) into the no-VLM p-box output through an additive adjustment sampled from empirical probability distributions (figs. S15 and S16). This highlights the influence of accelerating (or decelerating) land deformation on near-term coastal hazard projections and highlights the limitations of assuming linear subsidence rates in the context of evolving anthropogenic pressures. While the virtual tide gauge construction provides a valid first-order assessment of RSL change, it relies on several simplified assumptions that introduce limitations. First, the assumption of VLM as a constant, deterministic rate neglects temporal variability, including short-term nonlinear trends driven by anthropogenic processes. For instance, observed accelerations in VLM of 0.1 to 0.2 cm/year 2 in Pekalongan and Semarang indicate that linear rate assumptions may underestimate past and RSL rise in some subsiding regions (see Fig. 2 and figs. S3 and S4). Second, the use of satellite altimetry data—originally optimized for open-ocean conditions—to characterize coastal sea level anomalies introduces additional uncertainty. Nearshore, sea level is influenced by localized processes such as shelf circulation, riverine input, and tidal amplification, which may decouple coastal sea level behavior from offshore altimetric signals. Third, in projecting future RSL, the extrapolation of contemporary VLM rates assumes stationarity in subsidence processes, which in practice is unlikely to hold over multidecadal timescales due to shifting land use, hydrological conditions, or mitigation policies (see Discussion section for scenario-based land subsidence alternative). Nevertheless, in regions where long-term tide gauge observations are sparse or absent, the virtual tide gauge framework offers a practical and scalable alternative for generating spatially explicit RSL estimates, especially over the next few decades. Acknowledgments We thank M. Colet for helpful discussions and valuable insights that contributed to this research. We thank M. Silva for help with Fig. 3B . Funding: L.O.O. was supported by LDEO postdoctoral research fellowship. M.S., N.S., O.D., and W.Z. were supported by funding from the Department of Defense. Author contributions: Conceptualization: L.O.O. and M.S. Data curation: L.O.O., M.S., and P.K. Resources: L.O.O., M.S., A.A., P.K., and E.C. Methodology: L.O.O., M.S., P.K., A.T., and J.L.D. Formal analysis: L.O.O., M.S., P.K., J.L.D., and F.K. Software: L.O.O., M.S., and P.K. Investigation: L.O.O., M.S., P.K., A.A., A.T., F.K., E.C., N.S., O.D., W.Z., R.H.J., R.J.N., and P.S.J.M. Visualization: L.O.O. Funding acquisition: M.S. Validation: L.O.O., M.S., and P.K. Supervision: L.O.O., M.S., and F.K. Project administration: L.O.O. and M.S. Writing—original draft: L.O.O. Writing—review and editing: L.O.O., M.S., P.K., A.A., A.T., J.L.D., F.K., E.C., N.S., O.D., W.Z., R.H.J., S.D., R.J.N., and P.S.J.M. Competing interests: The authors declare that they have no competing interests. Data, code, and materials availability: The VLM and horizontal land motion data for Java Island generated for this study are available at https://doi.org/10.5281/zenodo.15786356 . The medium and low confidence projected RSL change and rates for five SSPs (SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) generated for this study are available at https://doi.org/10.5281/zenodo.15786356 . Sentinel-1 data used in this study are publicly available through the Alaska Satellite Facility and can be accessed at https://search.asf.alaska.edu . The satellite altimetry data for sea level rise are available from Copernicus Marine Environment Monitoring Service (CMEMS) and are available at https://data.marine.copernicus.eu/product/SEALEVEL_GLO_PHY_L4_MY_008_047/download . Access requires the creation of a free CMEMS user account, which is available to all users at no cost upon registration. The population for Java Island was estimated using the WorldPop dataset available at https://hub.worldpop.org/geodata/summary?id=77690 . The coastal elevation referenced to mean sea level uses the DeltaDTM dataset v1.1 available at https://doi.org/10.4121/21997565 . All other data and code needed to evaluate and reproduce the results in the paper are present in the paper and/or the Supplementary Materials. This study did not generate new materials. Supplementary Materials The PDF file includes: Figs. S1 to S21 Legends for table S1 to S3 Legend for data S1 References sciadv.aec0172_sm.pdf (3.2MB, pdf) Other Supplementary Material for this manuscript includes the following: Tables S1 to S3 sciadv.aec0172_tables_s1_to_s3.zip (38.3KB, zip) REFERENCES 1. Nicholls R. J., Cazenave A., Sea-level rise and its impact on coastal zones. Science 328, 1517–1520 (2010). [ DOI ] [ PubMed ] [ Google Scholar ] 2. Hinkel J., Lincke D., Vafeidis A. T., Perrette M., Nicholls R. J., Tol R. S., Levermann A., Coastal flood damage and adaptation costs under 21st century sea-level rise. Proc. Natl. Acad. Sci. U.S.A. 111, 3292–3297 (2014). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 3. M. Oppenheimer, B. C. Glavovic, J. Hinkel, R. van de Wal, A. K. Magnan, A. Abd-Elgawad, R. Cai, M. Cifuentes-Jara, R. M. DeConto, T. Ghosh, J. Hay, F. Isla, B. Marzeion, B. Meyssignac, Z. Sebesvari, Sea level rise and implications for low-lying islands, coasts and communities, in IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (Cambridge Univ. Press, 2019). [ Google Scholar ] 4. Dullaart J. C. M., Muis S., Bloemendaal N., Chertova M. V., Couasnon A., Aerts J. C. J. H., Accounting for tropical cyclones more than doubles the global population exposed to low-probability coastal flooding. Commun. Earth Environ. 2, 204 (2021). [ Google Scholar ] 5. Rahimi R., Tavakol-Davani H., Graves C., Gomez A., Valipour M. F., Compound inundation impacts of coastal climate change: Sea-level rise, groundwater rise and coastal precipitation. Water 12, 2776 (2020). [ Google Scholar ] 6. Barnard P. L., Befus K. M., Danielson J. J., Engelstad A. C., Erikson L. H., Foxgrover A. C., Hayden M. K., Hoover D. J., Leijnse T. W. B., Massey C., McCall R., Nadal-Caraballo N. C., Nederhoff K., O’Neill A. C., Parker K. A., Shirzaei M., Ohenhen L. O., Swarzenski P. W., Thomas J. A., van Ormondt M., Vitousek S., Vos K., Wood N. J., Jones J. M., Jones J. L., Projections of multiple climate-related coastal hazards for the US Southeast Atlantic. Nat. Clim. Chang. 15, 101–109 (2025). [ Google Scholar ] 7. Hinkel J., Nicholls R. J., Tol R. S. J., Wang Z. B., Hamilton J. M., Boot G., Vafeidis A. T., McFadden L., Ganopolski A., Klein R. J. T., A global analysis of erosion of sandy beaches and sea-level rise: An application of DIVA. Glob. Planet. Change 111, 150–158 (2013). [ Google Scholar ] 8. Vousdoukas M. I., Ranasinghe R., Mentaschi L., Plomaritis T. A., Athanasiou P., Luijendijk A., Feyen L., Sandy coastlines under threat of erosion. Nat. Clim. Chang. 10, 260–263 (2020). [ Google Scholar ] 9. Befus K. M., Barnard P. L., Hoover D. J., Finzi Hart J. A., Voss C. I., Increasing threat of coastal groundwater hazards from sea-level rise in California. Nat. Clim. Chang. 10, 946–952 (2020). [ Google Scholar ] 10. Moftakhari H. R., AghaKouchak A., Sanders B. F., Feldman D. L., Sweet W., Matthew R. A., Luke A., Increased nuisance flooding along the coasts of the United States due to sea level rise: Past and future. Geophys. Res. Lett. 42, 9846–9852 (2015). [ Google Scholar ] 11. Nicholls R. J., Adapting to sea-level rise. Resilience 13, 13–29 (2018). [ Google Scholar ] 12. Nicholls R. J., Lincke D., Hinkel J., Brown S., Vafeidis A. T., Meyssignac B., Hanson S. E., Merkens J. L., Fang J., A global analysis of subsidence, relative sea-level change and coastal flood exposure. Nat. Clim. Chang. 11, 338–342 (2021). [ Google Scholar ] 13. Gregory J. M., Griffies S. M., Hughes C. W., Lowe J. A., Church J. A., Fukimori I., Gomez N., Kopp R. E., Landerer F., Cozannet G. L., Ponte R. M., Concepts and terminology for sea level: Mean, variability and change, both local and global. Surv. Geophys. 40, 1251–1289 (2019). [ Google Scholar ] 14. FitzGerald D. M., Fenster M. S., Argow B. A., Buynevich I. V., Coastal impacts due to sea-level rise. Annu. Rev. Earth Planet. Sci. 36, 601–647 (2008). [ Google Scholar ] 15. Frederikse T., Landerer F., Caron L., Adhikari S., Parkes D., Humphrey V. W., Dangendorf S., Hogarth P., Zanna L., Cheng L., Wu Y. H., The causes of sea-level rise since 1900. Nature 584, 393–397 (2020). [ DOI ] [ PubMed ] [ Google Scholar ] 16. Erkens G., Bucx T., Dam R., De Lange G., Lambert J., Sinking coastal cities. Proc. IAHS 372, 189–198 (2015). [ Google Scholar ] 17. Shirzaei M., Freymueller J., Törnqvist T. E., Galloway D. L., Dura T., Minderhoud P. S., Measuring, modelling and projecting coastal land subsidence. Nat. Rev. Earth Environ. 2, 40–58 (2021). [ Google Scholar ] 18. Syvitski J. P. M., Kettner A. J., Overeem I., Hutton E. W. H., Hannon M. T., Brakenridge G. R., Day J., Vörösmarty C., Saito Y., Giosan L., Nicholls R. J., Sinking deltas due to human activities. Nat. Geosci. 2, 681–686 (2009). [ Google Scholar ] 19. Shirzaei M., Bürgmann R., Global climate change and local land subsidence exacerbate inundation risk to the San Francisco Bay Area. Sci. Adv. 4, eaap9234 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 20. Ohenhen L. O., Shirzaei M., Ojha C., Sherpa S. F., Nicholls R. J., Disappearing cities on US coasts. Nature 627, 108–115 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 21. Wöppelmann G., Marcos M., Vertical land motion as a key to understanding sea level change and variability. Rev. Geophys. 54, 64–92 (2016). [ Google Scholar ] 22. Marcos M., Wöppelmann G., Matthews A., Ponte R. M., Birol F., Ardhuin F., Coco G., Santamaría-Gómez A., Ballu V., Testut L., Chambers D., Coastal sea level and related fields from existing observing systems. Surv. Geophys. 40, 1293–1317 (2019). [ Google Scholar ] 23. Dada O. A., Almar R., Morand P., Bergsma E. W. J., Angnuureng D. B., Minderhoud P. S. J., Future socioeconomic development along the West African coast forms a larger hazard than sea level rise. Commun. Earth Environ. 4, 150 (2023). [ Google Scholar ] 24. B. Fox-Kemper, H. T. Hewitt, C. Xiao, Ocean, cryosphere and sea level change, in Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (Cambridge Univ. Press, 2021); 10.1017/9781009157896.011. [ DOI ] [ Google Scholar ] 25. Yulianto F., Wibowo M., Yananto A., Perdana D. H. F., Wiguna E. A., Prabowo Y., Rahili N., Nurwijayanti A., Iswari M. Y., Ratnasari E., Rusdiutomo A., Nugroho S., Purwoko A. S., Aziz H., Fachrudin I., Coastal vulnerability assessment using the machine learning tree-based algorithms modeling in the north coast of Java, Indonesia. Earth Sci. Inform. 16, 3981–4008 (2023). [ Google Scholar ] 26. Abidin H. Z., Andreas H., Gumilar I., Fukuda Y., Pohan Y. E., Deguchi T., Land subsidence of Jakarta (Indonesia) and its relation with urban development. Nat. Hazards 59, 1753–1771 (2011). [ Google Scholar ] 27. Chaussard E., Amelung F., Abidin H., Hong S. H., Sinking cities in Indonesia: ALOS PALSAR detects rapid subsidence due to groundwater and gas extraction. Remote Sens. Environ. 128, 150–161 (2013). [ Google Scholar ] 28. S. T. Lee, G. Sihombing, “The world’s fastest-sinking megacity has one last chance to save itself,” Bloomberg , 9 December 2023. 29. Abidin H. Z., Gumilar I., Andreas H., Murdohardono D., Fukuda Y., On causes and impacts of land subsidence in Bandung Basin, Indonesia. Environ. Earth Sci. 68, 1545–1553 (2013). [ Google Scholar ] 30. Susilo S., Salman R., Hermawan W., Widyaningrum R., Wibowo S. T., Lumban-Gaol Y. A., Meilano I., Yun S.-H., GNSS land subsidence observations along the northern coastline of Java, Indonesia. Sci. Data 10, 421 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 31. Sidiq T. P., Gumilar I., Abidin H. Z., Meilano I., Purwarianti A., Lestari R., Spatial distribution and monitoring of land subsidence using Sentinel-1 SAR data in Java, Indonesia. Appl. Sci. 15, 3732 (2025). [ Google Scholar ] 32. Abidin H. Z., Andreas H., Gumilar I., Wibowo I. R., On correlation between urban development, land subsidence and flooding phenomena in Jakarta. Proc. Int. Assoc. Hydrol. Sci. 370, 15–20 (2015). [ Google Scholar ] 33. Cigna F., Tapete D., Present-day land subsidence rates, surface faulting hazard and risk in Mexico City with 2014–2020 Sentinel-1 IW InSAR. Remote Sens. Environ. 253, 112161 (2021). [ Google Scholar ] 34. Sarah D., Hutasoit L. M., Delinom R. M., Sadisun I. A., Wirabuana T., A physical study of the effect of groundwater salinity on the compressibility of the Semarang-DemakAquitard, Java Island. Geosciences 8, 130 (2018). [ Google Scholar ] 35. Aditiya A., Ito T., Present-day land subsidence over Semarang revealed by time series InSAR new small baseline subset technique. Int. J. Appl. Earth Obs. Geoinf. 125, 103579 (2023). [ Google Scholar ] 36. Azeriansyah R., Ching K. E., Yuwono B. D., Impact of land subsidence-induced three-dimensional surface deformation on infrastructure in the Semarang-Demak Alluvial Plain, Indonesia. J. Earth Mar. Technol. 5, 111–122 (2025). [ Google Scholar ] 37. Rygus M., Bianchi M., Novellino A., Hussain E., Taufiq A., Rusli S. R., Meisina C., Permanent aquifer storage loss from long-term groundwater withdrawal: A case study of subsidence in Bandung (Indonesia). J. Hydrol. Reg. Stud. 57, 102129 (2025). [ Google Scholar ] 38. Shirzaei M., Rudolph M. L., Manga M., Deep and shallow sources for the Lusi mud eruption revealed by surface deformation. Geophys. Res. Lett. 42, 5274–5281 (2015). [ Google Scholar ] 39. Kopp R. E., Garner G. G., Hermans T. H. J., Jha S., Kumar P., Reedy A., Slangen A. B. A., Turilli M., Edwards T. L., Gregory J. M., Koubbe G., Levermann A., Merzky A., Nowicki S., Palmer M. D., Smith C., The framework for assessing changes to sea-level (FActS) v1.0: A platform for characterizing parametric and structural uncertainty in future global, relative, and extreme sea-level change. Geosci. Model Dev. 16, 7461–7489 (2023). [ Google Scholar ] 40. Ohenhen L. O., Shirzaei M., Land subsidence hazard and building collapse risk in the coastal city of Lagos, West Africa. Earths Future 10, e2022EF003219 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 41. L. Ley, Building on borrowed time: Rising seas and failing infrastructure in Semarang (Univ. of Minnesota Press, 2021). [ Google Scholar ] 42. Esteban M., Takagi H., Jamero L., Chadwick C., Avelino J. E., Mikami T., Fatma D., Yamamoto L., Thao N. D., Onuki M., Woodbury J., Valenzuela V. P. B., Crichton R. N., Shibayama T., Adaptation to sea level rise: Learning from present examples of land subsidence. Ocean Coast. Manag. 189, 104852 (2020). [ Google Scholar ] 43. Jiang H., Zhang J., Liu Y., Li J., Fang Z. N., Does flooding get worse with subsiding land? Investigating the impacts of land subsidence on flood inundation from Hurricane Harvey. Sci. Total Environ. 865, 161072 (2023). [ DOI ] [ PubMed ] [ Google Scholar ] 44. Abidin H. Z., Andreas H., Gumilar I., Sidiq T. P., Fukuda Y., Land subsidence in coastal city of Semarang (Indonesia): Characteristics, impacts and causes. Geomat. Nat. Hazards Risk 4, 226–240 (2013). [ Google Scholar ] 45. Minderhoud P. S. J., Middelkoop H., Erkens G., Stouthamer E., Groundwater extraction may drown mega-delta: Projections of extraction-induced subsidence and elevation of the Mekong delta for the 21st century. Environ. Res. Commun. 2, 111004 (2020). [ Google Scholar ] 46. Zhang M., Nicholls R. J., Wen J., AghaKouchak A., Bouma T. J., Darby S. E., Du S., Dai Z., Growing compound-flood risk, driven by both climate change and land subsidence, challenges flood risk reduction in major delta cities. One Earth 8, 101489 (2025). [ Google Scholar ] 47. Widodo J., Trihatmoko E., Setyaningrum N., Izumi Y., Handika R., Ardha M., Arief R., Sobue S., Nurlinda N., Pranantya P. A., Wiranu J. R., Khomarudin M. R., Technical and policy analysis: Time series of land subsidence for the evaluation of the Jakarta groundwater-free zone. Urban Sci. 9, 67 (2025). [ Google Scholar ] 48. Cao A., Esteban M., Valenzuela V. P. B., Onuki M., Takagi H., Thao N. D., Tsuchiya N., Future of Asian deltaic megacities under sea level rise and land subsidence: Current adaptation pathways for Tokyo, Jakarta, Manila, and Ho Chi Minh City. Curr. Opin. Environ. Sustain. 50, 87–97 (2021). [ Google Scholar ] 49. Phien-wej N., Giao P. H., Nutalaya P., Land subsidence in Bangkok, Thailand. Eng. Geol. 82, 187–201 (2006). [ Google Scholar ] 50. Naish T., Levy R., Hamling I., Hreinsdóttir S., Kumar P., Garner G. G., Kopp R. E., Golledge N., Bell R., Paulik R., Lawrence J., Denys P., Gillies T., Bengtson S., Howell A., Clark K., King D., Litchfield N., Newnham R., The significance of interseismic vertical land movement at convergent plate boundaries in probabilistic sea-level projections for AR6 scenarios: The New Zealand case. Earths Future 12, e2023EF004165 (2024). [ Google Scholar ] 51. Govorcin M., Bekaert D. P. S., Hamlington B. D., Sangha S. S., Sweet W., Variable vertical land motion and its impacts on sea level rise projections. Sci. Adv. 11, eads8163 (2025). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 52. Oelsmann J., Marcos M., Passaro M., Sanchez L., Dettmering D., Dangendorf S., Seitz F., Regional variations in relative sea-level changes influenced by nonlinear vertical land motion. Nat. Geosci. 17, 137–144 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 53. Fang J., Nicholls R. J., Brown S., Lincke D., Hinkel J., Vafeidis A. T., Du S., Zhao Q., Liu M., Shi P., Benefits of subsidence control for coastal flooding in China. Nat. Commun. 13, 6946 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 54. Ohenhen L. O., Zhai G., Lucy J., Werth S., Carlson G., Khorrami M., Onyike F., Sadhasivam N., Tiwari A., Ghobadi-Far K., Sherpa S. F., Lee J. C., Zehsaz S., Shirzaei M., Land subsidence risk to infrastructure in US metropolises. Nat. Cities 2, 543–554 (2025). [ Google Scholar ] 55. Shirzaei M., A wavelet-based multitemporal DInSAR algorithm for monitoring ground surface motion. IEEE Geosci. Remote Sens. Lett. 10, 456–460 (2013). [ Google Scholar ] 56. Shirzaei M., Bürgmann R., Fielding E. J., Application of wavelet-based InSAR time-series analysis for measuring surface deformation in Los Angeles. J. Geophys. Res. Solid Earth 120, 211–227 (2015). [ Google Scholar ] 57. Shirzaei M., Bürgmann R., Fielding E. J., c. Geophys. Res. Lett. 44, 2733–2742 (2017). [ Google Scholar ] 58. C. Werner, U. Wegmüller, T. Strozzi, A. Wiesmann, Gamma SAR and interferometric processing software, in ERS - Envisat Symposium (2000). [ Google Scholar ] 59. Lee J.-C., Shirzaei M., Novel algorithms for pair and pixel selection and atmospheric error correction in multitemporal InSAR. Remote Sens. Environ. 286, 113447 (2023). [ Google Scholar ] 60. Costantini M., A novel phase unwrapping method based on network programming. IEEE Trans. Geosci. Remote Sens. 36, 813–821 (1998). [ Google Scholar ] 61. M. Costantini, P. A. Rosen, Phase unwrapping techniques, in Proceedings of the IEEE International Geoscience and Remote Sensing Symposium (IGARSS) (IEEE, 1999), pp. 267–269. [ Google Scholar ] 62. Shirzaei M., Walter T. R., Estimating the effect of satellite orbital error using wavelet-based robust regression applied to InSAR deformation data. IEEE Trans. Geosci. Remote Sens. 49, 4600–4605 (2011). [ Google Scholar ] 63. Shirzaei M., Bürgmann R., Topography correlated atmospheric delay correction in radar interferometry using wavelet transforms. Geophys. Res. Lett. 39, L01405 (2012). [ Google Scholar ] 64. Blackwell E., Shirzaei M., Ojha C., Werth S., Tracking California’s sinking coast from space: Implications for relative sea-level rise. Sci. Adv. 6, eaba4551 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 65. R. F. Hanssen, Radar Interferometry: Data Interpretation and Error Analysis (Kluwer Academic Publishers, 2001). [ Google Scholar ] 66. Ohenhen L. O., Shirzaei M., Ojha C., Kirwan M. L., Hidden vulnerability of US Atlantic coast to sea-level rise due to vertical land motion. Nat. Commun. 14, 2038 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 67. Blewitt G., Hammond W. C., Kreemer C., Harnessing the GPS data explosion for interdisciplinary science. Eos 99, EO104623 (2018). [ Google Scholar ] 68. E. M. Mikhail, F. E. Ackermann, Observations and Least Squares (IEP Series in Civil Engineering, 1976). 69. Anantrasirichai N., Biggs J., Albino F., Bull D., The application of convolutional neural networks to detect slow, sustained deformation in InSAR time series. Geophys. Res. Lett. 46, 11850–11858 (2019). [ Google Scholar ] 70. Radman A., Akhoondzadeh M., Hosseiny B., Integrating InSAR and deep-learning for modeling and predicting subsidence over the adjacent area of Lake Urmia, Iran. GISci. Remote Sens. 58, 1413–1433 (2021). [ Google Scholar ] 71. Festa D., Novellino A., Hussain E., Bateson L., Casagli N., Confuorto P., del Soldato M., Raspini F., Unsupervised detection of InSAR time series patterns based on PCA and K-means clustering. Int. J. Appl. Earth Obs. Geoinf. 118, 103276 (2023). [ Google Scholar ] 72. Rygus M., Novellino A., Hussain E., Syafiudin F., Andreas H., Meisina C., A clustering approach for the analysis of InSAR time series: Application to the Bandung Basin (Indonesia). Remote Sens. 15, 3776 (2023). [ Google Scholar ] 73. Tiwari A., Shirzaei M., A novel machine learning and deep learning semi-supervised approach for automatic detection of InSAR-based deformation hotspots. Int. J. Appl. Earth Obs. Geoinf. 126, 103611 (2024). [ Google Scholar ] 74. Hotelling H., Analysis of a complex of statistical variables into principal components. J. Educ. Psychol. 24, 417–441 (1933). [ Google Scholar ] 75. M. Ester, H. P. Kriegel, J. Sander, X. Xu, A density-based algorithm for discovering clusters in large spatial databases with noise (in kdd, 1996), vol. 96, pp. 226–231. 76. D. J. Berndt, J. Clifford, Using dynamic time warping to find patterns in time series, in Proceedings of the 3rd International Conference on Knowledge Discovery and Data Mining (1994), vol 359–370. [ Google Scholar ] 77. Indonesia Geospasial, “Download Geological Shapefiles of All of Indonesia”; https://www.indonesia-geospasial.com [accessed 12 February 2025]. 78. D. Zanaga, R. Van De Kerchove, D. Daems, W. De Keersmaecker, C. Brockmann, G. Kirches, J. Wevers, O. Cartus, M. Santoro, S. Fritz, M. Lesiv, M. Herold, N.-E. Tsendbazar, P. Xu, F. Ramoino, O. Arino, ESA WorldCover 10 m 2021 v200 (Zenodo, 2022); 10.5281/zenodo.7254221. [ DOI ] 79. Copernicus Marine Environment Monitoring Service, “Global Ocean Gridded L4 Sea Surface Heights and Derived Variables Reprocessed 1993 Ongoing”; https://marine.copernicus.eu/ [accessed 12 March, 2025]. 80. Grégoire M., Alvera-Azcaráte A., Buga L., Capet A., Constantin S., D’ortenzio F., Doxaran D., Faugeras Y., Garcia-Espriu A., Golumbeanu M., González-Haro C., González-Gambau V., Kasprzyk J. P., Ivanov E., Mason E., Mateescu R., Meulders C., Olmedo E., Pons L., Pujol M. I., Sarbu G., Turiel A., Vandenbulcke L., Rio M. H., Monitoring Black Sea environmental changes from space: New products for altimetry, ocean colour and salinity. Potentialities and requirements for a dedicated in-situ observing system. Front. Mar. Sci. 9, 998970 (2023). [ Google Scholar ] 81. L. Ma, S. Li, D. Sun, NOAA JPSS Program Office, NOAA JPSS Visible Infrared Imaging Radiometer Suite (VIIRS) Level-3 Global Flood Map from NDE [Data set] (NOAA National Centers for Environmental Information, 2024). 82. Aribowo S., Husson L., Natawidjaja D. H., Authemayou C., Daryono M. R., Puji A. R., Valla P. G., Pamumpuni A., Wardhana D. D., de Gelder G., Djarwadi D., Lorcery M., Active back-arc thrust in north west Java, Indonesia. Tectonics 41, e2021TC007120 (2022). [ Google Scholar ] 83. Whincup P., Oktavianus R., Egy C., Jakarta groundwater: Victim of its own success. Hydrgeol. J. 32, 353–358 (2024). [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials Figs. S1 to S21 Legends for table S1 to S3 Legend for data S1 References sciadv.aec0172_sm.pdf (3.2MB, pdf) Tables S1 to S3 sciadv.aec0172_tables_s1_to_s3.zip (38.3KB, zip) Data Availability Statement The VLM and horizontal land motion data for Java Island generated for this study are available at https://doi.org/10.5281/zenodo.15786356 . The medium and low confidence projected RSL change and rates for five SSPs (SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) generated for this study are available at https://doi.org/10.5281/zenodo.15786356 . Sentinel-1 data used in this study are publicly available through the Alaska Satellite Facility and can be accessed at https://search.asf.alaska.edu . The satellite altimetry data for sea level rise are available from Copernicus Marine Environment Monitoring Service (CMEMS) and are available at https://data.marine.copernicus.eu/product/SEALEVEL_GLO_PHY_L4_MY_008_047/download . Access requires the creation of a free CMEMS user account, which is available to all users at no cost upon registration. The population for Java Island was estimated using the WorldPop dataset available at https://hub.worldpop.org/geodata/summary?id=77690 . The coastal elevation referenced to mean sea level uses the DeltaDTM dataset v1.1 available at https://doi.org/10.4121/21997565 . All other data and code needed to evaluate and reproduce the results in the paper are present in the paper and/or the Supplementary Materials. This study did not generate new materials. Articles from Science Advances are provided here courtesy of American Association for the Advancement of Science ACTIONS View on publisher site PDF (1.9 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top