Research ArticleGEOPHYSICS

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

See allHide authors and affiliations

Vol. 4, no. 3, eaap9234

## Abstract

The current global projections of future sea level rise are the basis for developing inundation hazard maps. However, contributions from spatially variable coastal subsidence have generally not been considered in these projections. We use synthetic aperture radar interferometric measurements and global navigation satellite system data to show subsidence rates of less than 2 mm/year along most of the coastal areas along San Francisco Bay. However, rates exceed 10 mm/year in some areas underlain by compacting artificial landfill and Holocene mud deposits. The maps estimating 100-year inundation hazards solely based on the projection of sea level rise from various emission scenarios underestimate the area at risk of flooding by 3.7 to 90.9%, compared with revised maps that account for the contribution of local land subsidence. Given ongoing land subsidence, we project that an area of 125 to 429 km2 will be vulnerable to inundation, as opposed to 51 to 413 km2 considering sea level rise alone.

## INTRODUCTION

Coastal flooding is likely to be the biggest socioeconomic impact of sea level rise (SLR) in the 21st century (16). 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 (1013).

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 (1618). 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 km2 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. ## RESULTS 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 km2 will be affected by RSLR. The corresponding values for the case of SLR only are 51 to 168 km2. Considering the H++ scenario, an area of 413 to 429 km2 will be vulnerable to inundation (fig. S15). Even if SLR was completely halted, LLS alone would put 45 km2 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. ## DISCUSSION 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).

## METHODS

### Multitemporal multitrack InSAR

The SAR data sets spanning from 13 July 2007 to 17 October 2010 include 32, 24, and 19 images in descending (incidence angle = 23°, heading angle = 193°) and ascending (incidence angle = 23°, heading angle = 350°) orbits of the Envisat C-band satellite and ascending (incidence angle = 34.5°, heading angle = 350°) orbit of ALOS L-band satellite, respectively (Fig. 1 and table S2). This period was chosen because data sets from three different geometries were available. Using these data sets, we generated 428, 151, and 77 interferograms, respectively. For the Envisat interferometric data set, the maximum temporal and perpendicular baselines are set to be 3 years and 300 m. The corresponding numbers are 3 years and 1500 m for the ALOS data set.

Each data set was processed separately to generate three accurate and high-resolution time series of surface deformation. To this end, we applied the wavelet-based InSAR (WabInSAR) algorithm (33, 34). The geometrical phase was calculated and removed using a 30-m DEM provided by the SRTM (35) (www2.jpl.nasa.gov/srtm/) and the satellite precise ephemeris data. The time series of complex interferometric phase noise was then calculated in the wavelet domain and analyzed in a statistical framework to identify elite (that is, less noisy) pixels. The absolute estimate of the phase change for elite pixels was obtained via an iterative 2D sparse phase unwrapping algorithm. Each unwrapped interferogram was corrected for the effect of orbital error (36) and the topography-correlated component of atmospheric delay (37). Using a robust regression, the phase changes from the large set of interferograms were inverted. This approach reduced the effects of outliers due to improper phase unwrapping. A high-pass filter, using continuous wavelet transforms, was applied to reduce the temporal component of the atmospheric delay. Lastly, we geocoded all data sets to obtain precise locations of elite pixels in a geographic reference frame. Figure S1 shows the line-of-sight (LOS) velocity associated with each data set.

Assume that {d0,,dNi}i=1,2,3 and {σ20,,σ2Ni}i=1,2,3 are the displacement time series and the associated variances for the data sets, where N is the number of acquisitions in each data set and i = 1,2,3 refers to Envisat descending, Envisat ascending, and ALOS ascending data sets, respectively. In the space domain and using a nearest neighbor algorithm, we oversampled data sets 2 and 3 over data set 1. In the time domain, the displacement and noise time series of one track were interpolated on the other. Assuming that, at time ta, da and σa2 are the displacement and associated SD, respectively, and that the respective values at time tb are db and σb2, then, at time step tc, the interpolated displacement and variance are given asdc=tctatbta(dbda)+daσc2=(tctatbtaσb)2+(tctatbtaσa)2+σa2(1)

