Tracking California’s sinking coast from space: Implications for relative sea-level rise

The ~1350-km-long coast of California is monitored by several SAR satellites and networks of GNSS stations (Fig. 1) because of the interest in crustal deformation data for assessing seismic hazards at plate boundary faults (32). The SAR datasets used here include acquisitions from Sentinel-1A/B C-band satellites in descending and ascending orbit geometry during 2014–2018 and ALOS L-band satellite in ascending orbit geometry during 2007–2011 (table S1). We also used 441 GNSS stations whose three-dimensional (3D) velocities are known within the inertial terrestrial reference frame (ITRF) (Fig. 1B). Given the broad extent of California’s coast, the area is covered by multiple overlapping InSAR frames. We performed a multitemporal SAR interferometric processing of frames in each viewing geometry of each satellite separately (see Methods). Then, we used the data at the overlapped zones of adjacent frames to perform an affine transformation (including a shift, a scale factor, and two rotations) and obtained a seamless map of line-of-sight (LOS) displacement velocities for each viewing geometry along the coast of California (figs. S1 to S3). A least-squares joint inversion of these extensive LOS velocity fields with observations of horizontal displacement velocity at the location of GNSS stations allowed for resolving the 3D deformation field at ~100-m ground resolution and millimeter per year precision level. To translate estimated vertical deformation rates into a continental reference frame, we used vertical velocities of GNSS stations within the 2012 North American reference frame (NA12) (33) and applied an affine transformation to transfer the estimated VLM into NA12.

Figure 2A shows the spatial distribution of the VLM rates along the coast of California within the NA12 reference frame. Negative values indicate land subsidence. The map includes nearly 35 million pixels at ~100 m × 100 m resolution. The map is characterized by broad-scale (~300 km) zones of subsidence and uplift. Subsidence occurs with rates exceeding −5 mm/year in the San Francisco Bay Area, −8 mm/year near Santa Cruz, and −2 mm/year in Southern California. There are two notably large regions of uplift of 3 to 5 mm/year, one north of the San Francisco Bay Area (latitudes >39°N) and the other in Central California between latitudes ~37°N and 35°N. We identify two relatively small zones of substantial uplift in Central and Southern California, which are associated with the Santa Clara and Orange County basins, respectively. We also see subsidence in the Ventura basin, likely related to groundwater extraction, which agrees with results from Hammond et al. (34).

Figure 2B shows the spatial distribution of the SD (σ) for the estimated VLM. Most of the values are below 2.5 mm/year. However, in Northern California (latitude >40°N), the estimated σ is larger, likely due to lower interferometric phase signal-to-noise ratio caused by surface vegetation. Also, despite adequate signal-to-noise ratio in the Orange County basin, the estimated σ is relatively large, possibly caused by the temporally variable rate of VLM attributed to the aquifer depletion and recharge process during the ALOS and Sentinel-1 observation periods (figs. S4 and S5).

Figure 2 (C and D) shows the comparison between estimated VLM and the vertical rates measured at GNSS station locations. The long-term rate of vertical displacement at GNSS stations is provided by the Nevada Geodetic Laboratory (33). The stations’ measurements span time frames of 1 to 22 years, with an average time span of ~8 years. We averaged VLM rates within 100 m of each station to compare to the GNSS rates (Fig. 2C). To avoid crowding and overlapping of the stations in Fig. 2C, we show a randomly selected subset (~15%) of GNSS stations. For the bivariate comparison (Fig. 2D), we only considered 159 reliable GNSS stations with σ <5 mm/year of vertical velocity. The two datasets agree with an average SD of difference of ~0.98 mm/year. However, the spatial patterns (Fig. 2C) also show several discrepancies. Notably, the stations in Northern California have a poorer agreement to the VLM map than those in the rest of the state. This is likely due to dense vegetation degrading the phase coherence. We also see disagreement over the Orange County basin, where the VLM map indicates an uplift signal at a rate of up to 5 mm/year, while GNSS long-term rates show slow subsidence. We will further discuss this signal and argue that the observed uplift is consistent with the current status of the aquifer system in Orange County.

