Research ArticleGEOLOGY

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

See allHide authors and affiliations

Science Advances  31 Jul 2020:
Vol. 6, no. 31, eaba4551
DOI: 10.1126/sciadv.aba4551

Abstract

Coastal vertical land motion affects projections of sea-level rise, and subsidence exacerbates flooding hazards. Along the ~1350-km California coastline, records of high-resolution vertical land motion rates are scarce due to sparse instrumentation, and hazards to coastal communities are underestimated. Here, we considered a ~100-km-wide swath of land along California’s coast and performed a multitemporal interferometric synthetic aperture radar (InSAR) analysis of large datasets, obtaining estimates of vertical land motion rates for California’s entire coast at ~100-m dimensions—a ~1000-fold resolution improvement to the previous record. We estimate between 4.3 million and 8.7 million people in California’s coastal communities, including 460,000 to 805,000 in San Francisco, 8000 to 2,300,00 in Los Angeles, and 2,000,000 to 2,300,000 in San Diego, are exposed to subsidence. The unprecedented detail and submillimeter accuracy resolved in our vertical land motion dataset can transform the analysis of natural and anthropogenic changes in relative sea-level and associated hazards.

INTRODUCTION

A substantial proportion of the world’s population lives in low-lying coastal areas (1). Of these, nearly 1 billion live in flood-prone areas lying below 10-m elevation (2). Coastal populations are expected to increase by the year 2050 due to coastward migration (3). Moreover, these coastal lands are subject to subsidence due to natural and anthropogenic processes (4). Land subsidence can increase flooding risk (5), wetland loss (6), and damage to infrastructure (5, 7) by increasing a region’s relative sea-level rise (RSLR) (8). In a coastal setting, such as California’s San Francisco Bay Area (9), low-elevation, highly populated coastal cities experiencing subsidence have an increased hazard of flooding, inundation, and related economic damages—especially if the area is financially disadvantaged, lacks proper infrastructure, or does not have adequate monitoring efforts in place (10). Furthermore, subsiding coastal wetlands, the U.S. Gulf Coast for example, are at increased risk of drowning and disappearing due to ocean flooding (11).

Natural drivers of subsidence and uplift, such as tectonics (1214) [except coseismic events (14)], isostasy (15), and sediment compaction (16, 17), often cause slow monotonic vertical land motion (VLM). Conversely, the VLM associated with anthropogenic processes such as surface water management or extraction of groundwater and hydrocarbons (6, 8, 18, 19) can be fast with temporally variable behavior. In some cases, depletion and recharge of aquifers can cause VLM at rates up to tens of centimeters per year VLM [e.g., (20, 21)]. Depending on the local geology, these processes can combine to generate overall uplift or subsidence in a landscape.

In California, the oblique convergence between the Pacific and North American plates is primarily accommodated by strike-slip motion along major faults such as the San Andreas fault system (22). However, the minor dip-slip component of these near-vertical faults can cause spatially and temporally variable VLM during earthquake cycles (22). Farther north, the migration of the Mendocino Triple Junction drives crustal thickening and uplift north of the junction and crustal thinning and subsidence to the south (23). California experiences slow subsidence due to glacial isostatic adjustment with rates up to ~2 mm/year at latitudes >37°N and up to ~1.5 mm/year at lower latitudes (24). Sediment compaction leads to subsidence in sedimentary basins, such as in the San Francisco Bay (9, 16), and droughts and groundwater withdrawal lead to subsidence in aquifer basins, such as in the Santa Clara Valley (2527). These combined effects significantly alter regional RSLR rates along California’s coast (9).

Despite the importance of monitoring land subsidence, the large-scale, high-resolution observations that would help characterize VLM and its related hazards are often not available. Instead, understanding of VLM in California is limited to sparsely distributed (~20 km average spacing) Global Navigation Satellite System (GNSS) stations. This coarse GNSS resolution reduces the accuracy of RSLR estimates; the lack of colocated GNSS stations with tide gauges leads to errors in estimates of global mean sea-level rise (28).

Spaceborne geodesy through combining observations obtained via interferometric synthetic aperture radar (InSAR) with those measured by GNSS stations offers a means to improve the accuracy, density, and spatial extent of VLM measurements (13, 26, 2931). Here, we present the first VLM rate map for the coast of California at ~100-m resolution and unprecedented precision (~1 mm/year) as a combination of the full archive of Advanced Land Observing Satellite (ALOS) L-band and Sentinel-1A/B C-band satellites with observations of horizontal velocities at GNSS stations. The ALOS archive covers 2007–2011, Sentinel-1A/B covers 2014–2018, and the GNSS data cover 1996–2018, with individual stations covering time frames ranging from 1 to 22 years and a mean observation time of 8 years. This map, composed of ~35 million pixels, shows a wide range of uplift and subsidence rates with both large- and small-scale patterns. We discuss these results in the context of the underlying drivers of VLM in California, exploring both the natural and anthropogenic effects. We also discuss the consequences of coastal land subsidence in California.

