Knowledge From the Afrikan World, For the World

PMC Articles

Global subsidence of river deltas

PMCID: PMC12823437

PMID: 41535473

Abstract

River deltas sustain dense human populations, major economic centres and vital ecosystems worldwide 1 , 2 . Rising sea levels and subsiding land threaten the sustainability of these valuable landscapes with relative sea-level rise and associated flood, land loss and salinization hazards 1 – 3 . Despite these risks, vulnerability assessments are impeded by the lack of contemporary, high-resolution, delta-wide subsidence observations 4 . Here we present spatially variable surface-elevation changes across 40 global deltas using interferometric synthetic aperture radar. Using this dataset, we quantify delta surface-elevation loss and show the prevalence and severity of subsidence in river deltas worldwide. Our analysis of three key anthropogenic drivers of delta elevation changes shows that groundwater storage has the strongest relative influence on vertical land motion in 10 of the 40 deltas. The other deltas are either influenced by multiple drivers or dominated by sediment flux or urban expansion. Furthermore, we find that contemporary subsidence surpasses absolute (geocentric) sea-level rise as the dominant driver of relative sea-level rise for most deltas over the twenty-first century. These findings suggest the need for targeted interventions addressing subsidence as an immediate and localized challenge, in parallel with broader efforts to mitigate and adapt to climate change-driven global sea-level rise. Spatially variable surface-elevation changes across 40 global deltas using interferometric synthetic aperture radar are reported.

Full Text