We highlight four low-lying metropolitan regions (Fig. 1C) that are affected by land subsidence: San Francisco Bay Area, Monterey Bay, Los Angeles, and San Diego. The VLM rate maps for each of these regions are illustrated in Fig. 3 (A to D). Note the different ranges of rates illustrated by the color bars and rate histograms. The vast majority of the San Francisco Bay perimeter is undergoing subsidence with rates reaching −5.9 mm/year, except for two zones of uplift associated with the Santa Clara and Livermore aquifers (Fig. 3A). Notably, the San Francisco International Airport is subsiding with rates faster than −2.0 mm/year. The Monterey Bay Area, including the city of Santa Cruz, is rapidly sinking without any zones of uplift (Fig. 3B). Rates of subsidence for this area reach −8.7 mm/year. The Los Angeles area shows subsidence along small coastal zones, but most of the subsidence is occurring inland (Fig. 3C). Behind and adjacent to the narrow zone of coastal subsidence, a large region of the coast is uplifting near Los Angeles. Farther inland, Los Angeles is surrounded by patchy, variable VLM rates. The variable nature of this VLM pattern could be due in part to the 2014 La Habra earthquake, which occurred on a shallow oblique fault to the northeast of Los Angeles (35), and in part due to water and hydrocarbon extraction. San Diego and the surrounding area are characterized entirely by subsidence (Fig. 3D), with rates reaching −2.7 mm/year.

The SAR datasets spanning 2007–2018 include acquisitions from Sentinel-1A/B C-band satellite in descending and ascending orbit geometry and ALOS L-band satellite in ascending orbit geometry (table S1). Figure 1B shows the footprints of different SAR frames, whose information is provided in table S1. Using these datasets, we generated several large sets of interferograms. For the Sentinel-1 interferometric dataset, the maximum temporal and perpendicular baselines are set to be 100 days and 200 m. The corresponding numbers are 3 years and 1500 m for the ALOS dataset. We also performed a multilooking factor to yield pixels of dimension ~80 m ×~80 m in the range and azimuth directions for ALOS, and ~70 m ×~70 m in the range and azimuth directions for Sentinel-1A/B. The combination of both the ALOS and Sentinel-1 datasets is necessary because we are interested in determining the multidecadal rate; broader temporal spread of data leads to a more robust outcome. Including the gaps, the combined time frame for ALOS and Sentinel-1 covers nearly 11 years.

Each dataset was processed separately to generate three accurate and high-resolution time series of surface deformation. To this end, we applied the Wavelet-based InSAR algorithm (51, 52). The geometrical phase was calculated and removed using a 30-m digital elevation model provided by the Shuttle Radar Topography Mission (53) (www2.jpl.nasa.gov/srtm/) and the satellite precise ephemeris data (54). We calculated the time series of complex interferometric phase noise in the wavelet domain and analyzed in a statistical framework to identify elite (i.e., less noisy) pixels. Next, we discarded the nonelite pixels and obtained the absolute estimate of the phase change for elite pixels via an iterative 2D sparse phase unwrapping algorithm. We corrected each unwrapped interferogram for the effect of orbital error (55) and the topography-correlated component of atmospheric delay (56). Using a robust regression, we inverted the phase changes from the large set of interferograms. This approach reduces the effects of outliers due to improper phase unwrapping. Using continuous wavelet transforms, we applied a high-pass filter to reduce the temporal component of the atmospheric delay as well as seasonal variations in the ground surface deformation due to local hydrology. Therefore, the deformation contributions with periods less than 1 year were removed. This operation is a reasonable approach because we are only interested in long-term deformation signals. Using a reweighted least-squares estimation, we calculated velocities along the LOS direction for each elite pixel as the slope of the best fitting line to the associated time series. Last, we geocoded all datasets to obtain precise locations of elite pixels in a geographic reference frame. Figures S6 to S8 show the LOS velocity associated with each frame.