The results presented here, in combination with other hydrogeological measurements, can contribute to understanding rates of discharge and recharge and the elastic and inelastic responses of aquifers to groundwater exploitation (6, 8, 19). It also can elucidate the effects of fossil fuel extraction and wastewater injection in regions of gas and oil production (6, 8). Furthermore, this VLM rate map provides valuable kinematic context to the overall understanding of California’s geology and tectonics (13, 14). Last, this study’s methodology is a cost-effective and time-efficient process for determining highly accurate VLM rates. The methods presented here can be applied to any region with adequate SAR and GNSS data coverage. The dataset of VLM rates generated here is available for download in community-wide subsets based on community outlines defined by the 2010 U.S. Census (see the Supplementary Materials).

RESULTS

Datasets and analysis

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.

Fig. 1 Study area and datasets.

(A) Location of California. (B) The interferometric datasets include acquisition from ALOS L-band (yellow rectangles) spanning 2007–2011 and acquisitions obtained in descending (green rectangles) and ascending (blue rectangles) orbits of Sentinel-1A/B satellites during 2014–2018. Table S1 summarizes the path and frame parameters of the SAR datasets. The GNSS dataset (black triangles) includes 811 stations with between 1 and 22 years of observation (and an average observation duration of ~8 years) (source: http://geodesy.unr.edu/). California’s active faults are shown in red (source: https://earthquake.usgs.gov/hazards/qfaults/). (C) Coastal elevation in California. Coastal zones, which are defined to be those with elevations less than 10 m, are shown in red. Segments of the coast with elevations higher than 10 m are colored by a yellow gradient (source: USGS NED).

VLM along the coast of California

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).

Fig. 2 VLM along the coast of California during 2007–2018.

(A) Estimated VLM rates in NA12 reference frame as a combination of SAR datasets and GNSS measurements. The map includes ~35 million pixels at ~100 m × ~ 100 m resolution. Black boxes indicate extent of study areas in Fig. 3 (A to D). (B) Standard deviation for each pixel of the estimated VLM. (C) Validating estimated VLM rates (open circles) against the GNSS vertical rates (filled circles). For this comparison, the observations of pixels within 100 m of each GNSS station are averaged. Transparent circles represent rates of 0 mm/year. (D) Bivariate plot comparing GNSS vertical rates (red) with estimated VLM rates (blue). Standard deviation of the difference between two datasets is 0.98 mm/year. Bars show SD (1-σ errors) for each dataset.

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).

Validation test and accuracy

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.

Major subsiding cities

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.

Fig. 3 VLMs for the four subregions during 2007–2018.

(A) San Francisco Bay Area, (B) Monterey Bay, (C) Los Angeles, and (D) San Diego. The location of each subregion is indicated in Fig. 2A. Note the different rate scales for each location. Overlain diamonds show vertical rates measured by local GNSS stations. Corresponding histograms show distributions of VLM rates for each panel with quarter millimeter bins. Note the differences in spread. The magenta vertical line marks a velocity of 0 mm/year.

DISCUSSION

We combined three extensive SAR datasets from different viewing geometries of L-band and C-band satellites with measurements of horizontal displacements rates at GNSS stations to create the first spatially continuous vertical deformation map for the full coast of California at millimeter accuracy and gridded at ~100-m spatial resolution. The results show VLM patterns that align with the known geologic and anthropogenic mechanisms acting on this landscape. Of particular importance is the subsidence underlying several large cities in Central and Southern California. In the following section, we discuss the limitations and implications of this study’s results.

Here, we assumed that the VLM rates are steady over the observation period of 2007–2018. This period is short and not continuous due to the gap between 2011 and 2014 between the two satellite acquisition periods. Short-term signals during the two time frames are interpreted as long-term geologic trends. The comparison with GNSS long-term vertical rates (Fig. 2, C and D) indicates that our assumption is reasonable, as long as the signal is not contaminated by transient effects associated with aquifer and reservoir compaction.