This recognized importance, which makes deltas indispensable, also increases their exposure to compounding climatic, environmental and anthropogenic threats. As low-lying landforms, with extensive areas less than 2 m above sea level, deltas are acutely susceptible to rising sea level, storm surge, land subsidence, shifting temperature and rainfall patterns, and other environmental pressures, which are amplified by climate change. These pressures degrade agricultural land; disrupt freshwater availability; exacerbate coastal and fluvial flooding; promote wetland loss, saltwater intrusion and shoreline retreat; and threaten infrastructure in deltas. Beyond direct physical impacts, the interplay of these hazards also creates potential cascading socioeconomic consequences. For example, land loss and freshwater scarcity may drive displacement and migration, heightening competition for dwindling resources and fuelling social tensions. Together, these intersecting climatic, environmental, human-driven pressures and multi-hazards render deltas the most fragile landscapes on Earth, with their low elevation and high urban exposure placing them at the forefront of climate and environmental risks (Extended Data Fig. 1).
Here, we present high-spatial-resolution datasets of surface-elevation change derived from Sentinel-1 synthetic aperture radar (SAR) interferometry across 40 deltas globally (Fig. 1). These datasets capture delta-wide temporal trends, subsidence rates and horizontal motion at 75 m resolution, spanning five continents and 29 countries. Our analysis encompasses all major river deltas with a population exceeding 3 million people, historically recognized sinking deltas and representatives of less-populated, understudied deltas of regional ecological and economic importance (Methods).
Source data
We measured the spatial patterns and rates of subsidence in 40 deltas by analysing the complete archive of the Sentinel-1 SAR dataset between 2014 and 2023 using advanced multitemporal interferometric SAR (InSAR) analysis (Methods). InSAR measures surface-elevation changes, capturing vertical land motion (VLM), sediment deposition and erosional processes. For consistency, to reflect both VLM and surface-elevation change in the deltas, we use the terms VLM or elevation gain or loss to describe net surface-elevation change across all delta environments, with positive values indicating uplift or elevation gain and negative values indicating subsidence or net elevation loss. Throughout this study, negative VLM is quoted with negative signs and references land subsidence rates, whereas only the absolute values are reported when presenting subsidence rates.
Our analysis shows that subsidence threatens deltas globally, with the delta-scale average rate of VLM on all deltas indicating subsidence (Fig. 1). In 12 out of 40 deltas, the average sinking rate is moderate, at less than 2 mm yr−1. By contrast, more than half of the deltas exhibit subsidence rates exceeding 3 mm yr−1, and in 13 of these deltas (Nile, Po, Vistula, Ceyhan, Brahmani, Mahanadi, Chao Phraya, Mekong, Red, Ciliwung, Brantas, Godavari and Yellow River), the average subsidence rates exceed the current estimates of global SLR (that is, about 4 mm yr−1). Among these, the Chao Phraya (Thailand), Brantas (Indonesia) and Yellow River (China) deltas show an average sinking rate of more than twice the current global SLR rate. To further highlight the severity of subsidence in deltas, we compared the subsidence with the regional geocentric SLR rates for the twenty-first century (2001–present). In 18 of the 40 deltas (the Nile, Po, Vistula, Ceyhan, Rioni, Brahmani, Mahanadi, Ganges–Brahmaputra, Godavari, Chao Phraya, Mekong, Red River, Ciliwung, Brantas, Amazon, Parana, Pearl and Yellow River), the average rate of local land subsidence is greater than the rate of regional geocentric SLR (Fig. 1 and Supplementary Table 1). However, in almost every delta (except Rio Grande) at least 1% of the delta area is subsiding faster than both global and geocentric sea levels (Fig. 1 and Supplementary Table 1).
Among all deltas, we find that at least 35% of the area is sinking, and in 38 deltas (excluding Neva and Fraser), more than 50% of the delta area is sinking (Fig. 2a). Of the 40 deltas, 19 show widespread subsidence patterns, with greater than 90% of the delta area affected by subsidence (for example, Mississippi, Niger, Nile, Rhine–Meuse, Po, Vistula, Brahmani, Mahanadi, Ganges–Brahmaputra, Chao Phraya, Mekong and Brantas deltas). Deltas with notable subsiding areas with greater than 50% of the delta area sinking faster than 5 mm yr−1 include the Chao Phraya (94% of delta area), Nile (80%), Brahmani (77%), Po (74%), Mahanadi (69%), Brantas (66%), Vistula (57%), Yellow River (53%) and Mekong (51%) deltas (Fig. 2a and Supplementary Table 1). In sum, we estimate that a total delta area of 460,370 km2 is exposed to subsidence. If we consider a global habitable geomorphic area of 710,000–855,000 km2 for deltas, approximately 54–65% of global delta areas are sinking just from the analysis of the 40 deltas. By region, South Asia, East Asia and Southeast Asia, with 17 representative deltas, have the greatest exposure to subsidence, with 274,000 km2 of delta area subsiding. Africa, South America, North America and Europe have total subsiding delta areas of 78,800 km2, 39,800 km2, 37,800 km2, and 30,000 km2, respectively. Seven large deltas—Ganges–Brahmaputra, Nile, Mekong, Yangtze, Amazon, Irrawaddy and Mississippi deltas—contribute about 57% of the total subsiding delta area, with a combined area of 265,000 km2. Coastal cities such as Alexandria (Nile), Bangkok (Chao Phraya), Dhaka and Kolkata (Ganges–Brahmaputra), Shanghai (Yangtze), Yangon (Irrawaddy), Cần Thá (Mekong), Thái Bình (Red River), Niigata (Chikuma-gawa), Jakarta (Ciliwung), Surabaya (Brantas) and Dongying (Yellow River) are experiencing subsidence at rates equal to or exceeding the delta-wide averages, indicative of the intensity of subsidence and elevation loss processes in cities on deltas.
a, Proportion of each delta exposed to different rates of subsidence. Note that only subsiding areas are represented in each bar, and areas of uplift within each delta are omitted to emphasize the extent of elevation loss. b–m, Spatial maps of VLM rates for the Fraser (Canada) (b), Mississippi (the USA) (c), Parana (Argentina) (d), Niger (Nigeria) (e), Nile (Egypt) (f), Po (Italy) (g), Ganges–Brahmaputra (India–Bangladesh) (h), Chao Phraya (Thailand) (i), Mekong (Vietnam) (j), Red River (Vietnam) (k), Pearl (China) (l), Yellow River (China) (m) deltas. Positive VLM (green–purple hues) suggests uplift or elevation gain, whereas negative VLM (yellow–orange–red hues) indicates land subsidence. The spatial VLM maps for the other 28 deltas are shown in Extended Data Figs. 2–4. Background image in b–m is Esri, streets-dark. Scale bars, 5 km (b); 50 km (c,f,h,j); 20 km (d,e,i,k,l,m); 10 km (g).
Source data
Furthermore, we observe non-uniform spatially variable VLM within individual deltas, reflecting the complex interplay of natural and anthropogenic processes (Fig. 2 and Extended Data Figs. 2–4). Although all deltas exhibit an overall trend of subsidence, localized and broad zones of uplift, which vary from 0 mm yr−1 to greater than 5 mm yr−1 are observed in some areas (Fig. 2b,d,k,m, and Extended Data Figs. 2e,f,i,j,l and 3c,f). In some deltas (for example, Wouri, Zambezi, Indus, Ciliwung and Yellow River), the observed uplift or elevation-gaining parts correlate with patterns of horizontal land motion (Extended Data Figs. 5–7). Possible mechanisms may include sediment redistribution processes potentially driven by river dynamics or growth faulting, either of which can cause localized zones of elevation gain even within a predominantly subsiding deltaic system. This highlights the necessity of comprehensive assessments and models of delta vulnerability to consider not only overall absolute subsidence rates but also the spatial heterogeneity of elevation change dynamics.
To quantify the relative contributions of anthropogenic factors to delta subsidence and elevation loss, we analysed the relationship between three main anthropogenic drivers—groundwater storage change, sediment flux alteration and urban expansion—and non-glacial isostatic adjustment VLM/subsidence rates across the 40 deltas (Methods and Supplementary Table 2).
Figure 3a and Extended Data Fig. 8 show the interplay of anthropogenic factors and their correlation with subsidence rates across the 40 deltas. Deltas experiencing groundwater storage (GWS) loss (indicative of groundwater extraction), negative sediment flux change (red and yellow hues; reflecting sediment reduction due to upstream human activities) and higher urban population growth tend to have higher rates of subsidence (for example, the Yellow River, Po, Nile, Chao Phraya and Mekong deltas). Conversely, deltas with GWS stability or gain (net increase in groundwater storage), positive sediment flux change (blue colours; sediment surplus) and limited urban expansion show lower subsidence rates (for example, Saloum, Amazon and Ogooué deltas).
a, Bubble plot showing the relationship between VLM rates and anthropogenic drivers across deltas. Plot shows VLM rate (mm yr−1) against GWS rate (mm yr−1). Bubble colours represent sediment flux change (%), in which positive values (blue colours) indicate increased sediment supply due to human activities (and thereby increased potential to gain elevation and compensate subsidence-induced elevation loss), whereas negative values (yellow–orange–red colours) indicate a decline in sediment availability. Bubble size indicates urban fraction change (%), with larger circles representing a greater urban expansion over the twenty-first century. The dashed line represents the MLR fit. See Extended Data Fig. 8 for individual pairwise relationships between each anthropogenic driver and VLM. b, Ternary plot of subsidence rates with nLIME scores.
Source data
The initial multilinear regression (MLR) model, which included interaction terms between the different anthropogenic factors, poorly captured subsidence dynamics on the deltas (R2 = 0.2 ± 0.1), as it failed to account for nonlinear interactions between the different processes (Fig. 3a). For instance, urban expansion not only directly increases infrastructure loading but also indirectly elevates groundwater demand, thereby compounding aquifer depletion and extraction-induced subsidence, which are synergistic effects that linear models cannot resolve.
To address these limitations, we used a random forest (RF) machine learning approach designed to capture nonlinear relationships and variable interactions. The RF model shows a moderate to strong relationship between the predictors (GWS, sediment flux and urban expansion) and VLM, achieving improved performance over the MLR model (R2 = 0.6 ± 0.1; RMSE (root mean square error) = 1.9 ± 0.1 mm yr−1; MAE (mean absolute error) = 1.4 ± 0.2 mm yr−1), and capturing complex, non-additive relationships between anthropogenic stressors and subsidence rates (Fig. 3a and Supplementary Fig. 1). However, we observe some underestimation at high subsidence rates (>8.0 mm yr−1) (Supplementary Fig. 1), which probably suggests that natural processes or other anthropogenic predictors (not considered in our analysis) may contribute to subsidence in these highly dynamic deltaic environments.
Note that the primary objective in our analysis is not to predict subsidence rates across deltas, but rather to identify and extract key features that explain the dynamic relationships between the three anthropogenic drivers and subsidence across these deltas. Feature importance analysis from the RF model identifies GWS as the dominant anthropogenic predictor of delta subsidence (0.5 ± 0.2), whereas sediment flux change (0.3 ± 0.2) and urbanization (0.3 ± 0.1) have secondary roles as subsidence rate predictors across these deltas (Fig. 3a and Supplementary Fig. 1b). However, the large standard deviations in feature importance values reflect substantial variability in predictor dominance across subsampled delta subsets, suggesting that the primary contributors to subsidence differ locally depending on the anthropogenic or geomorphic context. To resolve delta-specific mechanisms, we applied local interpretable model-agnostic explanations (LIME), which interprets individual predictions by approximating the RF model locally with simpler, interpretable functions. Deltas with low LIME model fidelity (R2 < 0.5) were excluded from this interpretative analysis, refining the dataset from 40 to 28 deltas (Methods). The low fidelity scores for some deltas could be due to unaccounted processes (natural and/or other anthropogenic) in our RF model. The retained 28 deltas show improved overall model performance (R2 = 0.7 ± 0.1; RMSE = 0.4 ± 0.1 mm yr−1; MAE = 0.3 mm yr−1), ensuring reliable interpretation of local feature importance. Normalized LIME feature importance scores (nLIME) showed substantial heterogeneity in predictor dominance (Supplementary Table 2). GWS emerged as the most significant factor across the different deltas (0.6 ± 0.3), whereas sediment flux change (0.3 ± 0.1) and urbanization (0.1 ± 0.1) exhibited lower but context-dependent impacts (Supplementary Fig. 1b).
To assess the dominant influence on land motion across individual deltas, the nLIME for each delta was mapped onto a ternary diagram (Fig. 3b). Of the 28 deltas, 35%, including the Mekong, Ganges–Brahmaputra, Rhine–Meuse, Fraser, Cauvery, Irrawaddy and Red River systems, cluster within the GWS portion of the diagram (nLIMEGWS > 0.7), suggesting that observed GWS changes in these deltas are the primary driver of subsidence among the three anthropogenic variables examined (Fig. 3b and Supplementary Table 2). The Chao Phraya and Yellow River deltas, with the highest average subsidence rates, plot near the centre of the ternary diagram, reflecting relatively balanced contributions from GWS, sediment flux and urban expansion. Sediment flux correlates most closely with elevation changes in deltaic systems, such as the Saloum, Mississippi, Amazon and Rio Grande deltas, suggesting that reduced sediment delivery may exacerbate land elevation loss in these deltas. The Nile, Po, Chikuma-gawa, Mahanadi, Kabani, Niger and Volta deltas exhibit mixed contributions from GWS, sediment flux changes and population change, with GWS slightly outweighing sediment deficits as predictors in the Nile and Po deltas, possibly reflecting reliance on aquifer-dependent irrigation. These findings are consistent with delta-specific studies that attribute accelerated subsidence in densely populated Asian deltas—Mekong, Ganges–Brahmaputra and Chao Phraya—to urbanization and unsustainable groundwater extraction for agriculture, industry and domestic use. Moreover, the Nile, Po and Mississippi deltas, which were historically sustained by seasonal floods that deposited sediments, are now documented to experience severe sediment deficits due to dams and levees, accelerating elevation loss.
To quantify the contributions of SLR and land subsidence in deltas, we evaluated their relative impact on the exposed delta populations. Our analysis shows that current average subsidence rates exceed geocentric SLR in 18 of the 40 deltas, including the Nile, Mekong, Red River, Ganges–Brahmaputra, Brahmani, Mahanadi, Chao Phraya, Ciliwung, Brantas and Yellow River deltas, affecting approximately 236 million people—a population about 50% larger than those residing in deltas in which the current rates of geocentric SLR outpace the subsidence rates (156.9 million) (Fig. 4a). This disparity is particularly pronounced for vulnerable populations occupying land below 1 m elevation. In these lowest elevation areas, subsidence dominates the contribution to RSLR in about two-thirds of the deltas, including Amazon, Fraser, Niger, Rhone, Vistula, Ganges–Brahmaputra, Mekong, Red River, Pearl, Yangtze and Godavari deltas (Fig. 4b). Of the 76 million people living in delta areas with an elevation below 1 m, 84% (63.7 million people) reside in rapidly sinking areas of the deltas (Fig. 4b). These observations are striking, revealing the current dominance of subsidence over geocentric SLR in global deltas. Moreover, the spatial heterogeneity of VLM creates localized extreme rates of subsidence within deltas, further exacerbating their vulnerability. Under the current trajectory, moderate emission scenarios (shared socioeconomic pathway 2-4.5 (SSP2-4.5)), current maximum subsidence rates in the deltas already surpass projected twenty-first-century SLR rates (no VLM). Through the end of the twenty-first century, current maximum subsidence rates in all 40 deltas exceed projected SLR rates (Fig. 4c). This disparity extends to the 95th percentile subsidence rates, representing widespread, high-magnitude sinking across the deltas. In 29 deltas, 95th percentile subsidence rates exceed the projected SLR rates by 2050, outpacing SLR by 1.1 (Niger delta) to 10.3 (Yellow River delta) times. By 2100, as the current maximum rate of SLR (SSP2-4.5) accelerates to 0.9 cm yr−1, current 95th percentile subsidence rates still dominate in 22 deltas, surpassing geocentric SLR by up to seven times. Even accounting for worst-case, high-emission scenarios (SSP5-8.5), subsidence will exceed projected SLR rates in all deltas (considering maximum subsidence) and in 23 deltas (considering 95th percentile subsidence) through 2050. By 2100, current maximum subsidence rates exceed projected SLR in 38 of 40 deltas, whereas 95th percentile subsidence rates remain dominant in seven deltas (Godavari, Chao Phraya, Mekong, Ciliwung, Brantas, Red River and Yellow River) (Supplementary Table 1).
Source data
To visualize disparities in adaptive capacity and risk, we mapped global deltas into a two-dimensional (2D) impact matrix defined by RSLR and ND-GAIN adaptation readiness scores (Fig. 5). This framework allows for a comparative assessment of deltas assuming that the adaptation readiness of the delta is reflected by the adaptation readiness of its country, categorizing them into four quadrants: (1) Unprepared Divers (high RSLR (>4 mm yr−1), low readiness (<0.52)); (2) Rising Ready (high RSLR, high readiness (>0.52)); (3) Latent Threats (low RSLR, low readiness); and (4) Safe Havens (low RSLR, high readiness). 65% of the deltas (26 out of 40 deltas), predominantly in low- and middle-income nations, fall into the Unprepared Divers group, in which nations have a diminished adaptive capacity and RSLR rates exceeding current global SLR (Fig. 5a). These challenges are compounded for indigenous communities, who primarily live in the lowest-lying delta areas; lack the resources needed to implement large-scale adaptation; and face relocation barriers due to cultural and subsistence ties despite escalating risks.
Source data
Most deltas in high-income countries, including the Yellow River (China), Vistula (Poland), Po (Italy), Rhine–Meuse (the Netherlands) and Mississippi (the USA) deltas, cluster in the Rising Ready group, demonstrating robust governance (Fig. 5a). For example, the integrated flood management approach of the Dutch delta, which combines ecological restoration with infrastructural fortifications, has become a model for coastal hazard resilience. However, some deltas even within this group face substantial gaps. For instance, the Mississippi delta has lost more than 5,000 km2 of land (mainly wetlands) since 1932 because of a lack of adaptation (for example, sediment diversion projects), whereas the Po delta struggles with salinization driven by agricultural groundwater extraction, highlighting how economic priorities can undermine adaptation even in high-income regions. Although RSLR exceeds global rates of SLR in most deltas, exceptions exist. The Latent Threats group includes the Saloum and Neva deltas, which exhibit relatively low RSLR and low adaptive capacity (Fig. 5a), indicating their unpreparedness and potential vulnerability to a future rise in sea level (Fig. 4c). The Rioni and the Fraser delta fall into the Safe Havens group, in which lower RSLR is coupled with higher adaptive capacity, indicative of low risk and preparedness for current and future sea-level changes. The Rioni Delta is the only delta in our sample exhibiting negative sea-level trends for the twenty-first century, in which long-term regional sea-level decline masks short-term fluctuations (Methods).
To examine the evolving risk landscape, we compared twentieth-century and present-day impact matrices (Fig. 5b). For our analysis, we used tide gauge data to estimate twentieth-century RSLR rates, which were available for only 15 of the 40 deltas. Our estimates show that 10 deltas previously classified as Latent Threats (low RSLR, low readiness) and Safe Havens (low RSLR, high readiness) groups during the twentieth century have transitioned to Unprepared Divers (high RSLR, low readiness) and Rising Ready (high RSLR, high readiness) groups in the twenty-first century (Fig. 5b). This shift highlights the accelerating contemporary RSLR trends, driven by land subsidence and SLR. Deltas such as the Mississippi, Ganges–Brahmaputra and Mekong show sustained increases in long-term RSLR rates above 4.0 mm yr−1 since the twentieth century, exacerbating vulnerabilities in these densely populated regions. Conversely, the Chao Phraya and the Rioni deltas showed a decline in RSLR and improved adaptive capacity in the twentieth century. However, although the Rioni Delta exhibited a more than 200% decline in RSLR, the Chao Phraya Delta still experiences high RSLR rates (12.3 mm yr–1). The pronounced decrease in RSLR for the twentieth century in the Rioni Delta probably reflects localized subsidence at the tide gauge station rather than a delta-wide RSLR trend (Methods). The greatest change in RSLR was observed in the Nile Delta, surging from 1 mm yr–1 in the twentieth century to more than 10 mm yr–1 in the twenty-first century (Fig. 5b). Moreover, we find that all deltas in low- and middle-income countries in the present-day Unprepared Divers groups, transitioned from the Latent Threats group, suggesting stagnant adaptive capacity despite worsening RSLR. By contrast, deltas such as the Yangtze (China), Pearl (China) and Vistula (Poland) shifted from Latent Threats to Rising Ready, demonstrating increased adaptation readiness due to economic growth, raising governance and institutional capacity to adapt, although RSLR has surged (Fig. 5b). Although deltas in the Rising Ready quadrant showed potential for robust adaptation policies, deltas in the Unprepared Divers remain trapped in cycles of reactive, underfunded responses.
Our estimates show that globally, 42,000 km2 of the delta area at present lies below sea level, containing a population of 10.2 million people (Extended Data Fig. 1). The 35 deltas with the greatest exposure included in this analysis are Nile, Mississippi, Rhine–Meuse, Mekong, Niger, Cauvery, Po, Red River, Vistula, Rhone, Amazon, Ganges–Brahmaputra, Chao Phraya, Kabani, Pearl, Rio Grande, Yangtze, Yellow River, Senegal, Indus, Saloum, Grijalva, Ceyhan/Seyhan, Rioni, Cross, Chikuma-gawa, Volta, Brantas, Neva, Wouri, Irrawaddy, Ogooué, Zambezi, Magdalena and Ciliwung (Extended Data Fig. 1). The cumulative delta area and population below sea level are 38,000 km2 and 10.1 million people, respectively, reaching within rounding errors of the global total exposure. Deltas such as the Danube, Orinoco and Shatt-el-Arab met the selection criteria but were excluded due to challenges associated with the SAR imaging and interferometric analysis (including spatial coverage gaps, excessive temporal baselines, poor coherence and limited data availability). The five supplementary deltas are Brahmani, Mahanadi, Godavari, Parana and Fraser deltas.
We analysed 132 SAR frames from the Sentinel-1A/B C-band satellite, spanning September 2016 to May 2023. The SAR datasets include 3,300 images obtained in single-orbit geometry (ascending or descending) for 13 deltas and 10,700 images obtained in both ascending and descending orbits for 27 deltas. See Supplementary Table 3 for the complete inventory of SAR images used in each delta. For each SAR dataset, we applied a multi-looking factor of 32:6 (range:azimuth) to improve the signal-to-noise ratio, obtaining an average pixel resolution of about 75 m. To minimize decorrelation errors, we also constrained the interferometric pairs to a maximum temporal and perpendicular baselines of 300 days and 80 m, respectively. For deltas requiring multi-frame coverage (for example, Amazon, Mississippi, Mekong, Ganges–Brahmaputra, Nile, Red River and Niger), we arranged in a mosaic form the overlapping adjacent frames along a single path before processing or post-processed deltas with coverage spanning multiple paths to ensure full spatial continuity across expansive deltas.
In the 27 deltas with overlapping spatiotemporal SAR satellite coverage and different orbit geometries (ascending and descending), we estimate the horizontal (east–west) and VLM components of deformation by jointly inverting the LOS time series of the ascending and descending tracks. To this end, we identified the co-located pixels of the LOS time series by resampling the pixels from the descending track onto the ascending track to obtain two co-located LOS displacement velocities . Given and their associated variances are the LOS displacement and variances for a given pixel, the model to combine the LOS velocities to generate a high-resolution map of the east–west (E) and VLM (U) displacements are given bywhere, C represents the unit vectors for projecting (E) and (U) displacements onto the LOS, which is a function of the heading angle of the satellite and incidence angles of each pixel. The solution to the model in equation (1) is given bywhere X represents the unknowns (E) and (U), G is the design matrix comprising the unit vectors for projecting the horizontal and vertical displacements onto the line of sight, L are the observations , and P is the weight matrix, which is inversely proportional to the observant variances . To obtain the parameter variance–covariance matrix (QXX), we use the concept of error propagation to calculate the associated parameter uncertainties given the observation errors as follows:
The distribution of the standard deviations (precision of the results) for all pixels (20.5 million) across the 40 deltas is shown in Supplementary Fig. 2. The standard deviation distribution shows that 99% of the pixels have a value <0.5 mm yr–1. We evaluated the accuracy of the results by comparing the averaged VLM rates of pixels within a radius of 100 m with more than 100 independent GNSS data (that is, stations that were not used in the reference frame transformation). The validation included 122 GNSS stations across 23 deltas with historical long-term records (spanning various periods before and/or including the InSAR observation window) and 81 GNSS stations across 15 deltas with time series covering at least 70% of the InSAR observation period (2014–2023) (Supplementary Fig. 3). We found a strong correlation (0.7–0.8), between GNSS and InSAR velocities, with an RMSE of 1.4 mm yr–1 for long-term rates (Supplementary Fig. 3a) and 1.2 mm yr–1 for rates within the InSAR observation period (Supplementary Fig. 3b). The improved agreement for temporally coincident measurements suggests that nonlinear subsidence behaviour contributes to some scatter when comparing historical GNSS rates to contemporary InSAR measurements, although the overall correlation remains strong in both cases. Note that some GNSS stations used for validation, while within the broader processed SAR frame, are outside the clipped delta boundaries. Note that the final delta extents were delineated using a tiered approach. Primary boundaries were derived from ref. , supplemented by ref.  for deltas not covered in the former. For extensive deltas in which the entire delta surface is not analysed (for example, the Ganges–Brahmaputra), boundaries were defined using the SAR spatial extent.
For each delta, we extracted the ensemble mean and standard deviation in GIA vertical velocity to correct observed deformation rates and isolate contemporary, non-GIA contributions to VLM. Supplementary Table 1 shows the mean GIA-induced VLM, the associated standard deviation and the per cent contribution of GIA to the total observed VLM magnitude for each delta. GIA accounts for the largest proportion and exceeds (>100%) the total VLM in the Neva (540%) and Fraser (455%) deltas, in which low observed VLM rates are substantially influenced by strong GIA uplift. Moderate GIA contributions (25–55%) are observed in five deltas, including the Rio Grande, Mississippi, Volta, Rhine and Ogooué deltas. Most of the deltas (55%) exhibit minimal GIA influence, with contributions under 10%, indicating that observed VLM is primarily governed by contemporary anthropogenic and natural processes such as groundwater withdrawal, sediment compaction, or tectonics. In 28–67% (accounting for uncertainty) of the deltas, the sign and approximate magnitude of observed and GIA-corrected VLM are consistent, implying limited distortion from GIA and the sustained expression of contemporary processes on the average local subsidence. By contrast, the Fraser and Neva deltas illustrate how substantial GIA-induced uplift in high-latitude, post-glacial regions can obscure contemporary subsidence processes through opposing vertical trends. In both cases, modest observed subsidence rates (Fraser −0.4 mm yr−1 and Neva −0.2 mm yr–1) are counteracted by substantial GIA uplift of 1.8 ± 2.3 mm yr−1 and 1.0 ± 0.3 mm yr−1, respectively.
We analysed the relationship between major anthropogenic pressures on global deltas to subsidence and elevation loss by quantifying the contributions of groundwater storage change, sediment flux alteration and urban expansion to the residual rates of sinking (after GIA correction) across the 40 deltas. These globally consistent datasets provide insights into human-induced impacts on land subsidence and elevation change in river deltas (Supplementary Table 2).
To isolate GWS change from TWS, we used the 1/4° global land data assimilation system Noah model to remove changes in SMS and SWE contributions and used the WaterGAP Global Hydrology Model (WGHM v.2.2d) (refs. ) to remove SWS contributions. The contribution from SWE was negligible in most deltas, given their prevailing arid and semi-arid climate (Fig. 1), although it was included to maintain consistency across all deltas. SWS components include contributions from rivers, lakes, wetlands and reservoir storage within the GRACE footprint for each delta. The residual signal following removal of SWS, SMS and SWE was interpreted as the GWS anomaly.
Supplementary Figs. 4 and 5 compare the time-invariant model (black curves) with the stochastic-seasonal model (red curves) for GRACE-derived GWS and RSLR from tide gauges in the Mississippi and Chao Phraya deltas. These plots show that a stochastic-seasonal process better represents the observed variability in the time series. The post-fit residuals of the time-invariant model show some systematic seasonal patterns, particularly during periods when seasonal amplitudes deviate from the assumed constant values (Supplementary Figs. 4b,d and 5b,d). By contrast, the stochastic model accommodates time-dependent variations in seasonal amplitudes, resulting in reduced (often near-zero) residuals (Supplementary Figs. 4b,d and 5b,d), demonstrating the advantage of the stochastic-seasonal model in capturing transient seasonal variations rather than fixed annual and semiannual cycles.
The GWS rates for each delta are summarized in Supplementary Table 2, and Fig. 3a and Extended Data Fig. 8a show the relationship with the subsidence rates. Negative GWS trends indicate mass depletion, primarily driven by groundwater extraction, whereas positive trends represent net groundwater accumulation due to recharge processes, reduced extraction or hydrological interventions. To evaluate the reliability of GRACE-derived GWS trends, we compared them with in situ groundwater level trends for 18 deltas (Supplementary Fig. 6). Groundwater levels were compiled from two publicly available sources: 13 deltas from ref.  and 5 deltas from the Global Groundwater Monitoring Network. Given the spatial scale discrepancy between GRACE (basin-wide) and well observations (point-scale), we emphasized agreement in trend direction rather than absolute magnitudes. Each site was categorized based on the sign of the GRACE and well trends, and a confusion matrix was constructed to assess consistency. The analysis yielded an overall classification accuracy of 88.9%, with six sites exhibiting positive–positive trends (PPT) and 10 showing negative–negative trends (NNT). Only two sites showed mixed behaviour (NPT or PNT), and no site exhibited fully opposing trends. Moreover, a high correlation (R = 0.7) was observed between the GRACE-based GWS and well-derived trends, further supporting the consistency of GRACE estimates at the basin scale despite localized variability in in situ measurements. Although the coarse spatial resolution of GRACE/GRACE-FO may not capture localized variations, its basin-scale sensitivity is well-suited to characterizing basin-wide groundwater trends. Moreover, the dominance of groundwater extraction in many deltas probably ensures that GWS trends are the primary signal captured.
We find a modest linear correlation (R = 0.5) between GWS and subsidence rate; however, a cubic regression model (R = 0.6) provides a better fit (Extended Data Fig. 8a).
The pristine and disturbed sediment flux, along with computed sediment flux changes for each delta, are summarized in Supplementary Table 2. A negative sediment flux change indicates a decline or loss in fluvial sediment supply (disturbed < pristine) due to human activities, whereas a positive sediment flux change reflects an increase or gain (disturbed > pristine). We acknowledge that this framework represents a simplified characterization of complex sediment delivery processes and may not capture all temporal variations in sediment supply. Furthermore, some concerns have been raised about potential errors in global sediment flux datasets, which we consider as a limitation in our analysis.
Figure 3a and Extended Data Fig. 8b show the relationship between sediment flux change and subsidence rates. Although a poor correlation (R < 0.4) is observed, we find that 62% of the deltas (25 out of 40) exhibit negative sediment flux change, indicating widespread human-induced reductions in sediment supply.
Supplementary Table 2 summarizes the urban fraction dataset (2000 and 2020) and the urban fraction change for each delta. Figure 3a and Extended Data Fig. 8c show the subsidence–urban expansion relationship across the 40 deltas. All deltas showed consistent urban expansion in the twenty-first century, ranging from relatively low increases (<1%) in the Ogooué river delta to significant increases (>400%) in the Indus delta. However, despite this rapid expansion, the Indus delta remains one of the least urbanized, with only 0.4% of its total area classified as urban in 2020. By contrast, the Ciliwung (Jakarta) and Neva (Saint Petersburg) deltas exhibit the highest urban fractions, exceeding 50%. A logarithmic fit best describes the full dataset and reveals a moderate but significant nonlinear inverse correlation (correlation, R = 0.38–0.51), indicating that deltas with significant urban land conversion tend to experience more pronounced land sinking (Extended Data Fig. 8c). Steadily urbanizing deltas, such as the Rio Grande and Rhine–Meuse, exhibit slower subsidence rates, whereas rapidly urbanizing deltas, such as the Brahmani and Yellow River deltas, show faster rates of land sinking. However, regional variability is evident, as some deltas deviate from the overall trend (for example, Indus and Cauvery deltas). When excluding outliers (the Indus and Cauvery deltas), subsidence and urban expansion exhibit a strong linear correlation across deltas (Extended Data Fig. 8c).
We also explored the relationship among the anthropogenic drivers (Extended Data Fig. 8d–f), finding a low (R = 0.1–0.3) correlation depending on the specific driver.
Given the nonlinear and interacting relationships among GWS, sediment flux alteration, urban expansion and residual land subsidence (after GIA correction) discussed above, a machine learning framework was implemented to model these complexities. First, we attempted a multilinear regression model, incorporating interaction terms between variables, formulated aswhere VLM is the predicted VLM, x0 is the intercept, X are the predictor variables (GWS, sediment flux alteration and urban expansion), x are the regression coefficients for each predictor variable, x represents the interaction effects between predictor variables and ϵ is the residual error term. However, this multilinear regression model yielded poor performance (correlation R = 0.38; R2 = 0.15; RMSE = 4.7 mm yr1) (Fig. 3a), demonstrating the inefficiency of linear models to capture these complex dependencies and the need for a machine learning model.
The RF regressor optimizes each decision tree using the mean square error (MSE) defined as a cost function to identify node splits and model performance during model training and testing:where VLM is the observed VLM rate for individual delta i, is the predicted VLM rate and N is the total number of observations. To assess uncertainty, we used Monte Carlo simulations to create multiple holdout fractions (0.1–0.5) across 100 iterations, randomly subsampling the 40 deltas for training and validation in each iteration. This random partitioning ensures that each delta is used in both training and validation phases across iterations, enhancing the robustness against overfitting and sampling bias. The final RF model predictions were obtained by averaging prediction estimates across all iterations. The final model performance was evaluated using the coefficient of determination (R2), RMSE and mean absolute error (MAE):where is the mean observed VLM rate, and the other variables are defined in equation (12). The feature importance If for input feature {X = X1, X2, X3} was computed using the following equation, based on the cumulative reduction in node, j impurity among all the trees:where N denotes the total number of trees and ΔI denotes the change in impurity.
Although RF effectively captures nonlinear relationships, its ensemble structure limits delta-specific interpretability. To resolve local insights into delta-specific subsidence drivers, we applied LIME, a technique within the field of explainable artificial intelligence (XAI). LIME approximates black-box models such as RF by fitting interpretable models to perturbed samples of the input data, allowing for local feature importance estimation. For each delta X, LIME approximates the RF prediction locally by using a linear surrogate model trained on perturbed instances around X. The explanation function is obtained by solving the following minimization problem:where ξ(X) is the local interpretable model for each delta X, g is the interpretable model, f is the RF model, is a proximity kernel, is the loss function measuring the differences between f and g, and Ω(g) penalizes complexity. This process was repeated for each delta, and deltas with low LIME model fidelity (R2 < 0.5) were excluded to ensure reliable interpretation (Supplementary Table 2). The final dataset for interpretation consisted of 30 deltas, in which LIME produced more consistent feature importance estimates. The feature importance scores from LIME are normalized to obtain normalized LIME (nLIME) scores:where ωf is the LIME-derived coefficient for feature f and F is set for all features. The nLIME scores provide an instance-specific (local) explanation rather than a global one to evaluate the relative contributions of GWS, sediment flux alteration and urban expansion in each delta. The nLIME values for each delta are summarized in Supplementary Table 2 and were analysed in a ternary diagram to visualize the heterogeneity in delta-specific subsidence and elevation-loss drivers (Fig. 3b).
Historical relative sea-level changes were obtained from the Revised Local Reference database of the Permanent Service for Mean Sea Level (https://psmsl.org), which provides monthly relative sea-level records from globally distributed tide gauge stations. These tide gauge records have undergone quality control procedures, including corrections for datum inconsistencies, jumps and spurious data points, and validation through comparisons with neighbouring tide gauge stations. For this study, we selected 20 tide gauge stations across 15 deltas (the Mississippi, Rio Grande, Fraser, Amazon, Chao Phraya, Mekong, Red River, Nile, Ganges–Brahmaputra, Vistula, Rhine–Meuse, Chikuma-gawa, Yangtze, Pearl and Rioni deltas), considering only stations within 100 m of the delta boundary and at least 5 years (twentieth century) of valid record. The RSLR rates for each delta were estimated by applying the stochastic-seasonal model (equations (6–8)) over the full observational record for each tide gauge. For deltas with multiple stations (for example, the Mississippi, Ganges–Brahmaputra and Rhine–Meuse deltas), individual station rates were averaged to provide a delta-wide estimate of twentieth-century RSLR. Note that the representativeness of the derived RSLR may vary for each delta following individual tide gauge characteristics (for example, is the station founded on bedrock or ‘floating’ in unconsolidated sediments, is the station GNSS corrected). Supplementary Fig. 7 shows the time series of relative sea level over the twentieth century for six representative deltas. Supplementary Table 1 provides a complete summary of the RSLR rates for the 15 deltas. The median twentieth-century RSLR trend across all deltas is 2.9 mm yr–1, with measured rates ranging from –0.5 mm yr–1 in the Amazon delta (indicating declining twentieth-century sea level) to a maximum rate of 1.5 cm yr–1 in the Chao Phraya Delta (Fig. 5b).
To estimate present-day (early twenty-first century) absolute (geocentric) SLR rates, we used the multi-mission satellite altimetry data from 2001 to present, obtained from Copernicus Marine Environment Monitoring Service (CMEMS). This dataset provides 1/8° (about 12.5 km) gridded monthly sea level anomalies (SLA) referenced to a 20-year mean baseline (1993–2012). SLA estimates are derived from optimal interpolation, merging the level 3 along-track measurement from multiple contemporaneous altimeter missions (Jason-3, Sentinel-3A, HY-2A, Saral/AltiKa, Cryosat-2, Jason-2, Jason-1, TOPEX/Poseidon, ENVISAT, GFO and ERS1/2) (https://marine.copernicus.eu/). Several necessary corrections have been applied to the raw altimetry data, including instrumental biases and drifts, geophysical, tidal and atmospheric corrections, to ensure accurate SLA estimates. Monthly mean sea-level anomalies were obtained for each delta by spatially averaging the altimetry grid points within a 100-m radius, culling outliers beyond the 95th percentile. Supplementary Fig. 8 shows the monthly SLA time series in six deltas. We estimated the twenty-first-century trends in sea-level anomalies, using equations (6–8). The altimetry-derived geocentric SLR rates for the twenty-first century show exacerbating regional SLR rates over global sea-level estimates (about 4 mm yr–1) for 45% of the deltas (18 out of 40) (Supplementary Table 1). Regional sea-level rates vary from 0.2 mm yr–1 in the Parana delta to 7.3 mm yr–1 over the Mississippi delta (Fig. 1 and Supplementary Table 1). However, a negative geocentric sea-level rate of −1.9 mm yr–1 was observed in the Rioni Delta (Black Sea) (Supplementary Table 1). This long-term sea-level decline in the twenty-first century persists in the background of short-term fluctuations (Supplementary Fig. 8d); a characteristic feature of Black Sea sea-level dynamics. This twenty-first-century decline in geocentric sea level for the Rioni Delta represents more than a 100% reduction compared with historical (twentieth-century) rates, even when accounting for average VLM across the delta. To investigate this anomaly, we estimated VLM at the Poti tide gauge (Rioni Delta) by differencing twenty-first-century RSLR rates obtained from the Poti tide gauge station from geocentric SLR. The resulting VLM rate of −6.7 mm yr–1 matches the average InSAR-derived VLM rate (−5.9 ± 0.7 mm yr–1) within 100 m of the tide gauge. This rapid subsidence rate at the coast of Poti represents localized conditions and highlights the need for caution when extrapolating point-based tide gauge measurements to infer delta-wide or city-wide subsidence and exposure. Note that satellite altimetry data, although highly valuable for global sea-level monitoring, were primarily optimized for open ocean conditions. Coastal environments naturally exhibit additional complexity due to processes such as shelf circulation, freshwater discharge and tidal amplification, which contribute to the inherent variability in nearshore sea-level measurements compared with offshore altimetric observations.
We use projected sea-level rates from the Intergovernmental Panel on Climate Change Sixth Assessment Report (AR6) to assess future SLR rates across all deltas. The sea-level rate projections integrate process-based models that account for the key contributors to climate-induced sea-level change, such as thermal expansion, ocean dynamics, and glacier and ice sheet mass loss, and consider uncertainties in global temperature change and their influence on sea-level drivers. We focus on the no-VLM 50th percentile (median) projected rates for 2050 (mid-twenty-first century) and 2100 (end of the twenty-first century) under shared socioeconomic pathway 2-4.5 (SSP2-4.5) and SSP5-8.5 scenario. SSP5-8.5 represents a high reference scenario associated with the highest emission levels (global atmospheric CO2 concentrations exceeding 800–1,100 ppm by 2100) and associated warming of 3.3–5.7 °C (refs. ). These projections provide an upper-bound reference scenario, capturing the potential worst-case outcome for future SLR. Figure 4c shows the comparison of projected SLR rates with observed land subsidence rates.

Sections

"[{\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig6\"], \"section\": \"Main\", \"text\": \"This recognized importance, which makes deltas indispensable, also increases their exposure to compounding climatic, environmental and anthropogenic threats. As low-lying landforms, with extensive areas less than 2\\u2009m above sea level, deltas are acutely susceptible to rising sea level, storm surge, land subsidence, shifting temperature and rainfall patterns, and other environmental pressures, which are amplified by climate change. These pressures degrade agricultural land; disrupt freshwater availability; exacerbate coastal and fluvial flooding; promote wetland loss, saltwater intrusion and shoreline retreat; and threaten infrastructure in deltas. Beyond direct physical impacts, the interplay of these hazards also creates potential cascading socioeconomic consequences. For example, land loss and freshwater scarcity may drive displacement and migration, heightening competition for dwindling resources and fuelling social tensions. Together, these intersecting climatic, environmental, human-driven pressures and multi-hazards render deltas the most fragile landscapes on Earth, with their low elevation and high urban exposure placing them at the forefront of climate and environmental risks (Extended Data Fig. 1).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig1\", \"Sec6\"], \"section\": \"Main\", \"text\": \"Here, we present high-spatial-resolution datasets of surface-elevation change derived from Sentinel-1 synthetic aperture radar (SAR) interferometry across 40 deltas globally (Fig. 1). These datasets capture delta-wide temporal trends, subsidence rates and horizontal motion at 75\\u2009m resolution, spanning five continents and 29 countries. Our analysis encompasses all major river deltas with a population exceeding 3 million people, historically recognized sinking deltas and representatives of less-populated, understudied deltas of regional ecological and economic importance (Methods).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM4\"], \"section\": \"Land subsidence in global deltas.\", \"text\": \"Source data\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Sec6\"], \"section\": \"Global analysis of delta subsidence\", \"text\": \"We measured the spatial patterns and rates of subsidence in 40 deltas by analysing the complete archive of the Sentinel-1 SAR dataset between 2014 and 2023 using advanced multitemporal interferometric SAR (InSAR) analysis (Methods). InSAR measures surface-elevation changes, capturing vertical land motion (VLM), sediment deposition and erosional processes. For consistency, to reflect both VLM and surface-elevation change in the deltas, we use the terms VLM or elevation gain or loss to describe net surface-elevation change across all delta environments, with positive values indicating uplift or elevation gain and negative values indicating subsidence or net elevation loss. Throughout this study, negative VLM is quoted with negative signs and references land subsidence rates, whereas only the absolute values are reported when presenting subsidence rates.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig1\", \"Fig1\", \"MOESM2\", \"Fig1\", \"MOESM2\"], \"section\": \"Global analysis of delta subsidence\", \"text\": \"Our analysis shows that subsidence threatens deltas globally, with the delta-scale average rate of VLM on all deltas indicating subsidence (Fig. 1). In 12 out of 40 deltas, the average sinking rate is moderate, at less than 2\\u2009mm\\u2009yr\\u22121. By contrast, more than half of the deltas exhibit subsidence rates exceeding 3\\u2009mm\\u2009yr\\u22121, and in 13 of these deltas (Nile, Po, Vistula, Ceyhan, Brahmani, Mahanadi, Chao Phraya, Mekong, Red, Ciliwung, Brantas, Godavari and Yellow River), the average subsidence rates exceed the current estimates of global SLR (that is, about 4\\u2009mm\\u2009yr\\u22121). Among these, the Chao Phraya (Thailand), Brantas (Indonesia) and Yellow River (China) deltas show an average sinking rate of more than twice the current global SLR rate. To further highlight the severity of subsidence in deltas, we compared the subsidence with the regional geocentric SLR rates for the twenty-first century (2001\\u2013present). In 18 of the 40 deltas (the Nile, Po, Vistula, Ceyhan, Rioni, Brahmani, Mahanadi, Ganges\\u2013Brahmaputra, Godavari, Chao Phraya, Mekong, Red River, Ciliwung, Brantas, Amazon, Parana, Pearl and Yellow River), the average rate of local land subsidence is greater than the rate of regional geocentric SLR (Fig. 1 and Supplementary Table 1). However, in almost every delta (except Rio Grande) at least 1% of the delta area is subsiding faster than both global and geocentric sea levels (Fig. 1 and Supplementary Table 1).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig2\", \"Fig2\", \"MOESM2\"], \"section\": \"Global analysis of delta subsidence\", \"text\": \"Among all deltas, we find that at least 35% of the area is sinking, and in 38 deltas (excluding Neva and Fraser), more than 50% of the delta area is sinking (Fig. 2a). Of the 40 deltas, 19 show widespread subsidence patterns, with greater than 90% of the delta area affected by subsidence (for example, Mississippi, Niger, Nile, Rhine\\u2013Meuse, Po, Vistula, Brahmani, Mahanadi, Ganges\\u2013Brahmaputra, Chao Phraya, Mekong and Brantas deltas). Deltas with notable subsiding areas with greater than 50% of the delta area sinking faster than 5\\u2009mm\\u2009yr\\u22121 include the Chao Phraya (94%\\u00a0of delta area), Nile (80%), Brahmani (77%), Po (74%), Mahanadi (69%), Brantas (66%), Vistula (57%), Yellow River (53%) and Mekong (51%) deltas (Fig. 2a and Supplementary Table 1). In sum, we estimate that a total delta area of 460,370\\u2009km2 is exposed to subsidence. If we consider a global habitable geomorphic area of 710,000\\u2013855,000\\u2009km2 for deltas, approximately 54\\u201365% of global delta areas are sinking just from the analysis of the 40 deltas. By region, South Asia, East Asia and Southeast Asia, with 17 representative deltas, have the greatest exposure to subsidence, with 274,000\\u2009km2 of delta area subsiding. Africa, South America, North America and Europe have total subsiding delta areas of 78,800\\u2009km2, 39,800\\u2009km2, 37,800\\u2009km2, and 30,000\\u2009km2, respectively. Seven large deltas\\u2014Ganges\\u2013Brahmaputra, Nile, Mekong, Yangtze, Amazon, Irrawaddy and Mississippi deltas\\u2014contribute about 57% of the total subsiding delta area, with a combined area of 265,000\\u2009km2. Coastal cities such as Alexandria (Nile), Bangkok (Chao Phraya), Dhaka and Kolkata (Ganges\\u2013Brahmaputra), Shanghai (Yangtze), Yangon (Irrawaddy), C\\u1ea7n Th\\u00e1 (Mekong), Th\\u00e1i B\\u00ecnh (Red River), Niigata (Chikuma-gawa), Jakarta (Ciliwung), Surabaya (Brantas) and Dongying (Yellow River) are experiencing subsidence at rates equal to or exceeding the delta-wide averages, indicative of the intensity of subsidence and elevation loss processes in cities on deltas.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig7\", \"Fig9\"], \"section\": \"Spatial pattern of VLM across global deltas.\", \"text\": \"a, Proportion of each delta exposed to different rates of subsidence. Note that only subsiding areas are represented in each bar, and areas of uplift within each delta are omitted to emphasize the extent of elevation loss. b\\u2013m, Spatial maps of VLM rates for the Fraser (Canada) (b), Mississippi (the USA) (c), Parana (Argentina) (d), Niger (Nigeria) (e), Nile (Egypt) (f), Po (Italy) (g), Ganges\\u2013Brahmaputra (India\\u2013Bangladesh) (h), Chao Phraya (Thailand) (i), Mekong (Vietnam) (j), Red River (Vietnam) (k), Pearl (China) (l), Yellow River (China) (m)\\u00a0deltas. Positive VLM (green\\u2013purple hues) suggests uplift or elevation gain, whereas negative VLM (yellow\\u2013orange\\u2013red hues) indicates land subsidence. The spatial\\u00a0VLM maps for the other 28 deltas are shown in Extended Data Figs. 2\\u20134. Background image in b\\u2013m is Esri, streets-dark. Scale bars, 5\\u2009km (b); 50\\u2009km (c,f,h,j); 20\\u2009km (d,e,i,k,l,m); 10\\u2009km (g).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM4\"], \"section\": \"Spatial pattern of VLM across global deltas.\", \"text\": \"Source data\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig2\", \"Fig7\", \"Fig9\", \"Fig2\", \"Fig7\", \"Fig8\", \"Fig10\", \"Fig12\"], \"section\": \"Global analysis of delta subsidence\", \"text\": \"Furthermore, we observe non-uniform spatially variable VLM within individual deltas, reflecting the complex interplay of natural and anthropogenic processes (Fig. 2 and Extended Data Figs. 2\\u20134). Although all deltas exhibit an overall trend of subsidence, localized and broad zones of uplift, which vary from 0\\u2009mm\\u2009yr\\u22121 to greater than 5\\u2009mm\\u2009yr\\u22121 are observed in some areas (Fig. 2b,d,k,m, and Extended Data Figs. 2e,f,i,j,l and 3c,f). In some deltas (for example, Wouri, Zambezi, Indus, Ciliwung and Yellow River), the observed uplift or elevation-gaining parts correlate with patterns of horizontal land motion (Extended Data Figs. 5\\u20137). Possible mechanisms may include sediment redistribution processes potentially driven by river dynamics or growth faulting, either of which can cause localized zones of elevation gain even within a predominantly subsiding deltaic system. This highlights the necessity of comprehensive assessments and models of delta vulnerability to consider not only overall absolute subsidence rates but also the spatial heterogeneity of elevation change dynamics.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Sec6\", \"MOESM2\"], \"section\": \"Anthropogenic drivers of delta subsidence\", \"text\": \"To quantify the relative contributions of anthropogenic factors to delta subsidence and elevation loss, we analysed the relationship between three main anthropogenic drivers\\u2014groundwater storage change, sediment flux alteration and urban expansion\\u2014and non-glacial isostatic adjustment VLM/subsidence rates across the 40 deltas (Methods and Supplementary Table 2).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig3\", \"Fig13\"], \"section\": \"Anthropogenic drivers of delta subsidence\", \"text\": \"Figure 3a and Extended Data Fig. 8 show the interplay of anthropogenic factors and their correlation with subsidence rates across the 40 deltas. Deltas experiencing groundwater storage (GWS) loss (indicative of groundwater extraction), negative sediment flux change (red and yellow hues; reflecting sediment reduction due to upstream human activities) and higher urban population growth tend to have higher rates of subsidence (for example, the Yellow River, Po, Nile, Chao Phraya and Mekong deltas). Conversely, deltas with GWS stability or gain (net increase in groundwater storage), positive sediment flux change (blue colours; sediment surplus) and limited urban expansion show lower subsidence rates (for example, Saloum, Amazon and Ogoou\\u00e9 deltas).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig13\"], \"section\": \"Anthropogenic drivers of land subsidence and elevation loss in global deltas.\", \"text\": \"a, Bubble plot showing the relationship between VLM rates and anthropogenic drivers across deltas. Plot shows VLM rate (mm\\u2009yr\\u22121) against GWS rate (mm\\u2009yr\\u22121). Bubble colours represent sediment flux change (%), in which positive values (blue colours) indicate increased sediment supply due to human activities (and thereby increased potential to gain elevation and compensate subsidence-induced elevation loss), whereas negative values (yellow\\u2013orange\\u2013red colours) indicate a decline in sediment availability. Bubble size indicates urban fraction change (%), with larger circles representing a greater urban expansion over the twenty-first century. The dashed line represents the MLR fit. See Extended Data Fig. 8 for individual pairwise relationships between each anthropogenic driver and VLM. b, Ternary plot of subsidence rates with nLIME scores.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM4\"], \"section\": \"Anthropogenic drivers of land subsidence and elevation loss in global deltas.\", \"text\": \"Source data\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig3\"], \"section\": \"Anthropogenic drivers of delta subsidence\", \"text\": \"The initial multilinear regression (MLR) model, which included interaction terms between the different anthropogenic factors, poorly captured subsidence dynamics on the deltas (R2\\u2009=\\u20090.2\\u2009\\u00b1\\u20090.1), as it failed to account for nonlinear interactions between the different processes (Fig. 3a). For instance, urban expansion not only directly increases infrastructure loading but also indirectly elevates groundwater demand, thereby compounding aquifer depletion and extraction-induced subsidence, which are synergistic effects that linear models cannot resolve.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig3\", \"MOESM1\", \"MOESM1\"], \"section\": \"Anthropogenic drivers of delta subsidence\", \"text\": \"To address these limitations, we used a random forest (RF) machine learning approach designed to capture nonlinear relationships and variable interactions. The RF model shows a moderate to strong relationship between the predictors (GWS, sediment flux and urban expansion) and VLM, achieving improved performance over the MLR model (R2\\u2009=\\u20090.6\\u2009\\u00b1\\u20090.1; RMSE (root mean square error)\\u2009=\\u20091.9\\u2009\\u00b1\\u20090.1\\u2009mm\\u2009yr\\u22121; MAE (mean absolute error)\\u2009=\\u20091.4\\u2009\\u00b1\\u20090.2\\u2009mm\\u2009yr\\u22121), and capturing complex, non-additive relationships between anthropogenic stressors and subsidence rates (Fig. 3a and Supplementary Fig. 1). However, we observe some underestimation at high subsidence rates (>8.0\\u2009mm\\u2009yr\\u22121) (Supplementary Fig. 1), which probably suggests that natural processes or other anthropogenic predictors (not considered in our analysis) may contribute to subsidence in these highly dynamic deltaic environments.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig3\", \"MOESM1\", \"Sec6\", \"MOESM2\", \"MOESM1\"], \"section\": \"Anthropogenic drivers of delta subsidence\", \"text\": \"Note that the primary objective in our analysis is not to predict subsidence rates across deltas, but rather to identify and extract key features that explain the dynamic relationships between the three anthropogenic drivers and subsidence across these deltas. Feature importance analysis from the RF model identifies GWS as the dominant anthropogenic predictor of delta subsidence (0.5\\u2009\\u00b1\\u20090.2), whereas sediment flux change (0.3\\u2009\\u00b1\\u20090.2) and urbanization (0.3\\u2009\\u00b1\\u20090.1) have secondary roles as subsidence rate predictors across these deltas (Fig. 3a and Supplementary Fig. 1b). However, the large standard deviations in feature importance values reflect substantial variability in predictor dominance across subsampled delta subsets, suggesting that the primary contributors to subsidence differ locally depending on the anthropogenic or geomorphic context. To resolve delta-specific mechanisms, we applied local interpretable model-agnostic explanations (LIME), which interprets individual predictions by approximating the RF model locally with simpler, interpretable functions. Deltas with low LIME model fidelity (R2\\u2009<\\u20090.5) were excluded from this interpretative analysis, refining the dataset from 40 to 28 deltas (Methods). The low fidelity scores for some deltas could be due to unaccounted processes (natural and/or other anthropogenic) in our RF model. The retained 28 deltas show improved overall model performance (R2\\u2009=\\u20090.7\\u2009\\u00b1\\u20090.1; RMSE\\u2009=\\u20090.4\\u2009\\u00b1\\u20090.1\\u2009mm\\u2009yr\\u22121; MAE\\u2009=\\u20090.3\\u2009mm\\u2009yr\\u22121), ensuring reliable interpretation of local feature importance. Normalized LIME feature importance scores (nLIME) showed substantial heterogeneity in predictor dominance (Supplementary Table 2). GWS emerged as the most significant factor across the different deltas (0.6\\u2009\\u00b1\\u20090.3), whereas sediment flux change (0.3\\u2009\\u00b1\\u20090.1) and urbanization (0.1\\u2009\\u00b1\\u20090.1) exhibited lower but context-dependent impacts (Supplementary Fig. 1b).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig3\", \"Fig3\", \"MOESM2\"], \"section\": \"Anthropogenic drivers of delta subsidence\", \"text\": \"To assess the dominant influence on land motion across individual deltas, the nLIME for each delta was mapped onto a ternary diagram (Fig. 3b). Of the 28 deltas, 35%, including the Mekong, Ganges\\u2013Brahmaputra, Rhine\\u2013Meuse, Fraser, Cauvery, Irrawaddy and Red River systems, cluster within the GWS portion of the diagram (nLIMEGWS\\u2009>\\u20090.7), suggesting that observed GWS changes in these deltas are the primary driver of subsidence among the three anthropogenic variables examined (Fig. 3b and Supplementary Table 2). The Chao Phraya and Yellow River deltas, with the highest average subsidence rates, plot near the centre of the ternary diagram, reflecting relatively balanced contributions from GWS, sediment flux and urban expansion. Sediment flux correlates most closely with elevation changes in deltaic systems, such as the Saloum, Mississippi, Amazon and Rio Grande deltas, suggesting that reduced sediment delivery may exacerbate land elevation loss in these deltas. The Nile, Po, Chikuma-gawa, Mahanadi, Kabani, Niger and Volta deltas exhibit mixed contributions from GWS, sediment flux changes and population change, with GWS slightly outweighing sediment deficits as predictors in the Nile and Po deltas, possibly reflecting reliance on aquifer-dependent irrigation. These findings are consistent with delta-specific studies that attribute accelerated subsidence in densely populated Asian deltas\\u2014Mekong, Ganges\\u2013Brahmaputra and Chao Phraya\\u2014to urbanization and unsustainable groundwater extraction for agriculture, industry and domestic use. Moreover, the Nile, Po and Mississippi deltas, which were historically sustained by seasonal floods that deposited sediments, are now documented to experience severe sediment deficits due to dams and levees, accelerating elevation loss.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig4\", \"Fig4\", \"Fig4\", \"Fig4\", \"MOESM2\"], \"section\": \"Relative impacts of SLR and subsidence\", \"text\": \"To quantify the contributions of SLR and land subsidence in deltas, we evaluated their relative impact on the exposed delta populations. Our analysis shows that current average subsidence rates exceed geocentric SLR in 18 of the 40 deltas, including the Nile, Mekong, Red River, Ganges\\u2013Brahmaputra, Brahmani, Mahanadi, Chao Phraya, Ciliwung, Brantas and Yellow River deltas, affecting approximately 236 million people\\u2014a population about 50% larger than those residing in deltas in which the current rates of geocentric SLR outpace the subsidence rates (156.9 million) (Fig. 4a). This disparity is particularly pronounced for vulnerable populations occupying land below 1\\u2009m elevation. In these lowest elevation areas, subsidence dominates the contribution to RSLR in about two-thirds of the deltas, including Amazon, Fraser, Niger, Rhone, Vistula, Ganges\\u2013Brahmaputra, Mekong, Red River, Pearl, Yangtze and Godavari deltas (Fig. 4b). Of the 76 million people living in delta areas with an elevation below 1\\u2009m, 84% (63.7 million people) reside in rapidly sinking areas of the deltas (Fig. 4b). These observations are striking, revealing the current dominance of subsidence over geocentric SLR in global deltas. Moreover, the spatial heterogeneity of VLM creates localized extreme rates of subsidence within deltas, further exacerbating their vulnerability. Under the current trajectory, moderate emission scenarios (shared socioeconomic pathway 2-4.5 (SSP2-4.5)), current maximum subsidence rates in the deltas already surpass projected twenty-first-century SLR rates (no VLM). Through the end of the twenty-first century, current maximum subsidence rates in all 40 deltas exceed projected SLR rates (Fig. 4c). This disparity extends to the 95th percentile subsidence rates, representing widespread, high-magnitude sinking across the deltas. In 29 deltas, 95th percentile subsidence rates exceed the projected SLR rates by 2050, outpacing SLR by 1.1 (Niger delta) to 10.3 (Yellow River delta) times. By 2100, as the current maximum rate of SLR (SSP2-4.5) accelerates to 0.9\\u2009cm\\u2009yr\\u22121, current 95th percentile subsidence rates still dominate in 22 deltas, surpassing geocentric SLR by up to seven times. Even accounting for worst-case, high-emission scenarios (SSP5-8.5), subsidence will exceed projected SLR rates in all deltas (considering maximum subsidence) and in 23 deltas (considering 95th percentile subsidence) through 2050. By 2100, current maximum subsidence rates exceed projected SLR in 38 of 40 deltas, whereas 95th percentile subsidence rates remain dominant in seven deltas (Godavari, Chao Phraya, Mekong, Ciliwung, Brantas, Red River and Yellow River) (Supplementary Table 1).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM4\"], \"section\": \"Relative contributions of land subsidence and SLR in global deltas.\", \"text\": \"Source data\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig5\", \"Fig5\"], \"section\": \"Adaptive capacity in vulnerable deltas\", \"text\": \"To visualize disparities in adaptive capacity and risk, we mapped global deltas into a two-dimensional (2D) impact matrix defined by RSLR and ND-GAIN adaptation readiness scores (Fig. 5). This framework allows for a comparative assessment of deltas assuming that the adaptation readiness of the delta is reflected by the adaptation readiness of its country, categorizing them into four quadrants: (1) Unprepared Divers (high RSLR (>4\\u2009mm\\u2009yr\\u22121), low readiness (<0.52)); (2) Rising Ready (high RSLR, high readiness (>0.52)); (3) Latent Threats (low RSLR, low readiness); and (4) Safe Havens (low RSLR, high readiness). 65% of the deltas (26 out of 40 deltas), predominantly in low- and middle-income nations, fall into the Unprepared Divers group, in which nations have a diminished adaptive capacity and RSLR rates exceeding current global SLR (Fig. 5a). These challenges are compounded for indigenous communities, who primarily live in the lowest-lying delta areas; lack the resources needed to implement large-scale adaptation; and face relocation barriers due to cultural and subsistence ties despite escalating risks.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM4\"], \"section\": \"RSLR and adaptive capacity in global deltas.\", \"text\": \"Source data\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig5\", \"Fig5\", \"Fig4\", \"Sec6\"], \"section\": \"Adaptive capacity in vulnerable deltas\", \"text\": \"Most deltas in high-income countries, including the Yellow River (China), Vistula (Poland), Po (Italy), Rhine\\u2013Meuse (the Netherlands) and Mississippi (the USA) deltas, cluster in the Rising Ready group, demonstrating robust governance (Fig. 5a). For example, the integrated flood management approach of the Dutch delta, which combines ecological restoration with infrastructural fortifications, has become a model for coastal hazard resilience. However, some deltas even within this group face substantial gaps. For instance, the Mississippi delta has lost more than 5,000\\u2009km2 of land (mainly wetlands) since 1932 because of a lack of adaptation (for example, sediment diversion projects), whereas the Po delta struggles with salinization driven by agricultural groundwater extraction, highlighting how economic priorities can undermine adaptation even in high-income regions. Although RSLR exceeds global rates of SLR in most deltas, exceptions exist. The Latent Threats group includes the Saloum and Neva deltas, which exhibit relatively low RSLR and low adaptive capacity (Fig. 5a), indicating their unpreparedness and potential vulnerability to a future rise in sea level (Fig. 4c). The Rioni and the Fraser delta fall into the Safe Havens group, in which lower RSLR is coupled with higher adaptive capacity, indicative of low risk and preparedness for current and future sea-level changes. The Rioni Delta is the only delta in our sample exhibiting negative sea-level trends for the twenty-first century, in which long-term regional sea-level decline masks short-term fluctuations (Methods).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig5\", \"Fig5\", \"Sec6\", \"Fig5\", \"Fig5\"], \"section\": \"Adaptive capacity in vulnerable deltas\", \"text\": \"To examine the evolving risk landscape, we compared twentieth-century and present-day impact matrices (Fig. 5b). For our analysis, we used tide gauge data to estimate twentieth-century RSLR rates, which were available for only 15 of the 40 deltas. Our estimates show that 10 deltas previously classified as Latent Threats (low RSLR, low readiness) and Safe Havens (low RSLR, high readiness) groups during the twentieth century have transitioned to Unprepared Divers (high RSLR, low readiness) and Rising Ready (high RSLR, high readiness) groups in the twenty-first century (Fig. 5b). This shift highlights the accelerating contemporary RSLR trends, driven by land subsidence and SLR. Deltas such as the Mississippi, Ganges\\u2013Brahmaputra and Mekong show sustained increases in long-term RSLR rates above 4.0\\u2009mm\\u2009yr\\u22121 since the twentieth century, exacerbating vulnerabilities in these densely populated regions. Conversely, the Chao Phraya and the Rioni deltas showed a decline in RSLR and improved adaptive capacity in the twentieth century. However, although the Rioni Delta exhibited a more than 200% decline in RSLR, the Chao Phraya Delta still experiences high RSLR rates (12.3\\u2009mm\\u2009yr\\u20131). The pronounced decrease in RSLR for the twentieth century in the Rioni Delta probably reflects localized subsidence at the tide gauge station rather than a delta-wide RSLR trend (Methods). The greatest change in RSLR was observed in the Nile Delta, surging from 1\\u2009mm\\u2009yr\\u20131 in the twentieth century to more than 10\\u2009mm\\u2009yr\\u20131 in the twenty-first century (Fig. 5b). Moreover, we find that all deltas in low- and middle-income countries in the present-day Unprepared Divers groups, transitioned from the Latent Threats group, suggesting stagnant adaptive capacity despite worsening RSLR. By contrast, deltas such as the Yangtze (China), Pearl (China) and Vistula (Poland) shifted from Latent Threats to Rising Ready, demonstrating increased adaptation readiness due to economic growth, raising governance and institutional capacity to adapt, although RSLR has surged (Fig. 5b). Although deltas in the Rising Ready quadrant showed potential for robust adaptation policies, deltas in the Unprepared Divers remain trapped in cycles of reactive, underfunded responses.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig6\", \"Fig6\"], \"section\": \"Selection of global river deltas\", \"text\": \"Our estimates show that globally, 42,000\\u2009km2 of the delta area at present lies below sea level, containing a population of 10.2 million people (Extended Data Fig. 1). The 35 deltas with the greatest exposure included in this analysis are Nile, Mississippi, Rhine\\u2013Meuse, Mekong, Niger, Cauvery, Po, Red River, Vistula, Rhone, Amazon, Ganges\\u2013Brahmaputra, Chao Phraya, Kabani, Pearl, Rio Grande, Yangtze, Yellow River, Senegal, Indus, Saloum, Grijalva, Ceyhan/Seyhan, Rioni, Cross, Chikuma-gawa, Volta, Brantas, Neva, Wouri, Irrawaddy, Ogoou\\u00e9, Zambezi, Magdalena and Ciliwung (Extended Data Fig. 1). The cumulative delta area and population below sea level are 38,000\\u2009km2 and 10.1 million people, respectively, reaching within rounding errors of the global total exposure. Deltas such as the Danube, Orinoco and Shatt-el-Arab met the selection criteria but were excluded due to challenges associated with the SAR imaging and interferometric analysis (including spatial coverage gaps, excessive temporal baselines, poor coherence and limited data availability). The five supplementary deltas are Brahmani, Mahanadi, Godavari, Parana and Fraser deltas.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM2\"], \"section\": \"SAR dataset\", \"text\": \"We analysed 132 SAR frames from the Sentinel-1A/B C-band satellite, spanning September 2016 to May 2023. The SAR datasets include 3,300 images obtained in single-orbit geometry (ascending or descending) for 13 deltas and 10,700 images obtained in both ascending and descending orbits for 27 deltas. See Supplementary Table 3 for the complete inventory of SAR images used in each delta. For each SAR dataset, we applied a multi-looking factor of 32:6 (range:azimuth) to improve the signal-to-noise ratio, obtaining an average pixel resolution of about 75\\u2009m. To minimize decorrelation errors, we also constrained the interferometric pairs to a maximum temporal and perpendicular baselines of 300\\u2009days and 80\\u2009m, respectively. For deltas requiring multi-frame coverage (for example, Amazon, Mississippi, Mekong, Ganges\\u2013Brahmaputra, Nile, Red River and Niger), we arranged in a mosaic form the overlapping adjacent frames along a single path before processing or post-processed deltas with coverage spanning multiple paths to ensure full spatial continuity across expansive deltas.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Equ1\"], \"section\": \"SAR interferometric analysis\", \"text\": \"In the 27 deltas with overlapping spatiotemporal SAR satellite coverage and different orbit geometries (ascending and descending), we estimate the horizontal (east\\u2013west) and VLM components of deformation by jointly inverting the LOS time series of the ascending and descending tracks. To this end, we identified the co-located pixels of the LOS time series by resampling the pixels from the descending track onto the ascending track to obtain two co-located LOS displacement velocities . Given  and their associated variances  are the LOS displacement and variances for a given pixel, the model to combine the LOS velocities to generate a high-resolution map of the east\\u2013west (E) and VLM (U) displacements are given bywhere, C represents the unit vectors for projecting (E) and (U) displacements onto the LOS, which is a function of the heading angle of the satellite and incidence angles of each pixel. The solution to the model in equation (1) is given bywhere X represents the unknowns (E) and (U), G is the design matrix comprising the unit vectors for projecting the horizontal and vertical displacements onto the line of sight, L are the observations , and P is the weight matrix, which is inversely proportional to the observant variances . To obtain the parameter variance\\u2013covariance matrix (QXX), we use the concept of error propagation to calculate the associated parameter uncertainties given the observation errors as follows:\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM1\", \"MOESM1\", \"MOESM1\", \"MOESM1\"], \"section\": \"SAR interferometric analysis\", \"text\": \"The distribution of the standard deviations (precision of the results) for all pixels (20.5 million) across the 40 deltas is shown in Supplementary Fig. 2. The standard deviation distribution shows that 99% of the pixels have a value <0.5\\u2009mm\\u2009yr\\u20131. We evaluated the accuracy of the results by comparing the averaged VLM rates of pixels within a radius of 100\\u2009m with more than 100 independent GNSS data (that is, stations that were not used in the reference frame transformation). The validation included 122 GNSS stations across 23 deltas with historical long-term records (spanning various periods before and/or including the InSAR observation window) and 81 GNSS stations across 15 deltas with time series covering at least 70% of the InSAR observation period (2014\\u20132023) (Supplementary Fig. 3). We found a strong correlation (0.7\\u20130.8), between GNSS and InSAR velocities, with an RMSE of 1.4\\u2009mm\\u2009yr\\u20131 for long-term rates (Supplementary Fig. 3a) and 1.2\\u2009mm\\u2009yr\\u20131 for rates within the InSAR observation period (Supplementary Fig. 3b). The improved agreement for temporally coincident measurements suggests that nonlinear subsidence behaviour contributes to some scatter when comparing historical GNSS rates to contemporary InSAR measurements, although the overall correlation remains strong in both cases. Note that some GNSS stations used for validation, while within the broader processed SAR frame, are outside the clipped delta boundaries. Note that the final delta extents were delineated using a tiered approach. Primary boundaries were derived from ref.\\u2009, supplemented by ref.\\u2009 for deltas not covered in the former. For extensive deltas in which the entire delta surface is not analysed (for example, the Ganges\\u2013Brahmaputra), boundaries were defined using the SAR spatial extent.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM2\"], \"section\": \"GIA influence on VLM\", \"text\": \"For each delta, we extracted the ensemble mean and standard deviation in GIA vertical velocity to correct observed deformation rates and isolate contemporary, non-GIA contributions to VLM. Supplementary Table 1 shows the mean GIA-induced VLM, the associated standard deviation and the per cent contribution of GIA to the total observed VLM magnitude for each delta. GIA accounts for the largest proportion and exceeds (>100%) the total VLM in the Neva (540%) and Fraser (455%) deltas, in which low observed VLM rates are substantially influenced by strong GIA uplift. Moderate GIA contributions (25\\u201355%) are observed in five deltas, including the Rio Grande, Mississippi, Volta, Rhine and Ogoou\\u00e9 deltas. Most of the deltas (55%) exhibit minimal GIA influence, with contributions under 10%, indicating that observed VLM is primarily governed by contemporary anthropogenic and natural processes such as groundwater withdrawal, sediment compaction, or tectonics. In 28\\u201367% (accounting for uncertainty) of the deltas, the sign and approximate magnitude of observed and GIA-corrected VLM are consistent, implying limited distortion from GIA and the sustained expression of contemporary processes on the average local subsidence. By contrast, the Fraser and Neva deltas illustrate how substantial GIA-induced uplift in high-latitude, post-glacial regions can obscure contemporary subsidence processes through opposing vertical trends. In both cases, modest observed subsidence rates (Fraser \\u22120.4\\u2009mm\\u2009yr\\u22121 and Neva \\u22120.2\\u2009mm\\u2009yr\\u20131) are counteracted by substantial GIA uplift of 1.8\\u2009\\u00b1\\u20092.3\\u2009mm\\u2009yr\\u22121 and 1.0\\u2009\\u00b1\\u20090.3\\u2009mm\\u2009yr\\u22121, respectively.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM2\"], \"section\": \"Anthropogenic drivers datasets\", \"text\": \"We analysed the relationship between major anthropogenic pressures on global deltas to subsidence and elevation loss by quantifying the contributions of groundwater storage change, sediment flux alteration and urban expansion to the residual rates of sinking (after GIA correction) across the 40 deltas. These globally consistent datasets provide insights into human-induced impacts on land subsidence and elevation change in river deltas (Supplementary Table 2).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig1\"], \"section\": \"Groundwater storage change\", \"text\": \"To isolate GWS change from TWS, we used the 1/4\\u00b0 global land data assimilation system Noah model to remove changes in SMS and SWE contributions and used the WaterGAP Global Hydrology Model (WGHM v.2.2d) (refs.\\u2009) to remove SWS contributions. The contribution from SWE was negligible in most deltas, given their prevailing arid and semi-arid climate (Fig. 1), although it was included to maintain consistency across all deltas. SWS components include contributions from rivers, lakes, wetlands and reservoir storage within the GRACE footprint for each delta. The residual signal following removal of SWS, SMS and SWE was interpreted as the GWS anomaly.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM1\", \"MOESM1\", \"MOESM1\", \"MOESM1\", \"MOESM1\", \"MOESM1\"], \"section\": \"Groundwater storage change\", \"text\": \"Supplementary Figs. 4 and 5 compare the time-invariant model (black curves) with the stochastic-seasonal model (red curves) for GRACE-derived GWS and RSLR from tide gauges in the Mississippi and Chao Phraya deltas. These plots show that a stochastic-seasonal process better represents the observed variability in the time series. The post-fit residuals of the time-invariant model show some systematic seasonal patterns, particularly during periods when seasonal amplitudes deviate from the assumed constant values (Supplementary Figs. 4b,d and 5b,d). By contrast, the stochastic model accommodates time-dependent variations in seasonal amplitudes, resulting in reduced (often near-zero) residuals (Supplementary Figs. 4b,d and 5b,d), demonstrating the advantage of the stochastic-seasonal model in capturing transient seasonal variations rather than fixed annual and semiannual cycles.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM2\", \"Fig3\", \"Fig13\", \"MOESM1\"], \"section\": \"Groundwater storage change\", \"text\": \"The GWS rates for each delta are summarized in Supplementary Table 2, and Fig. 3a and Extended Data Fig. 8a show the relationship with the subsidence rates. Negative GWS trends indicate mass depletion, primarily driven by groundwater extraction, whereas positive trends represent net groundwater accumulation due to recharge processes, reduced extraction or hydrological interventions. To evaluate the reliability of GRACE-derived GWS trends, we compared them with in situ groundwater level trends for 18 deltas (Supplementary Fig. 6). Groundwater levels were compiled from two publicly available sources: 13 deltas from ref.\\u2009 and 5 deltas from the Global Groundwater Monitoring Network. Given the spatial scale discrepancy between GRACE (basin-wide) and well observations (point-scale), we emphasized agreement in trend direction rather than absolute magnitudes. Each site was categorized based on the sign of the GRACE and well trends, and a confusion matrix was constructed to assess consistency. The analysis yielded an overall classification accuracy of 88.9%, with six sites exhibiting positive\\u2013positive trends (PPT) and 10 showing negative\\u2013negative trends (NNT). Only two sites showed mixed behaviour (NPT or PNT), and no site exhibited fully opposing trends. Moreover, a high correlation (R\\u2009=\\u20090.7) was observed between the GRACE-based GWS and well-derived trends, further supporting the consistency of GRACE estimates at the basin scale despite localized variability in in situ measurements. Although the coarse spatial resolution of GRACE/GRACE-FO may not capture localized variations, its basin-scale sensitivity is well-suited to characterizing basin-wide groundwater trends. Moreover, the dominance of groundwater extraction in many deltas probably ensures that GWS trends are the primary signal captured.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig13\"], \"section\": \"Groundwater storage change\", \"text\": \"We find a modest linear correlation (R\\u2009=\\u20090.5) between GWS and subsidence rate; however, a cubic regression model (R\\u2009=\\u20090.6) provides a better fit (Extended Data Fig. 8a).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM2\"], \"section\": \"Sediment flux alteration\", \"text\": \"The pristine and disturbed sediment flux, along with computed sediment flux changes for each delta, are summarized in Supplementary Table 2. A negative sediment flux change indicates a decline or loss in fluvial sediment supply (disturbed\\u2009<\\u2009pristine) due to human activities, whereas a positive sediment flux change reflects an increase or gain (disturbed\\u2009>\\u2009pristine). We acknowledge that this framework represents a simplified characterization of complex sediment delivery processes and may not capture all temporal variations in sediment supply. Furthermore, some concerns have been raised about potential errors in global sediment flux datasets, which we consider as a limitation in our analysis.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig3\", \"Fig13\"], \"section\": \"Sediment flux alteration\", \"text\": \"Figure 3a and Extended Data Fig. 8b show the relationship between sediment flux change and subsidence rates. Although a poor correlation (R\\u2009<\\u20090.4) is observed, we find that 62% of the deltas (25 out of 40) exhibit negative sediment flux change, indicating widespread human-induced reductions in sediment supply.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM2\", \"Fig3\", \"Fig13\", \"Fig13\", \"Fig13\"], \"section\": \"Urban expansion\", \"text\": \"Supplementary Table 2 summarizes the urban fraction dataset (2000 and 2020) and the urban fraction change for each delta. Figure 3a and Extended Data Fig. 8c show the subsidence\\u2013urban expansion relationship across the 40 deltas. All deltas showed consistent urban expansion in the twenty-first century, ranging from relatively low increases (<1%) in the Ogoou\\u00e9 river delta to significant increases (>400%) in the Indus delta. However, despite this rapid expansion, the Indus delta remains one of the least urbanized, with only 0.4% of its total area classified as urban in 2020. By contrast, the Ciliwung (Jakarta) and Neva (Saint Petersburg) deltas exhibit the highest urban fractions, exceeding 50%. A logarithmic fit best describes the full dataset and reveals a moderate but significant nonlinear inverse correlation (correlation, R\\u2009=\\u20090.38\\u20130.51), indicating that deltas with significant urban land conversion tend to experience more pronounced land sinking (Extended Data Fig. 8c). Steadily urbanizing deltas, such as the Rio Grande and Rhine\\u2013Meuse, exhibit slower subsidence rates, whereas rapidly urbanizing deltas, such as the Brahmani and Yellow River deltas, show faster rates of land sinking. However, regional variability is evident, as some deltas deviate from the overall trend (for example, Indus and Cauvery deltas). When excluding outliers (the Indus and Cauvery deltas), subsidence and urban expansion exhibit a strong linear correlation across deltas (Extended Data Fig. 8c).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig13\"], \"section\": \"Urban expansion\", \"text\": \"We also explored the relationship among the anthropogenic drivers (Extended Data Fig. 8d\\u2013f), finding a low (R\\u2009=\\u20090.1\\u20130.3) correlation depending on the specific driver.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig3\"], \"section\": \"RF analysis for identifying anthropogenic drivers of subsidence and elevation loss\", \"text\": \"Given the nonlinear and interacting relationships among GWS, sediment flux alteration, urban expansion and residual land subsidence (after GIA correction) discussed above, a machine learning framework was implemented to model these complexities. First, we attempted a multilinear regression model, incorporating interaction terms between variables, formulated aswhere VLM is the predicted VLM, x0 is the intercept, X are the predictor variables (GWS, sediment flux alteration and urban expansion), x are the regression coefficients for each predictor variable, x represents the interaction effects between predictor variables and \\u03f5 is the residual error term. However, this multilinear regression model yielded poor performance (correlation R\\u2009=\\u20090.38; R2\\u2009=\\u20090.15; RMSE\\u2009=\\u20094.7\\u2009mm\\u2009yr1) (Fig. 3a), demonstrating the inefficiency of linear models to capture these complex dependencies and the need for a machine learning model.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Equ12\"], \"section\": \"RF analysis for identifying anthropogenic drivers of subsidence and elevation loss\", \"text\": \"The RF regressor optimizes each decision tree using the mean square error (MSE) defined as a cost function to identify node splits and model performance during model training and testing:where VLM is the observed VLM rate for individual delta i,  is the predicted VLM rate and N is the total number of observations. To assess uncertainty, we used Monte Carlo simulations to create multiple holdout fractions (0.1\\u20130.5) across 100 iterations, randomly subsampling the 40 deltas for training and validation in each iteration. This random partitioning ensures that each delta is used in both training and validation phases across iterations, enhancing the robustness against overfitting and sampling bias. The final RF model predictions were obtained by averaging prediction estimates across all iterations. The final model performance was evaluated using the coefficient of determination (R2), RMSE and mean absolute error (MAE):where  is the mean observed VLM rate, and the other variables are defined in equation (12). The feature importance If for input feature {X\\u2009=\\u2009X1,\\u2009X2,\\u2009X3} was computed using the following equation, based on the cumulative reduction in node, j impurity among all the trees:where N denotes the total number of trees and \\u0394I denotes the change in impurity.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM2\", \"MOESM2\", \"Fig3\"], \"section\": \"RF analysis for identifying anthropogenic drivers of subsidence and elevation loss\", \"text\": \"Although RF effectively captures nonlinear relationships, its ensemble structure limits delta-specific interpretability. To resolve local insights into delta-specific subsidence drivers, we applied LIME, a technique within the field of explainable artificial intelligence (XAI). LIME approximates black-box models such as RF by fitting interpretable models to perturbed samples of the input data, allowing for local feature importance estimation. For each delta X, LIME approximates the RF prediction locally by using a linear surrogate model trained on perturbed instances around X. The explanation function is obtained by solving the following minimization problem:where \\u03be(X) is the local interpretable model for each delta X, g is the interpretable model, f is the RF model,  is a proximity kernel,  is the loss function measuring the differences between f and g, and \\u03a9(g) penalizes complexity. This process was repeated for each delta, and deltas with low LIME model fidelity (R2\\u2009<\\u20090.5) were excluded to ensure reliable interpretation (Supplementary Table 2). The final dataset for interpretation consisted of 30 deltas, in which LIME produced more consistent feature importance estimates. The feature importance scores from LIME are normalized to obtain normalized LIME (nLIME) scores:where \\u03c9f is the LIME-derived coefficient for feature f and F is set for all features. The nLIME scores provide an instance-specific (local) explanation rather than a global one to evaluate the relative contributions of GWS, sediment flux alteration and urban expansion in each delta. The nLIME values for each delta are summarized in Supplementary Table 2 and were analysed in a ternary diagram to visualize the heterogeneity in delta-specific subsidence and elevation-loss drivers (Fig. 3b).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Equ6\", \"Equ8\", \"MOESM1\", \"MOESM2\", \"Fig5\"], \"section\": \"Historical, current and projected SLR rates\", \"text\": \"Historical relative sea-level changes were obtained from the Revised Local Reference database of the Permanent Service for Mean Sea Level (https://psmsl.org), which provides monthly relative sea-level records from globally distributed tide gauge stations. These tide gauge records have undergone quality control procedures, including corrections for datum inconsistencies, jumps and spurious data points, and validation through comparisons with neighbouring tide gauge stations. For this study, we selected 20 tide gauge stations across 15 deltas (the Mississippi, Rio Grande, Fraser, Amazon, Chao Phraya, Mekong, Red River, Nile, Ganges\\u2013Brahmaputra, Vistula, Rhine\\u2013Meuse, Chikuma-gawa, Yangtze, Pearl and Rioni deltas), considering only stations within 100\\u2009m of the delta boundary and at least 5\\u2009years (twentieth century) of valid record. The RSLR rates for each delta were estimated by applying the stochastic-seasonal model (equations (6\\u20138)) over the full observational record for each tide gauge. For deltas with multiple stations (for example, the Mississippi, Ganges\\u2013Brahmaputra and Rhine\\u2013Meuse deltas), individual station rates were averaged to provide a delta-wide estimate of twentieth-century RSLR. Note that the representativeness of the derived RSLR may vary for each delta following individual tide gauge characteristics (for example, is the station founded on bedrock or \\u2018floating\\u2019 in unconsolidated sediments, is the station GNSS corrected). Supplementary Fig. 7 shows the time series of relative sea level over the twentieth century for six representative deltas. Supplementary Table 1 provides a complete summary of the RSLR rates for the 15 deltas. The median twentieth-century RSLR trend across all deltas is 2.9\\u2009mm\\u2009yr\\u20131, with measured rates ranging from \\u20130.5\\u2009mm\\u2009yr\\u20131 in the Amazon delta (indicating declining twentieth-century sea level) to a maximum rate of 1.5\\u2009cm\\u2009yr\\u20131 in the Chao Phraya Delta (Fig. 5b).\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"MOESM1\", \"Equ6\", \"Equ8\", \"MOESM2\", \"Fig1\", \"MOESM2\", \"MOESM2\", \"MOESM1\"], \"section\": \"Historical, current and projected SLR rates\", \"text\": \"To estimate present-day (early twenty-first century) absolute (geocentric) SLR rates, we used the multi-mission satellite altimetry data from 2001 to present, obtained from Copernicus Marine Environment Monitoring Service (CMEMS). This dataset provides 1/8\\u00b0 (about 12.5\\u2009km) gridded monthly sea level anomalies (SLA) referenced to a 20-year mean baseline (1993\\u20132012). SLA estimates are derived from optimal interpolation, merging the level 3 along-track measurement from multiple contemporaneous altimeter missions (Jason-3, Sentinel-3A, HY-2A, Saral/AltiKa, Cryosat-2, Jason-2, Jason-1, TOPEX/Poseidon, ENVISAT, GFO and ERS1/2) (https://marine.copernicus.eu/). Several necessary corrections have been applied to the raw altimetry data, including instrumental biases and drifts, geophysical, tidal and atmospheric corrections, to ensure accurate SLA estimates. Monthly mean sea-level anomalies were obtained for each delta by spatially averaging the altimetry grid points within a 100-m radius, culling outliers beyond the 95th percentile. Supplementary Fig. 8 shows the monthly SLA time series in six deltas. We estimated the twenty-first-century trends in sea-level anomalies, using equations (6\\u20138). The altimetry-derived geocentric SLR rates for the twenty-first century show exacerbating regional SLR rates over global sea-level estimates (about 4\\u2009mm\\u2009yr\\u20131) for 45% of the deltas (18 out of 40) (Supplementary Table 1). Regional sea-level rates vary from 0.2\\u2009mm\\u2009yr\\u20131 in the Parana delta to 7.3\\u2009mm\\u2009yr\\u20131 over the Mississippi delta (Fig. 1 and Supplementary Table 1). However, a negative geocentric sea-level rate of \\u22121.9\\u2009mm\\u2009yr\\u20131 was observed in the Rioni Delta (Black Sea) (Supplementary Table 1). This long-term sea-level decline in the twenty-first century persists in the background of short-term fluctuations (Supplementary Fig. 8d); a characteristic feature of Black Sea sea-level dynamics. This twenty-first-century decline in geocentric sea level for the Rioni Delta represents more than a 100% reduction compared with historical (twentieth-century) rates, even when accounting for average VLM across the delta. To investigate this anomaly, we estimated VLM at the Poti tide gauge (Rioni Delta) by differencing twenty-first-century RSLR rates obtained from the Poti tide gauge station from geocentric SLR. The resulting VLM rate of \\u22126.7\\u2009mm\\u2009yr\\u20131 matches the average InSAR-derived VLM rate (\\u22125.9\\u2009\\u00b1\\u20090.7\\u2009mm\\u2009yr\\u20131) within 100\\u2009m of the tide gauge. This rapid subsidence rate at the coast of Poti represents localized conditions and highlights the need for caution when extrapolating point-based tide gauge measurements to infer delta-wide or city-wide subsidence and exposure. Note that satellite altimetry data, although highly valuable for global sea-level monitoring, were primarily optimized for open ocean conditions. Coastal environments naturally exhibit additional complexity due to processes such as shelf circulation, freshwater discharge and tidal amplification, which contribute to the inherent variability in nearshore sea-level measurements compared with offshore altimetric observations.\"}, {\"pmc\": \"PMC12823437\", \"pmid\": \"41535473\", \"reference_ids\": [\"Fig4\"], \"section\": \"Historical, current and projected SLR rates\", \"text\": \"We use projected sea-level rates from the Intergovernmental Panel on Climate Change Sixth Assessment Report (AR6) to assess future SLR rates across all deltas. The sea-level rate projections integrate process-based models that account for the key contributors to climate-induced sea-level change, such as thermal expansion, ocean dynamics, and glacier and ice sheet mass loss, and consider uncertainties in global temperature change and their influence on sea-level drivers. We focus on the no-VLM 50th percentile (median) projected rates for 2050 (mid-twenty-first century) and 2100 (end of the twenty-first century) under shared socioeconomic pathway 2-4.5 (SSP2-4.5) and SSP5-8.5 scenario. SSP5-8.5 represents a high reference scenario associated with the highest emission levels (global atmospheric CO2 concentrations exceeding 800\\u20131,100\\u2009ppm by 2100) and associated warming of 3.3\\u20135.7\\u2009\\u00b0C (refs.\\u2009). These projections provide an upper-bound reference scenario, capturing the potential worst-case outcome for future SLR. Figure 4c shows the comparison of projected SLR rates with observed land subsidence rates.\"}]"

Metadata

"{\"issue-copyright-statement\": \"\\u00a9 Springer Nature Limited 2026\"}"