Global climate change and local land subsidence exacerbate inundation risk to the San Francisco Bay Area

Coastal flooding is likely to be the biggest socioeconomic impact of sea level rise (SLR) in the 21st century (1–6). Several U.S. states (including Texas, Louisiana, and Florida, as well as U.S. territories such as Puerto Rico and the Virgin Islands) have just experienced devastation caused by the unprecedented flooding after Hurricane Harvey, Hurricane Irma, and Hurricane Maria. Storm intensity, associated rainfall, and storm surges affecting the coastal area were likely amplified by the elevated ocean temperature caused by ongoing global climate change (7). Moreover, increasing rates of mean global SLR, coinciding with the global rise in Earth’s temperature, increase the susceptibility of coastal areas to flooding. The increase in the rate of global SLR, from ~1.7 mm/year in the late 20th century to ~3.1 mm/year in the early 21st century, is mostly attributed to accelerated ice mass loss at major glaciers and polar ice caps and to thermal expansion of ocean water (8, 9). Although the contribution of these sources to SLR is fairly well understood, the rate of relative SLR (RSLR) can vary significantly due to a number of other processes such as isostatic adjustments, ocean currents, earthquakes, and volcanic episodes, as well as local land subsidence (LLS) associated with sediment and aquifer-system compaction (10–13).

The LLS may substantially increase the flooding risk associated with storms and SLR (10, 14, 15) by contributing to the rate of RSLR, depending on the properties of sediments, grain size, fluid content, loading history, and rate of pumping and production of subsurface fluids (16–18). Despite its importance, however, accurate and comprehensive measurements of LLS are generally unavailable, relying mostly on sparsely distributed geodetic and tide gauge measurements. Space geodetic measurements using the global navigation satellite system (GNSS) and interferometric synthetic aperture radar (InSAR) can provide much improved spatial coverage and resolution to characterize LLS (10, 14, 15, 19). Here, we report on accurate observations of LLS in the San Francisco Bay Area (SFBA) in a continental reference frame from a combination of GNSS and InSAR observations at an unprecedented resolution and accuracy. This allows us to update projections of RSLR in the SFBA and accordingly revise inundation risk maps.

The spatial pattern of SLR in the northeast Pacific is generally affected by climate patterns (such as El Niño and La Niña), changes in gravitational pull due to large glaciers and ice melting, isostatic rebound caused by ice sheet loss, and tectonic processes, including interseismic land uplift and subsidence (8, 9). The regional vertical land motion due to current tectonic and nontectonic deformation along the northeast Pacific shorelines varies from uplift of 1.5 to 3 mm/year north of Cape Mendocino to subsidence at a rate of ~1 mm/year to the south (8, 9, 20). Projections of SLR for the northeast Pacific suggested by the Intergovernmental Panel on Climate Change (IPCC) in 2007 (8, 9) predict a rise of 4 to 30 cm by 2030, 12 to 61 cm by 2050, and 42 to 167 cm by 2100, relative to 2000. Considering different emission scenarios, so-called representative concentration pathways (RCPs), and the upgraded rate of Antarctica ice sheet loss in the 21st century (21), the likely ranges of projected SLR can be much higher than that computed by IPCC (table S1) (2). The projections provided by the IPCC have already been used to identify areas of potential inundation in coastal regions. These maps are the basis for planning flood resilience strategies for large coastal urban areas such as the SFBA, the greater Los Angeles Basin, and New York (1). Heberger et al. (22) consider a 100-year 1.4-m SLR along the entire coast of California and estimate nearly $100 billion in property losses, 8750 km² of inundated wetlands, and 480,000 displaced people. These estimates are highly sensitive to the assumed SLR projection (23) and contributions from spatially variable LLS. Here, we use the SFBA as a case example to develop improved inundation scenarios, integrating high-resolution digital topography, detailed and accurate estimates of coastal LLS, and probabilistic projections of SLR under a wide range of emission and global warming scenarios.