Forward projection of these rates is critical for their accurate application to predictions of future land subsidence. However, the question remains as to which functional form should be used to obtain a robust estimate of future VLM and uncertainties. The VLM rate due to glacial isostatic adjustment can be considered steady over a century (36). Also, the contribution from soil compaction, depending on the history, net deposition rate, and local geology, can be regarded as near steady (16). Thus, the contemporary rates of VLM associated with these two causes can be projected forward using a linear model, and the errors remain negligible [e.g., (9)]. However, compaction of hydrocarbon reservoirs is a function of demand driven by the price, and it may increase as oil prices rise (37). Also, land subsidence due to groundwater extraction is driven by demand for fresh water and recharge and is likely to increase during the 21st century as the frequency of the meteorological (less rain) and agricultural (less soil moisture) drought increases (38). The contributions in VLM due to reservoir and aquifer compaction are often nonlinear, and models capturing various physical and socioeconomic parameters are needed to obtain accurate VLM projections. In future work, we plan to explore this question further and to develop a projection of VLM and its associated implications for RSLR and storm surge flooding for the year 2100 for the study area.

Several highly populated areas experience subsidence or variable subsidence and uplift (Fig. 3). The subsidence around the banks of the San Francisco Bay (Fig. 3A) is largely influenced by the tectonics of the region, primarily the San Andreas Fault (13). It also is influenced by sediment compaction of the landfill and Bay Mud deposits that comprise much of the underlying ground for locations such as the San Francisco International Airport and Redwood City (39). The uplift signals in the Santa Clara and Livermore valleys indicate the response of aquifers to recharge processes (25, 26). However, this uplift is likely an overestimate of the rebound from aquifer recharge as the satellite acquisition time frames do not cover the entirety of the 2012–2015 California drought. The InSAR time frames considered here miss the years 2011–2014, during which the Santa Clara Valley aquifer system was recorded as subsiding before beginning its rebound in late 2014 (27). Apart from these two zones of uplift, the overall sense of subsidence matches results from Hammond et al. (40), which were derived from Global Positioning System (GPS) data covering 1996–2016, encompassing the drought period. This discrepancy in the Santa Clara and Livermore basins further supports the idea that these uplift signals are associated with aquifer recharge.

The regional variability of VLM in and around Los Angeles (Fig. 3C) is congruent with the degree of spatially variable VLM measured by Brooks et al. (29) and reflects the region’s extensive groundwater usage and fossil fuel extraction. Recharge of local aquifers likely accounts for the expansive zone of rapid uplift surrounding Santa Ana (41). Figures S4 and S5 show the regional groundwater trend as recorded by the California Department of Water Resources wells at 100- to 600-m depth in the Coastal Plain of Orange County aquifer system during the satellite data acquisition time frames (2007–2011 and 2014–2018 for ALOS and Sentinel-1, respectively). Several of the wells during 2007–2011 and most of the wells during 2014–2018 show increases in groundwater levels during these time frames, which agrees with the region of uplift encompassing this area (Fig. 3C and figs. S4 and S5). Again, disagreement with other long-term vertical rates derived for this region (34, 40, 42, 43) reenforces that this uplift is indeed an aquifer recharge signal. Tectonics also likely contribute to this region’s deformation pattern. The shearing, contraction, and rotation of the Western Transverse Ranges have been found to cause basin subsidence, as in the Ventura basin, where faulting occurs at rates on the order of millimeters per year (44, 45). Conversely, transtensional faulting patterns have also been found to cause subsidence, as in the San Bernardino basin, where interseismic strain results in subsidence at rates up to 0.3 mm/year (43). The superposition of these faulting signals with anthropogenic drivers likely accounts for the complexity of the VLM signal in the greater Los Angeles region.

In the Monterey Bay region (Fig. 3B), the likely primary drivers of this subsidence are sediment compaction of the Salinas, Elkhorn, and Pajaro River sediments and tectonic activity along the Monterey Bay Fault Zone and the nearby San Andreas Fault (46). The river basins’ subsidence is in disagreement with the results of Hammond et al. (40), which shows uplift in the center of Monterey Bay. This difference could be groundwater related, but the remediation of this disagreement will require further analysis of regional well data. In San Diego (Fig. 3D), the coastline is composed primarily of quaternary alluvial sediments except for the Point Loma peninsula, which lies on Upper Cretaceous marine sedimentary units. Farther inland, the geology switches to both marine and nonmarine sedimentary units from the Pliocene and Eocene, respectively (47). Compaction of the quaternary alluvium may be responsible for the subsidence trend for this part of the coast. The widespread subsidence observed in the San Diego area agrees with GPS-derived results for this region (40).