To set up the interpolation in time domain, we first generated a vector with length M including acquisition dates of all three SAR data sets. This vector is the base for temporal oversampling, namely, all three data sets were oversampled on this new vector.

Assume {y0, y1, …, yM}i = 1,2,3 and {σ20, …, σ2M}i = 1,2,3as the interpolated time series of displacement and variance for a given pixel. Below, we set up a Kalman filter–based framework to combine these three InSAR time series with horizontal velocities of GNSS data sets (E-W and N-S components) to generate a seamless, high-resolution, and accurate time series of east (E), west (W), and U (vertical) motions. The GNSS velocities were obtained from Bürgmann et al. (26) and represented a 200-station subset of the regional BAVU velocity field (38) that relied on data from both continuous and campaign GNSS measurements. We split this data set into two parts: One group of stations was used to combine with the InSAR data sets (so-called tie points), and the other part was used to validate the results (so-called checkpoints) (Fig. 1). The E and N displacements of the GNSS tie points were interpolated on the location of the elite pixels from data set 1, using a Kriging interpolation approach and inverse distance weighting (39).

The measurement model at any time step tk relates the LOS observations to the corresponding E, N, and U displacementsytk1=Ce1Etk+Cn1Ntk+Cu1Utk+e1ytk2=Ce2Etk+Cn2Ntk+Cu2Utk+e2ytk3=Ce3Etk+Cn3Ntk+Cu3Utk+e3(2)where C represents the unit vectors projecting the 3D displacements onto the LOS and is a function of heading and incidence angles and e is the measurement noise equal to σ, which is assumed to be normally distributed.

The dynamic model is used to integrate GNSS and InSAR data and applies a statistical smoothing in time (40) and is given asEtk+1=(tk+1tk)×VeG+Etk+εeNtk+1=(tk+1tk)×VnG+Ntk+εn(3)where V represents the interpolated E or N velocities of the GNSS tie points and ε is the system dynamic noise estimated via Eq. 2. Grewal and Andrews (40) have provided a recursive solution to Eqs. 2 and 3. Because this approach is an extension to an existing InSAR time series algorithm (that is, WabInSAR), hereafter, we shall call it WabInSAR-3D.

To validate the WabInSAR-3D approach, we used a synthetic test. To this end, we simulated 3D deformation time series as a combination of a periodic function with a period of 1 year, a Heaviside function with step size of 5 and a linear trend with slopes of 3, 2, and 1 in E, N, and U directions, respectively (fig. S1). Using the geometrical parameters of the three SAR data sets, we projected the simulated 3D displacement field onto the LOS of Envisat descending and ascending and ALOS ascending orbit track, similar to our data set over the SFBA. The simulated LOS time series were inverted to recover the 3D displacement field. We first applied the method presented by Ozawa and Ueda (41), which jointly inverted data from different look angles and solved for E and U time series. The top row in fig. S1 shows the results from implementing the Ozawa and Ueda (41) approach together with the simulated E, N, and U time series. Using this approach, no N displacement was retrieved because they assumed that its contribution was negligible due to the near-polar orbiting satellites. For the E component, the periodicity, trend, and jump were recovered accurately, owing to multiple look and heading angles. However, the estimated U component is relatively uncertain, the mismatch between simulated and recovered signal increases near the jump, and the error propagates throughout the signal. The bottom row in fig. S1 presents results from applying the WabInSAR-3D. The estimated E component is accurate and matches the simulated signal slightly better than that of Ozawa and Ueda (41). For the N component, WabInSAR-3D is able to at least retrieve the trend and slight periodicity, but the jump is lost. The estimated U component here shows a significant improvement over the other approach and retrieves the periodicity, trend, and jump accurately. Now that we have validated WabInSAR-3D using a simulated data set, we applied it to the SAR and GNSS data sets over the SFBA.

The E, N, and U velocities relative to GNSS station LUTZ are shown in Fig. 2 (A to C), which are obtained from applying the WabInSAR-3D approach to 68 SAR images acquired by Envisat and ALOS satellites, together with horizontal velocities of 26 GNSS tie points, which are randomly selected. We obtained an SD of 2.19 and 2.34 mm/year, comparing the estimated E and N velocities, with horizontal velocities of 29 GNSS checkpoints not used in the combination. Using observations of four continuous GNSS stations of the BARD network relative to station LUTZ, the time series of E, N, and U displacement were also validated (fig. S3).