Once the accurate LOS velocity of each frame was obtained, we followed the procedure of Ojha et al. (57) to mosaic SAR frames that belong to the same satellite track of a sensor, generating three large scale maps of LOS displacements in Sentinel-1 ascending and descending as well as ALOS ascending geometry. To this end, we used pixels within overlapping areas of adjacent frames as tie points and applied an affine transformation that includes a constant translation and two rotations. We determined the parameters of this transformation using measurements of collocated pixels within the overlap zone. Applying this transformation to the second frame aligns its LOS velocity with that of the first one. This procedure yields a seamless map of the LOS velocity field as long as the contributions from atmospheric errors are appropriately corrected. We removed the extraneous pixels caused by the overlap to avoid double counting in the intersection. The newly joined frames then became the starting frame for the next neighbor. We continued this process until all frames of a given track were oriented with respect to the starting frame.

We repeated this process for each of the three satellite tracks, thus creating three spatially continuous LOS velocity maps for the coast of California. To reduce error propagation from the compilation step, mitigate long-wavelength errors due to ionospheric and residual orbital errors, and transform LOS velocities from a local to a global stable reference frame, we used the displacement velocities of 441 GNSS stations located within our study area (33). Using the local incidence and satellite heading angles, we projected the 3D displacement velocities from GNSS onto the satellite LOS direction (58). Next, we applied the affine transformation to align the spatially continuous maps to the GNSS observations projected on LOS. At this step, LOS velocities associated with each viewing geometry were referred to a global reference frame, here NA12. This procedure is similar to that used for aligning individual GNSS solutions to the reference frame when adjusting large GNSS networks. Figures S1 to S3 show maps of LOS velocities, which are a seamless combination of the individual frames.

Next, we implemented a Kriging interpolation approach with inverse distance weighting (59) to interpolate the LOS velocities of Sentinel-1 tracks and GNSS horizontal velocities on the location of pixels within the ALOS track. Thus, for each pixel, we obtained three LOS observations and two GNSS horizontal velocities (E-W and N-S components), all in the NA12 frame.

Let {y_sa, y_sd, y_aa} and {σ²_sa, σ²_sd, σ²_aa} be the interpolated LOS velocities and variances for a given pixel, where subscripts indicate Sentinel-1 ascending (sa), Sentinel-1 descending (sd), and ALOS ascending (aa). The stochastic model to combine these three LOS velocities with horizontal velocities of GNSS datasets to generate a seamless, high-resolution, and accurate map of east (E), north (N), and vertical (U) motions is given bywhere C represents the unit vectors projecting the 3D displacements onto the LOS (58), E_G and N_G are the interpolated GNSS velocities in the E-W and N-S directions, and ε is the observations error with N(0, σ) probability distribution (σ is the SD). Note that in Eq. 1, E and N are unknown while E_G and N_G are observed. This kind of stochastic model is known as a unified least-squares adjustment because it lifts the barrier between observation and parameters (60). We assigned a σ of 5 mm/year for the LOS velocities and used the σ provided by Nevada Geodetic Laboratory (33) for the GNSS E-W and N-S components. The weight matrix is inversely proportional to σ², and the solution to Eq. 1 is given bywhere

The parameters variance-covariance matrix is $Q_{\mathrm{XX}}={\mathrm{\sigma }}_0^2{(A^T\mathrm{PA})}^{-1}$, where ${\mathrm{\sigma }}_0^2=\frac{r^T\mathit{\text{Pr}}}{n-u}$, r = L − AX, n is the number of equations, and u is the number of unknowns. By including both ALOS and Sentinel-1 datasets, we increase our redundancy, thus refining errors. Figures S13 to S15 show the residual maps for the three LOS velocity maps. Systematic patterns with opposite signs in ALOS and Sentinel-1 descending residual maps are likely due to remaining orbital and ionospheric errors, which are adjusted when applying Eq. 1 to solve for the 3D velocity field. This is an advantage of combining data from different viewing geometries, which further highlights importance including ALOS data in our analysis. Once the 3D displacement for each pixel was calculated, we performed an affine transformation to realign further the vertical velocities obtained through Eq. 2 with the stable North America reference frame (NA12). We used the vertical velocities of GNSS stations that lie in the area of InSAR coverage, provided by the Nevada Geodetic Laboratory (33), shown in Fig. 1B. The resulting VLM rate and SDs are shown in Fig. 2 (A and B). The comparison against the GNSS vertical observation is shown in Fig. 2C. We obtained an SD of 0.98 mm/year for the difference between two datasets (Fig. 2D).