Subsidence poses risks to infrastructure and ecosystems (7, 18). Figure 4 illustrates subsidence on a by-county basis (Fig. 4A), the smaller communities and their percentages of area affected by subsidence (Fig. 4B), and estimates of the number of people (assuming uniform population distribution) in each community that will be exposed to subsidence (Fig. 4C). We find that the communities of San Francisco, Los Angeles, and San Diego each experience large percentages of subsiding area. These three communities also have the highest estimated population exposure to subsidence, with likely ranges (67% probability) of ~460,000 to 805,000, ~8500 to 2,299,000, and ~202,000 to 2,259,000 people in the San Francisco, Los Angeles, and San Diego communities, respectively. In total, we estimate that between ~4.3 million and 8.7 million people in California’s coastal communities are exposed to subsidence (table S2). This exposure analysis does not consider proximity to water bodies, the presence of infrastructure, or other factors that make land subsidence hazardous. Rather, it serves as a population-based guide for which coastal communities are most in need of more detailed risk assessments.

Fig. 4 Exposure to land subsidence along the coast of California during 2007–2018.

(A) VLM rates are saturated to emphasize negative values, indicating land subsidence. Black polygons map California’s counties. (B) Percentage of each communities’ area (black polygons) affected by subsidence. California community shape outlines are defined by the 2010 U.S. Census and are subsections of the greater California counties. (C) Population affected by subsidence for each California community (black polygons) and the three communities with the highest estimated exposed population: San Francisco, Los Angeles, and San Diego.

RSLR, which is of importance to low-lying coastal communities and ecosystems, often deviates from the regional estimates of sea-level rise (48). This study, by measuring contemporary land subsidence rates at an unprecedented resolution, shows that much of California’s coast is subsiding. This will cause increases in RSLR rates for those sections of the coastline. The VLM will likely accelerate throughout the 21st century due to compaction of aquifers, reservoirs, and soil associated with increased groundwater extraction linked to intensified and frequent droughts (49), demand for fossil energy (37), and coastal land reclamation to relieve the pressure from coastward population migration (3). Without consideration of VLM, future projections of SLR, including those considering extreme scenarios with higher rates of Antarctica and Greenland ice loss during the second half of the 21st century (50), may be underestimating the future inundation hazards to coastal communities. This study emphasizes the urgency with which the flood resiliency plans must adapt to scenarios in which coastal land elevation drops rapidly. Understanding the key processes driving relative sea-level change enables policy makers to prioritize the risk reduction and adaptation interventions and to better identify the communities and ecosystems most vulnerable to flooding.

METHODS

SAR datasets

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.

SAR interferometric analysis

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.

SAR frame compilation

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.

3D velocity inversion

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 {ysa, ysd, yaa} and {σ2sa, σ2sd, σ2aa} 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 byysa=CesaE+CnsaN+CusaU+εsaysd=CesdE+CnsdN+CusdU+εsdyaa=CeaaE+CnaaN+CuaaU+εaaEG=E+εeNG=N+εn(1)where C represents the unit vectors projecting the 3D displacements onto the LOS (58), EG and NG 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 EG and NG 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 σ2, and the solution to Eq. 1 is given byX=(ATPA)1ATPL(2)whereX=(E N U)T,A=(CesaCnsaCusaCesdCnsdCusdCeaaCnsdCuaa100010),P=(1(σsa)2000001(σsd)2000001(σaa)2000001(σe)2000001(σn)2)

The parameters variance-covariance matrix is QXX=σ02(ATPA)1, where σ02=rTPrnu, r = LAX, 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).

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/31/eaba4551/DC1

https://creativecommons.org/licenses/by-nc/4.0/

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: M.S. and E.B. are supported by the National Aeronautics and Space Administration (NASA) grant 80NSSC170567. C.O. and S.W. are supported by the National Aeronautics and Space Administration grant NNX17AD98G. Author contributions: All authors made significant and extensive contributions to the work presented in this manuscript. E.B. and M.S. designed the experiment. M.S. and C.O. processed the SAR dataset and performed the SAR interferometric analysis. E.B., S.W., and M.S. analyzed the groundwater data. E.B. and M.S. compiled the SAR frames, performed the 3D inversion analysis, derived the exposure metrics, and wrote the paper. E.B. created the figures and the Supplementary Materials. All authors reviewed and edited the manuscript, figures, and the Supplementary Materials. Competing interests: The authors declare that they have no competing interests. Data and materials availability: GNSS velocities are provided by Nevada Geodetic Lab (http://geodesy.unr.edu/). SAR datasets are obtained from the Alaska Satellite Facilities (www.asf.alaska.edu). U.S. Geological Survey provided LiDAR data shown in Fig. 1B. The VLM data necessary to evaluate the conclusions in the paper are available for download on a community-wide basis at the following online repository: https://drive.google.com/drive/folders/1NUpEiRjGTVXwZr61p2K-GPuV_0A6pJnz?usp=sharing. The format is further explained in the Supplementary Materials. 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

Stay Connected to Science Advances

Navigate This Article