VLM data
High-resolution VLM data are based on the Virginia Tech Earth Observation and Innovation (EOI) Lab VLM product, with spatially continuous coverage for the Pacific, Atlantic and Gulf coasts of the USA27,28,51,52,53. The dataset provides VLM measurements at millimetre-level accuracy and a resolution of about 50 m within a 100-km strip along the coasts of the USA. For each coast, the VLM rates were determined by integrating SAR images from Sentinel-1 A/B and ALOS-1 satellites between 2007 and 2020 (see Supplementary Table 23 for satellite frames used for each coast) with observations of horizontal and vertical velocities at global navigation satellite system (GNSS) stations. To produce the spatially continuous surface-deformation map, InSAR line-of-sight (LOS) displacements were generated for the numerous SAR frames along the coasts.
We use GAMMA software to process SAR datasets54,55 and the wavelet-based InSAR (WabInSAR) algorithm to perform post-processing and multitemporal analysis56,57,58,59,60 (Supplementary Fig. 4a). To this end, thousands of high-quality interferograms were generated and several wavelet‐based analyses were applied to the interferograms to denoise the pixels and reduce the effects of spatially uncorrelated DEM error56,57 and topographically correlated atmospheric phase delay57. Next, the velocity along the LOS direction for each pixel is calculated as the slope of the best-fitting line to the associated time series using a reweighted least-squares estimation. Last, the numerous SAR frames are mosaiced following Ojha et al.61 and a stochastic model, which combines the LOS velocities with the GNSS datasets, was adopted to generate a high-resolution map of the VLM rate (Supplementary Fig. 4b).
To implement the stochastic model, we resampled the LOS velocities of Sentinel-1 tracks onto the ALOS track and interpolated the GNSS velocities on the pixels within the ALOS track using a Kriging interpolation technique with inverse distance weighting. Thus, we obtain several (5) observations per pixel for each coast, including LOS observations and GNSS velocities. Let {y0, y1,…, ym} and \(\{{\sigma }_{0}^{2},{\sigma }_{1}^{2},\ldots {\sigma }_{m}^{2}\}\) be the interpolated LOS velocities and variances, respectively, for a given pixel, in which subscripts 0, 1,…, m indicate the available satellite observations (Sentinel-1/ALOS-1) and orbits (ascending/descending) for a given US coast (Atlantic/Gulf/Pacific), as defined in Supplementary Table 23. The stochastic model to combine the LOS velocities with the velocities of GNSS datasets to generate a seamless, high-resolution and accurate map of east (E), north (N) and vertical (U) motions is given by equation (1):
$$\begin{array}{c}{y}_{0,1,\ldots m}={C}_{{\rm{e}}}^{0,1,\ldots m}E+{C}_{{\rm{n}}}^{0,1,\ldots m}N+{C}_{{\rm{u}}}^{0,1,\ldots m}U+{\varepsilon }^{0,1,\ldots m}\\ {E}_{{\rm{G}}}=E+{\varepsilon }^{{\rm{e}}}\\ \begin{array}{c}{N}_{{\rm{G}}}=N+{\varepsilon }^{{\rm{n}}}\\ {U}_{{\rm{G}}}=U+{\varepsilon }^{{\rm{u}}}\end{array}\end{array}$$
(1)
in which C represents the unit vectors projecting 3D displacements onto the LOS, which is a function of the satellite heading and incidence angles, ε is the observation errors equal to the standard deviations (σ), E, N and U are the unknowns and EG, NG and UG are the observed interpolated east, north and up GNSS velocities, respectively. The solution to equation (1) is given by equation (2):
$$X={\left({A}^{{\rm{T}}}PA\right)}^{-1}{A}^{{\rm{T}}}PL$$
(2)
in which X represents the unknowns, A is the Green’s function given by the unit vectors (C), L is the observation and P is the weight matrix, which is inversely proportional to the observant variance (σ2). The parameters variance-covariance matrix is \({Q}_{XX}=\frac{{r}^{{\rm{T}}}Pr}{{\rm{df}}}{({A}^{{\rm{T}}}PA)}^{-1}\), in which r = L − AX and df is the degrees of freedom. The standard deviations (precision of the results) for each pixel on the Atlantic/Gulf/Pacific coasts are shown in Supplementary Fig. 5a. The spatial distribution of the standard deviation shows that most values are below 3 mm per year for the US Atlantic and Gulf coasts. However, there are a few hotspots of high standard deviation around the Chesapeake Bay area (US Atlantic coast) and around the coast of Florida (US Gulf coast). We note higher estimated standard deviation values in the US Pacific coast, specifically in Northern California and the Orange County basin27 (Supplementary Fig. 5a). Generally, higher standard deviation values represent areas of lower precision. The observed lower precision in some pixels may be attributed to lower interferometric phase signal-to-noise ratio caused by surface vegetation, nonlinearity in the rates between the ALOS and Sentinel-1 observation periods owing to aquifer recharge and depletion, a limited number of GNSS stations used for the adjustment and a comparatively higher standard deviation of the GNSS station in the particular regions20,27,28. Furthermore, we validate the VLM rates using 756 GNSS stations (US Atlantic coast: 218; US Gulf coast: 157; US Pacific coast: 381) from the Nevada Geodetic Laboratory62 and Shirzaei et al.48. To perform the validation, we computed the average InSAR VLM rates within a 200-m radius around each GNSS station for comparison with the corresponding GNSS vertical rates (Extended Data Figs. 2–4). We obtained a standard deviation of 1.5 mm per year and a mean difference of less than 0.3 mm per year for the US Atlantic, Gulf and Pacific coasts (Supplementary Fig. 5b–d).
Spatial analysis of the complied VLM map (Fig. 1 and Extended Data Figs. 2–4) reveals extensive coastal areas with subsidence rates of more than 3 mm per year. Figure 1e highlights the spatially variable VLM for the 32 major coastal cities selected for this study: US Atlantic coast: Boston, MA; New York City, NY; Jersey City, NJ; Atlantic City, NJ; Virginia Beach, VA; Wilmington, NC; Myrtle Beach, SC; Charleston, SC; Savannah, GA; Jacksonville, FL; Miami, FL; US Gulf coast: Naples, FL; Mobile, AL; Biloxi, MS; New Orleans, LA; Slidell, LA; Lake Charles, LA; Port Arthur, TX; Galveston, TX; Texas City, TX; Freeport, TX; Corpus Christi, TX; US Pacific coast: Richmond, CA; Oakland, CA; San Francisco, CA; South San Francisco, CA; Foster City, CA; Santa Cruz, CA; Long Beach, CA; Huntington Beach, CA; Newport Beach, CA; San Diego, CA.
We find subsidence rates greater than 2 mm per year in 24 out of 32 major cities along the US Atlantic, Gulf and Pacific coasts, with notable subsidence rates (>5 mm per year) in cities such as Charleston (city number 8), Biloxi (city number 14), Galveston (city number 20) and Corpus Christi (city number 22) (Fig. 1e). On the US Pacific coast, we observe lower rates of land subsidence compared with the Atlantic and Gulf coasts, with some cities characterized by marked uplift (such as Richmond: city number 23; Long Beach: city number 29; Huntington Beach: city number 30; and Newport Beach: city number 31).
Subsidence along the coast is driven by natural and human processes and is a notable contributor to relative sea-level change2,19,48,49,63,64,65. Earlier studies suggested that complex processes drive observed subsidence along the US coasts27,66,67,68,69. These drivers include a combination of natural and anthropogenic processes, such as GIA, compaction of sediments, groundwater withdrawal, hydrocarbon extraction, surficial drainage/dewatering activities and regional tectonic activities. On a broad scale, disentangling the contribution of natural-driven and anthropogenic-driven processes is important for developing effective strategies to mitigate or adapt to the impacts of subsidence in low-lying coastal cities. On the one hand, in cities in which subsidence is a result of GIA and other natural processes, effective subsidence mitigation will probably involve an adaptive response, such as raised structures and infrastructure and flood-protection measures. On the other hand, for anthropogenic processes, proactive policy interventions and mitigation measures to reduce and control resulting subsidence, such as reducing groundwater and oil and gas extraction or changes in land use, may be helpful in sinking cities. As GIA is the main natural driver, we used the GIA ICE-6G-D model70 to estimate the GIA contributions at the SAR pixels and subtracted its effect from the observed VLM to assess the non-GIA contributions to the estimated VLM along US coasts (Supplementary Fig. 2). The relative reduction of subsidence by 46%, 4% and 20% for the Atlantic, Gulf and Pacific coasts, respectively, suggests that the effect of GIA on subsidence is dominant primarily along the US Atlantic coast and minimal for the Gulf and Pacific coasts (Supplementary Fig. 2c–e). Although the median rates of subsidence are reduced for all 32 major coastal cities, several areas with subsidence rates greater than 2 mm per year remain apparent in more than half of the selected cities, such as Boston, Atlantic City, Charleston, Biloxi, New Orleans, Texas City, San Francisco, Foster City and San Diego (Supplementary Fig. 3).
Coastal cities selection and elevation data
To select the 32 cities for analysis, we considered 41 major US coastal cities with VLM and LiDAR DEMs data. We conducted a preliminary analysis to determine the exposed area of each city, considering the IPCC localized (relative) SLR projections and mean high water (MHW) of the nearest tide gauge. Next, we screened the cities on the basis of the largest exposed area and selected ten cities from each US coastal region as follows: US Atlantic coast: Boston, MA; New York City, NY; Atlantic City, NJ; Virginia Beach, VA; Wilmington, NC; Myrtle Beach, SC; Charleston, SC; Savannah, GA; Jacksonville, FL; Miami, FL; US Gulf coast: Naples, FL; Mobile, AL; Biloxi, MS; New Orleans, LA; Slidell, LA; Lake Charles, LA; Port Arthur, TX; Texas City, TX; Freeport, TX; Corpus Christi, TX; US Pacific coast: Richmond, CA; Oakland, CA; San Francisco, CA; South San Francisco, CA; Foster City, CA; Santa Cruz, CA; Long Beach, CA; Huntington Beach, CA; Newport Beach, CA; San Diego, CA.
On the Pacific coast, we focused only on future inundation hazards for cities in California. The absence of coastal cities from Oregon and Washington is a result of the lack of high-resolution VLM data for the US northwest coast (Fig. 1) and the complexities of future inundation hazards in the region driven by earthquake and tsunami hazards. The aftermath of earthquake and tsunami hazards can cause substantial subsidence followed by inundation from tsunami waves. Evaluating such hazards requires a probabilistic analysis of future earthquake and tsunami hazards beyond this study.
Also, we added two cities, Jersey City, NJ and Texas City, TX that were located near other selected cities (New York City, NY and Galveston, TX, respectively) and are also important urban centres in their respective regions.
We use LiDAR DEM for the coastal-elevation data. The high-resolution LiDAR DEMs hosted by the National Oceanic and Atmospheric Administration (NOAA) Office for Coastal Management are available for the coastal USA71. In this study, we used a 3-m × 3-m grid resolution for the 32 cities, except Savannah (GA), Jacksonville (FL), Miami (FL) and all cities on the Pacific coast, which were obtained at a 5-m × 5-m grid resolution (Supplementary Table 23). All DEMs for each city use the North American Vertical Datum of 1988 (NAVD 88) vertical datum. Details on the implementation, vertical and horizontal accuracy, errors and temporal range are available with the data download71.
Population, properties and racial demographic data
We estimate the population and property datasets for each city using the Topologically Integrated Geographic Encoding and Referencing (TIGER) system demographic and economic data records available from the US Census Bureau (https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-data.2010.html). The dataset provides population and property estimates for each city in the USA, subdivided into census blocks based on the 2010 census data. We used the 2010 dataset because it is the most recent census data available from the US Census Bureau. The racial demographic dataset is based on the decennial Census data (https://data.census.gov/cedsci/advanced) corresponding to the 2010 census. For this study, we selected eight races: ‘White’, ‘Black or African American’, ‘American Indian and Alaska Native’, ‘Asian’, ‘Native Hawaiian and Other Pacific Islander’, ‘Hispanic or Latino’ and others (‘Some Other Race’ alone and ‘Two or More Races’), as defined by the decennial data.
Sea-level projections
We use the localized sea-level projections from the IPCC Sixth Assessment Report17,34. The projections consider the contributions to future sea levels from sterodynamic effects (ocean steric and ocean dynamic effects), ice sheets (Antarctic and Greenland ice sheets), land water storage, glacier and ice cap surface mass balance, thermal expansion and the IPCC estimates of total VLM based on tide-gauge observations—reflecting the sum of GIA and other VLM processes. To prevent double counting of VLM, we acquired the SLR projections without the effect of VLM for our analysis (geocentric SLR). The database provides projections of sea level at tide-gauge stations worldwide under five SSP scenarios (SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0 and SSP5-8.5). SSP1-1.9 limits warming to 1.5 °C above 1850–1900 levels by 2100, implying net-zero CO2 emissions around the middle of the century. SSP1-2.6 keeps warming below 2.0 °C relative to 1850–1900, with projected net-zero emissions in the second half of the century. The SSP2-4.5 scenario projects best-estimate warming of approximately 2.7 °C by the end of the twenty-first century relative to 1850–1900. SSP3-7.0 is a medium to high reference scenario resulting from no further climate policy with particularly high non-CO2 and aerosol emissions and a warming of 2.8–4.6 °C. SSP5-8.5 is a high reference scenario with the highest emission levels (above the current emissions trajectory) and warming of 3.3–5.7 °C. In this study, we apply the 17th (lower bound), 50th (median) and 83rd (upper bound) percentile projections under SSP2-4.5, which represents the current emissions trajectory. Supplementary Table 24 summarizes the tide-gauge stations used for the SLR projections in each city.
High-tide estimates
We used MHW at tide gauges to estimate the high-tide events for each coastal city. The tide-gauge measurements were obtained from NOAA tide and currents data72, using the NAVD 88 datum, consistent with the elevation datum and mean measurement for the present epoch. Tide-gauge stations used for each city were selected on the basis of the proximity to the city, which provides localized data crucial for accurate evaluation of the current exposure to high tide in the urban areas (Supplementary Table 24).
Inundation model
Using a bathtub model2,51,73,74,75 (see Supplementary Fig. 4c), we projected the inundation hazards for 32 cities on the US coasts. The input data for the inundation model are as follows:
-
1.
3-m or 5-m grid LiDAR DEM for each city.
-
2.
About 50-m resolution VLM data for each city.
-
3.
MHW levels at tide gauges adjacent to each city (Supplementary Table 24).
-
4.
IPCC geocentric SLR projections at the stations adjacent to each city (Supplementary Table 24).
To provide a comprehensive exposure assessment, we incorporate two temporal scales—the current (2020) and projected exposure (2050). First, the current exposure to high tide is assessed using MHW levels. Subsequently, projected exposure is evaluated by considering both VLM and geocentric SLR. Thus, the projected exposure represents further exposure, providing a baseline of current exposure against which future scenarios can be compared. To implement the inundation model, first, we resample the VLM rates on the LiDAR DEM. Next, we modify the elevation model to account for VLM projections, assuming a linear VLM rate from the base year of the DEM to the target years of 2020 and 2050 (refs. 34,48,51). Last, we evaluated the current (2020) scenario by subtracting the modified DEM height, which accounts for VLM projections up to the year 2020, from MHW levels (equation (3)). Subsequently, for the 2050 scenario, we apply SLR projections by subtracting the modified DEM height, updated for VLM projections to the year 2050, from both the geocentric SLR projection height and the MHW levels (equation (3)). Areas with a projected height below zero are inundated. This simplified static model is useful for local-scale simulations of inundated locations hydrologically connected to the coast74. However, it may overestimate or underestimate inundated areas on the coast owing to the reduced complexity of the model74. To reduce the errors associated with this approach, we implemented connected-component analysis to remove solitary grid cells from the inundation model, which represents topographically isolated low regions. Furthermore, we first present our inundation model as an undefended inundation map that does not account for the presence of levees or sea walls and we introduce and discuss the possible implications of flood-control structures on the impacts of relative SLR. For the defended scenario, we account for existing levees and sea walls by modifying the DEM height at the location of flood-control structures above the threshold for potential flooding.
To account for all error sources in the input data, we consider the uncertainties in the DEM, VLM and SLR datasets. Specifically, we propagate the 17th and 83rd percentiles for geocentric SLR projections, ±1 standard deviation for VLM and errors inherent in the DEM (equation (3)). These measures provide an error bound for the inundation analysis, ensuring robust estimation of the uncertainties associated with the projections.
$$\begin{array}{l}{{\rm{Inun}}}_{{\rm{med}}}={{\rm{DEM}}}_{{\rm{mod}}}-({{\rm{SLR}}}_{50}+{\rm{MHW}})\\ {{\rm{Inun}}}_{{\rm{low}},{\rm{up}}}={{\rm{Inun}}}_{{\rm{med}}}\pm (\sqrt{{{\rm{DEM}}}_{{\rm{err}}}^{2}+{((t-{t}_{0}){{\rm{VLM}}}_{{\rm{SD}}})}^{2}+{(({{\rm{SLR}}}_{83}-{{\rm{SLR}}}_{17})/2)}^{2}})\end{array}$$
(3)
in which Inunlow, Inunmed and Inunup represent the median, lower and upper bounds, respectively, for the inundation models. DEMmod is the modified DEM height, updated using the VLM projections. DEMerr is the vertical accuracy of the DEM. t represents the projection target years of 2020 or 2050. t0 represents the base year of the DEM. VLMSD is one standard deviation from the VLM data. MHW represents mean high water. SLR17, SLR50 and SLR83 represent the 17th, 50th and 83rd percentiles, respectively, from the geocentric SLR projections. Note that, for evaluating the current/baseline scenario (2020), SLR17, SLR50 and SLR83 are zero.
Socioeconomic exposure analysis
We used the TIGER demographic and economic data to assess the population and property exposure, which estimates the total population and properties subdivided into census blocks. We consider a census block inundated if greater than 20% of its area is inundated and assign the population and property for that block as exposed population or properties. To select the 20% threshold for the exposure of each census block, we conducted an empirical analysis across six representative cities. The distribution of exposed areas within these census blocks followed an extreme-value distribution. Statistical metrics revealed a median value of exposed area ranging from 18% to 23% for the six cities (Supplementary Fig. 6). Furthermore, the distribution exhibited a sharp decline beyond 10% (Supplementary Fig. 6). Therefore, a 20% criterion was established as a suitable threshold for quantifying the exposed population and properties. To quantify the home-value exposure, we used the ZIP Code ZHVI as a metric for housing cost. The estimated home-value exposure was calculated by multiplying the number of exposed properties within each city by the corresponding ZHVI (https://www.zillow.com/research/data/). We adjusted the ZHVI for the recent economic inflation using the mid-2021 housing price. The Zillow home-value data for each ZIP code used in this study is reported in Supplementary Table 25.
To investigate the disparate sociodemographic and socioeconomic impacts of relative SLR on vulnerable groups, we focused on analysing both racial and economic disparities in the exposed communities. To examine the racial disparities in the exposed communities, we considered eight races as defined by the US census decennial data: ‘White’, ‘Black or African American’, ‘American Indian and Alaska Native’, ‘Asian’, ‘Native Hawaiian and Other Pacific Islander’, ‘Hispanic or Latino’, ‘Some Other Race’ alone and ‘Two or More Races’. Racial minoritized groups are defined as individuals identifying as Black or African American; American Indian or Alaska Native; Asian; Native Hawaiian or Other Pacific Islander; Hispanic or Latino; and two or more groups38. The analysis shows an overrepresentation of the white population on the Atlantic and Pacific coasts, whereas minoritized populations are overrepresented on the Gulf coast. On the Atlantic coast, the white population makes up 55.1–71.4% of the exposed population by 2050, which is higher than their 38.3% share in the total population (Extended Data Fig. 7 and Supplementary Table 13). Minoritized groups make up 43.0% of the total population on the Gulf coast while accounting for 64.2–71.5% of the exposed population by 2050 (Extended Data Fig. 7 and Supplementary Table 14). On the Pacific coast, white residents comprise 57.6–70.9% of the exposed population by 2050, despite making up only 41.8% of the total population (Extended Data Fig. 7 and Supplementary Table 15). In cities such as Jersey City, New Orleans, Port Arthur and Oakland, minority groups are disproportionately represented among the exposed population (Supplementary Tables 13–15).
To assess the impacts of relative SLR on economic inequality, we used the property value as a proxy for economic status. Using a Kolmogorov–Smirnov statistical method, we compare the median home values in regions exposed to relative SLR by 2050 against those in each city, using an alpha value of 0.05 for statistical significance. Supplementary Tables 16–18 summarize the statistical test for the 32 cities. Note that 9–22 cities (considering lower to upper bound relative SLR projections) were excluded from the analysis because of limitations imposed by the central limit theorem. Across the 14 cities examined (considering median relative SLR projection), we find statistically significant economic disparities in 12 cities (Extended Data Fig. 8). In eight of these cities, we find that the median home value for the exposed population is higher than the total home value in the cities (that is, exposed properties are overvalued). However, in Atlantic City, New Orleans, Port Arthur and Foster City, we find that the median exposed-home values are lower than the overall median home value within the cities, highlighting their disproportionate economic exposure.
Subsidence hazard exposure analysis for levees
The polygon features for the levees across the US coasts were obtained from the United States Army Corps of Engineers (USACE)76. To determine the exposure to subsidence for the levees, we extracted the VLM rate for each point along the polygon feature. The subsidence exposure for levees in five cities (Miami, FL; New Orleans, LA; Port Arthur, TX; Freeport, TX; and Foster City, CA) are shown in Extended Data Fig. 9.