To transfer the vertical velocities from the central local reference frame relative to GNSS station LUTZ into a stable North America reference frame (NA12), we used vertical velocities of four permanent GNSS stations of the Plate Boundary Observatory (PBO) network that lie in the area of InSAR coverage (Fig. 1A), provided by Blewitt et al. (24), and applied a conformal translation. We found that a scale factor of 1 and uniform shift of 0.25 mm/year were sufficient to firmly tie the results from the combined InSAR and GNSS measurements to the continental reference frame. We obtained an SD of 0.52 mm/year for the transferred U velocity map in the NA12 reference frame to the four PBO stations.

### Compaction model

Using the 1D compaction model suggested by Kooi and de Vries (17), we ran simulations to investigate the impact of sediment thickness and accumulation rate on the present-day compaction rates, using the procedure detailed by Meckel et al. (12). For the Holocene Bay mud deposits of variable thickness (https://pubs.usgs.gov/mf/0976/plate-1.pdf), we found that rate variations are negligible over the next 100 years for different material and geotechnical parameters chosen. The most rapid subsidence (≫5 mm/year) was found in areas of anthropogenic landfill overlying thick Bay mud deposits, such as the northwest corner of Treasure Island where ~10-m-thick fill overlies >20-m-thick Bay mud (42). Here, we can draw on leveling data capturing the early phases of secondary compaction (42). Most of the compaction of the fill material was completed within less than 5 years after placement in 1936–1937. By 1960, the secondary compaction of the underlying Bay mud continued at a rate of ~20 mm/year (42). Subsidence rates in the northwest of Treasure Island relative to nearby bedrock on Yerba Buena Island had decayed to ~10 mm/year by 2000 (43) and still amounted to ~6 to 10 mm/year in our 2007–2010 data set. In these areas of recent fill emplacement, our linear extrapolation may overestimate the expected LLS in 2100 by 5 to 14%, assuming a sediment thickness of 10 to 30 m.

Correction (23 March 2018): The authors noted an incorrect funding source, Point Blue Conservation Science, in error. The correct funding details, NOAA National Estuarine Research System Science Collaborative, were added to the Acknowledgments.

## SUPPLEMENTARY MATERIALS

table S1. Projected SLR (in meters) for the Golden Gate tide gauge in San Francisco.

table S2. SAR data set used in this study, including Envisat data acquired in descending orbit track 70 and ascending orbit track 478, as well as ALOS SAR data obtained in ascending orbit, frame 740 and track 222.

fig. S1. Validation test using synthetic data sets.

fig. S2. LOS velocity associated with each data set.

fig. S3. Validation test of InSAR time series with GNSS measurements.

fig. S4. The distribution and LIDAR DEMs used in this study, as well as the associated sources.

fig. S5. One hundred–year projected uncertainty of the land subsidence from the InSAR- and GNSS-derived subsidence rates.

fig. S6. Present-day compaction rates.

fig. S7. Inundation map at 2030 given the lower bound of the likely range of SLR projection under RCP 4.5 scenario (table S1).

fig. S8. Inundation map at 2030 given the upper bound of the likely range of SLR projection under RCP 4.5 scenario (table S1).

fig. S9. Inundation map at 2050 given the lower bound of the likely range of SLR projection under RCP 4.5 scenario (table S1).

fig. S10. Inundation map at 2050 given the upper bound of the likely range of SLR projection under RCP 4.5 scenario (table S1).

fig. S11. Inundation map at 2100 given the lower bound of the likely range of SLR projection under RCP 2.6 scenario (table S1).

fig. S12. Inundation map at 2100 given the lower bound of the likely range of SLR projection under RCP 4.5 scenario (table S1).

fig. S13. Inundation map at 2100 given the upper bound of the likely range of SLR projection under RCP 4.5 scenario (table S1).

fig. S14. Inundation map at 2100 given the lower bound of the likely range of SLR projection under RCP 8.5 scenario (table S1).

fig. S15. Inundation map at 2100 given the SLR projection under H++ scenario (table S1).

Reference (44)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

Acknowledgments: 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.
View Abstract