Selection of global river deltas
We selected 40 deltas globally, prioritizing 35 deltaic systems with the greatest exposed area and population currently below sea level, supplemented by five less-exposed deltas of local and regional significance and previously identified risks9. To assess the 35 deltas with the greatest exposure among global river deltas, we used 955 delineated delta boundaries in ref. 6 and identified coastal delta elevation below sea level using the DeltaDTM dataset v.1.1 (ref. 11) resampled to 3 arcseconds (100 m) and referenced to mean sea level51. Global delta population was estimated by aggregating 100 m resolution WorldPop population count for each delta, which is calibrated to the 2020 national population estimates from the United Nations population data52.
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.
The final selection of 40 deltas spans five continents (Asia, Africa, Europe, North America and South America) and 29 countries, encompassing deltas with noted and emerging environmental, geophysical and social vulnerabilities9,24, historically sinking river deltas2 and densely populated coastal megacities3,4,53.
SAR dataset
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.
SAR interferometric analysis
We processed each SAR frame or single-path multiple-frame coverage to generate high-spatial resolution maps of surface deformation for the 40 deltas using a multitemporal wavelet-based InSAR (WabInSAR) algorithm54,55,56,57. First, we generated 59,000 high-quality interferograms from the coregistered SAR images using GAMMA software58,59, with an interferogram pair selection algorithm57 optimized through dyadic downsampling and Delaunay triangulation. To minimize phase errors and to maximize the pixel density associated with dynamic surface changes over deltas (for example, flooding, vegetation growth or soil saturation), we screened the initial set of interferograms based on their coherence stability to exclude interferograms with high coherence variability, while maintaining a 50% temporal baseline coverage. The final selection retained about 55,000 interferometric pairs (93%) for further analysis. Moreover, we implemented a statistical framework to discard noisy pixels with average coherence less than 0.7 for distributed scatterers and amplitude dispersion of greater than 0.35 for permanent scatterers57. Next, we used a minimum cost flow phase unwrapping algorithm optimized for sparse coherent pixels60,61 to estimate the absolute phase changes of the elite (less noisy) pixels in each interferogram. We corrected all unwrapped interferograms for the effects of residual orbital error62 and minimized the effects of topography-correlated components of atmospheric phase delay and spatially uncorrelated DEM error by applying a suite of wavelet-based filters54. Last, we estimated the time series, velocities and standard deviation for each geocoded elite pixel along the line of sight (LOS) of the satellite using a reweighted least-squares optimization55. The standard deviation of the LOS velocity corresponds to the uncertainty of the regression slope derived from the least-squares fit. For each delta, the reference point was selected as the pixel corresponding to a global navigation satellite systems (GNSS) station within the processed SAR frame when available. In areas without GNSS stations, a preliminary reference point was randomly selected from pixels with average temporal coherence >0.85. Following initial processing, the reference point was refined by visually identifying stable ground features (for example, bedrock outcrops and deep-foundation structures) and low displacement variability (standard deviation <1 mm yr–1), then reprocessing with this final reference point. For large deltas requiring overlapping SAR frame coverage, the LOS velocities were arranged in a mosaic form to ensure seamless spatial representation across the entire delta.
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 tracks63,64,65. 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 \(\{{\mathrm{LOS}}_{\mathrm{ASC}},{\mathrm{LOS}}_{\mathrm{DES}}\}\). Given \(\{{\mathrm{LOS}}_{\mathrm{ASC}},\,{\mathrm{LOS}}_{\mathrm{DES}}\}\) and their associated variances \(\{{\sigma }_{\mathrm{ASC}}^{2},{\sigma }_{\mathrm{DES}}^{2}\}\) 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 by
$$[\begin{array}{c}{{\rm{L}}{\rm{O}}{\rm{S}}}_{{\rm{A}}{\rm{S}}{\rm{C}}}\\ {{\rm{L}}{\rm{O}}{\rm{S}}}_{{\rm{D}}{\rm{E}}{\rm{S}}}\end{array}]=[\begin{array}{cc}{C}_{{\rm{A}}{\rm{S}}{\rm{C}}}^{E} & {C}_{{\rm{A}}{\rm{S}}{\rm{C}}}^{U}\\ {C}_{{\rm{D}}{\rm{E}}{\rm{S}}}^{E} & {C}_{{\rm{D}}{\rm{E}}{\rm{S}}}^{U}\end{array}]\,[\begin{array}{c}E\\ U\end{array}]$$ (1)
... continue reading