Tide gauge measurements in the SFBA dating back to 1850 (Fig. 1) show a positive trend, indicating an RSLR at an average rate of ~2mm/year. To characterize the LLS and its contribution to RSLR in the SFBA, we first combine large sets of SAR data with repeated measurements of horizontal GNSS positions to obtain accurate, high-resolution, and multitemporal measurements of the three-dimensional (3D) surface deformation field across the SFBA relative to a local reference point (GNSS station LUTZ). The SAR data sets, spanning from 13 July 2007 to 17 October 2010, include 32, 24, and 19 images in descending and ascending orbits of the Environmental Satellite (Envisat) C-band satellite and ascending orbit of the Advanced Land Observation Satellite (ALOS) L-band satellite, respectively (Fig. 1B and table S2). This period is chosen because data sets from three different viewing geometries are available. Using these data sets, we generate more than 650 differential interferograms (see Methods). Least squares joint inversion of this large interferometric data set with observations from a dense GNSS network allows for resolving the 3D deformation field at 100-m ground resolution and 1 to 10 mm/year precision level.

Figure 2 (A to C) shows the obtained E, N, and U velocities relative to GNSS station LUTZ at the location of 379,286 elite pixels. The maps of E and N velocities are dominated by the long wavelength signal due to shear-strain accumulation across the San Andreas Fault system and short-wavelength deformation associated with shallow creep along the Hayward Fault. The comparison with horizontal velocities of 29 GNSS checkpoints relative to LUTZ yields an SD of 2.19 and 2.34 mm/year for the difference between estimated E and N velocities, respectively. For the comparison, we average the measurements of pixels within 250-m distance of a given GNSS station. The time series of E, N, and U displacements with respect to LUTZ are further validated using observations of four continuous GNSS stations (locations shown in Fig. 1A), for which we find good agreement between the GNSS- and InSAR-derived time series (fig. S3).

To put the vertical deformation estimates into a continental framework, we use data from four continuous GNSS stations (inverted triangles in Fig. 1A) with high-quality time series in a stable North American reference frame (NA12) (24). We apply a 1D conformal transformation (shift and scale factor) to transfer the results from the local combined InSAR and GNSS measurements into NA12 (see Methods for details). We find an SD of 0.52 mm/year for the difference between our estimated U velocities and those observed at the continuous GNSS stations provided by Blewitt et al. (24).

The map of vertical velocities shows a complex pattern of uplift and subsidence across the SFBA. Most of the Pacific shorelines and areas adjacent to the San Francisco Bay are subject to subsidence at less than ~2 mm/year. Portions of Treasure Island, San Francisco, San Francisco International Airport, and Foster City are subsiding as fast as 10 mm/year (Figs. 2D). Santa Clara Valley to the south of the Bay is characterized by uplift at 1 to 2 mm/year (Fig. 2C), likely due to rising groundwater levels during this period (25). We find that most of the subsiding pixels are located on Quaternary substrate such as Holocene Bay mud deposits or man-made landfills subject to long-term compaction (17, 26).

We combine the LLS map in the continental reference frame (Fig. 2D) with high-resolution topographic light detection and ranging (LiDAR) data to update projections of RSLR rate and develop revised inundation risk maps. To this end, we apply an approach similar to that of Knowles (23) but consider various probabilistic projections for the likely range (that is, 67% probability) of SLR under different emission and global warming scenarios (table S1). We also consider a recent scenario, referred to as H++, under which higher rates of Antarctica and Greenland ice loss will develop throughout the second half of the 21st century (21). The baseline for these SLR projections is the average relative sea level during the period 1991–2009 (2). We compiled more than 850 1 km × 1 km tiles of a DEM generated from airborne LiDAR measurements with an average spatial resolution of 1 to 2 m and estimated error of 0.1 to 0.2 m, covering the entire SFBA coastal areas (fig. S4). The DEM heights are given in the North American Vertical Datum of 1988. The DEMs and estimated LLS rates are interpolated on a regular grid of 2 m. We examine two sets of inundation forecast scenarios, one of which considers only the effect of SLR and the other includes the effect of LLS in addition to SLR. To obtain an estimate of LLS by 2100, we extrapolate the observed 2007–2010 LLS rates using a linear function. The formal SD associated with this extrapolation is less than 5 cm in 2100 (fig. S5), well below SLR projection and DEM errors; thus, we do not account for them when we evaluate the uncertainties of inundation hazard. However, the assumption of constant compaction rates over the 21st century requires some discussion. For the Holocene Bay mud deposits of variable thickness, which formed since the post ice-age SLR, these first-order compaction calculations (fig. S6) (12, 17) show that rates can be expected to be effectively constant over the next 100 years, independent of the material and geotechnical parameters chosen.

The potential inundation maps due to SLR alone and a combination of LLS and SLR for the upper bound of the likely range of RCP 2.6 and RCP 8.5 scenarios are shown in Fig. 3 (A and B). The former scenario corresponds to the goals of the 2015 Paris agreement, whereas the latter considers the scenario that no effort is made to mitigate and remove emissions. Corresponding results for the other SLR scenarios, provided in table S1, are shown in figs. S7 to S15. Figure 3C also summarizes the potential inundated areas for the ranges of considered SLR projections with and without accounting for LLS. Figure 4 shows examples of areas subject to significant inundation by 2100. The extents of areas that will be inundated because of either LLS or SLR alone, as well as their combined effect, are shown, amplifying the extent of inundated area. We find that SLR alone poses a significant inundation risk to coastal urban areas and infrastructures, as well as ecologically valuable marshes and wetlands (27). Considering the lower and upper bounds of likely ranges of various RCP-projected SLR scenarios and LLS, we estimate that, in 2100, an area of 98 to 218 km² will be affected by RSLR. The corresponding values for the case of SLR only are 51 to 168 km². Considering the H++ scenario, an area of 413 to 429 km² will be vulnerable to inundation (fig. S15). Even if SLR was completely halted, LLS alone would put 45 km² at risk. Thus, a much larger area will be affected by inundation once the effect of LLS is taken into account, especially for the more modest SLR scenarios. Considering the full range of various RCP scenarios and the H++ one, we estimate that the maps of potential inundation based solely on SLR projections for the SFBA underestimate the affected land area by 3.7 to 90.9%. These estimates are conservative, in that high tides, storm events, and high-precipitation periods can substantially raise the water level above that obtained from projections of SLR.

Major consequences of exacerbated inundation risk for coastal areas include saltwater contamination of surface and underground waters, accelerated coastal erosion, wetland losses, and increased flooding, which may lead to unprecedented socioeconomic impacts (1, 3, 5, 6). Heberger et al. (28) estimate that, in California by 2100, more than 480,000 people and $100 billion worth of property will be exposed to flood risk caused by SLR. Worldwide in 2005, more than 40 million people lived in coastal areas prone to 100-year flood risk; this number will grow more than threefold by 2070, and the value of property exposed to flooding will increase to ~9% of the projected global gross domestic product, with the United States, Japan, and Netherlands being the countries with most exposure (3). As we demonstrate, coastal subsidence can significantly increase this exposure. The framework presented here to account for coastal land subsidence is transferable to other coastal cities and can be used to inform policy decisions affecting coastal activities. Global climate change is affecting the future inundation risks both through accelerating ice sheet melting (increasing the rate of SLR) and through more intense droughts, leading to unprecedented groundwater overdraft (29) and associated localized coastal land subsidence (30).

In an era in which climate change and SLR pose unprecedented threats to the coastal environment, urbanization, and population, Earth observation data, such as those provided by the European Space Agency Copernicus Sentinel-1A and Sentinel-1B satellites (launched on 3 April 2014 and 25 April 2016, respectively) and by the NISAR (NASA-ISRO SAR) satellite (the United States and India joint mission due to launch by 2021), will continue to provide crucial data to inform policy decisions (31). These data can be used to improve flood resilience plans for coastal megacities around the world, most of which are located in river deltas and lowlands subject to subsidence and with significant populations exposed to flood risk (32).

Funding: The contribution of M.S. was funded by NSF grants EAR-1357079 and EAR-1735630 and NASA grants 80NSSC170567 and NNX17AD98G, and R.B. acknowledges support by NASA grant NNX12AQ32G and by NOAA National Estuarine Research System Science Collaborative. Radar data were provided by the Alaska Satellite Facilities. We thank UNAVCO for the operation and maintenance of the PBO GNSS stations and data archiving. These data and services are provided by the UNAVCO Facility with support from the NSF and NASA under NSF Cooperative agreement no. EAR-0735156. Additional GNSS data for this study come from the BARD network, doi: 10.7932/BARD, operated by the University of California, Berkeley Seismological Laboratory, which is archived at the Northern California Earthquake Data Center (NCEDC), doi: 10.7932/NCEDC. The LiDAR data were collected through a collaboration among the California Ocean Protection Council, the National Oceanic and Atmospheric Administration Coastal Services Center, and the U.S. Geological Survey. Author contributions: M.S. led the analysis. R.B. assisted with the interpretation of the results. All authors contributed to the theory and to the writing of the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.