Time-varying land subsidence detected by radar altimetry: California, Taiwan and north China

Contemporary applications of radar altimetry include sea-level rise, ocean circulation, marine gravity, and icesheet elevation change. Unlike InSAR and GNSS, which are widely used to map surface deformation, altimetry is neither reliant on highly temporally-correlated ground features nor as limited by the available spatial coverage, and can provide long-term temporal subsidence monitoring capability. Here we use multi-mission radar altimetry with an approximately 23 year data-span to quantify land subsidence in cropland areas. Subsidence rates from TOPEX/POSEIDON, JASON-1, ENVISAT, and JASON-2 during 1992–2015 show time-varying trends with respect to displacement over time in California’s San Joaquin Valley and central Taiwan, possibly related to changes in land use, climatic conditions (drought) and regulatory measures affecting groundwater use. Near Hanford, California, subsidence rates reach 18 cm yr−1 with a cumulative subsidence of 206 cm, which potentially could adversely affect operations of the planned California High-Speed Rail. The maximum subsidence rate in central Taiwan is 8 cm yr−1. Radar altimetry also reveals time-varying subsidence in the North China Plain consistent with the declines of groundwater storage and existing water infrastructure detected by the Gravity Recovery And Climate Experiment (GRACE) satellites, with rates reaching 20 cm yr−1 and cumulative subsidence as much as 155 cm.

(NCP). The three croplands areas are affected by historical and ongoing subsidence 3,[11][12][13][14][15][16] , where groundwater extraction is the leading cause of subsidence 2,6,12,16 , and recent droughts have exacerbated groundwater depletion and its consequences in the SJV and the NCP 13,17 . Subsidence poses risks to existing high speed railways in CT, the NCP, and one under construction in the SJV 4 . Subsidence patterns and rates in much of the NCP, SJV and CT have been extensively studied using precision leveling, extensometry, InSAR, and GNSS 4,6,14,15,18 .
A synthesis of vertical displacement rates (VDRs) measured at 1,499 GNSS stations (Table S1 and Fig. S1, Sec. S1, Supplementary Information [SI]) between 1993 and 2015 across the Central Valley (including the SJV and the southern portion of the Sacramento Valley), California, show that the distribution of the GNSS sites within the valley are too sparse by themselves to map areal subsidence in adequate spatial detail (Fig. 1). Three coalescing subsidence bowls are evident: a southern bowl south of Bakersfield corresponding to the Arvin-Maricopa historical subsidence bowl 19 , a central bowl between Bakersfield and Fresno roughly centered on the Tulare-Hanford area and roughly corresponding to the Tulare-Wasco historical subsidence bowl 19 , and a northern bowl north of Fresno roughly centered on the Madera-Mendota area near the Delta-Mendota Canal 15 , which is located along the western edge of the San Joaquin Valley (Fig. 1). In general, the GNSS-derived subsidence rates are relatively large in the Tulare-Hanford area (maximum: 9.8 cm yr −1 ), and the rates decay southward and northward in the SJV to near zero. In the mountainous regions surrounding the Central Valley, the VDRs are largely positive, suggesting uplifts possibly caused by tectonics and reduced mass loading due to the current drought in California 20 . The mean subsidence rate in the central subsidence bowl (from Bakersfield to Fresno) is 2.86 cm yr −1 , compared with 0.62 cm yr −1 in the northern subsidence bowl and smaller mean rates in the Sacramento Valley. Note that a negative vertical displacement rate determined in this paper (in all relevant figures and tables) corresponds to a positive subsidence rate. To emphasize land subsidence in this paper, we replace negative displacement rates by subsidence rates in all descriptions.
Mean subsidence rates during 2002-14 in CT derived from annual precision leveling ranged from zero to about 8 cm yr −1 (Fig. 2a, Sec. S2, SI). Despite its high accuracy, precision leveling can be costly 4 and generally is not repeated frequently enough to resolve inter-annual and annual signals. About 310,000 pumping wells in Changhua and Yunlin Counties pump groundwater to meet agricultural and industrial water demands (Fig. 2b). Two major coalescing subsidence bowls are evident-on either side of the Zhuoshui River, and subsidence rates as much as 7 cm yr −1 threaten the operation of the Taiwan High Speed Rail 4 (THSR; Fig. 2a) and sustainable development in the region. For example, the prevalence of fine-grained sediments here makes this area more susceptible to compaction. Though pumping-well distribution is even, subsidence rates are not, suggesting a significant spatial variation of hydrogeological properties 6 .
The NCP has a population of 437 million and an area of 300,000 km 2 , and is one of many regions in China with critical subsidence problems 16,21 . In the NCP, groundwater provides more than 60% of fresh water supplies 22 and is the main water source for agricultural irrigation 23 . Near Tianjin, the mean subsidence rates range from 0.80 to 5.60 cm yr −1 , with a maximum rate of 16 cm yr −1 and cumulative subsidence of 3.90 m during 1965-85 18,24 . The GRACE satellite gravimetry data have detected mass losses that are attributed to severe groundwater storage declines in the SJV and NCP 13,22 .
In this paper, we explore satellite altimetry as a new remote-sensing subsidence-mapping method and demonstrate its utility in cropland areas where the differential and persistent scatterer (PS)-based InSAR methods are limited by sparse temporally-coherent, stable radar reflectors 3 ; however, improvements in identifying PS scatterers over agricultural regions have been reported 25,26 . We also demonstrate the temporal and along-track spatial detail of the radar altimetry method compared to GNSS (Fig. S1) and leveling (Fig. S2) which can be cost prohibitive to map subsidence in similar temporal and spatial detail. Satellite altimetry has been principally used to derive marine gravity/bathymetry, and to study ocean circulations, sea level rise, ice sheet elevation and water level changes 27-29 . Because only a few studies have used satellite altimetry to measure solid earth deformation and a standard numerical procedure to measure vertical displacement has not yet been established, there are concerns about the accuracies of altimeter-derived rates 30,31 . However, this study shows that, accurate VDRs at the sub-cm yr −1 level can be achieved using satellite altimetry by using more than two decades of near-continuous land surface elevation measurements from the improved coverage offered from multi-mission altimeter data with 10-day or monthly sampling: TOPEX/POSEIDON (TP, 1992(TP, -2002

Results
Subsidence rates derived from satellite radar altimetry (TP and EN) agree well with rates determined from precision leveling in CT and GNSS in the southern SJV (Fig. 3). Height changes derived from different radar altimetry missions also agree well. For example, in the SJV at the crossover of TP-043 and EN-684 (Fig. 3a), the time-series height changes derived from J1 (the follow-on mission for TP) and EN (Fig. S6, 2002-10) and the resulting subsidence rates (about 6 cm yr −1 ) show good agreement. Three major subsiding bowls identified in the SJV from the GNSS measurements ( Fig. 1) are reflected in the TP and EN along-track subsidence profiles. Along TP-043, the largest altimeter-derived subsidence rate near Hanford (A2, Fig. 3a) in the central subsidence bowl is about 10 cm yr −1 , and consistent with the GNSS-derived rates (Table S1). At GNSS stations P056 and P566, subsidence rates are <4 cm yr −1 and agree with the altimeter-derived rates to 1 cm yr −1 . The largest subsidence rate from EN is ~13 cm yr −1 along EN-611 (A1, Fig. 3b). In CT, the two principal subsidence bowls identified by leveling (Fig. 2a) are reflected in the TP and EN along-track subsidence profiles, and the patterns and rates are consistent with those from leveling ( Fig. 3d-f).
Subsidence patterns in the SJV and CT are time-varying and could be influenced by changing groundwater use related to land-use, climatic and regulatory factors. The rates in CT also vary with time and space. In Box 4 (Fig. 4b), the subsidence rates increased in recent years. In Box 5 (Fig. 4b), between 23.60°N and 23.75°N, the changes in the subsidence rates from TP to J2 show the extent of subsidence shrinks (from green lines to red lines) in recent years. The location of maximum subsidence shifts northward to about 23.72°N and the extent of subsidence increases northward from 23.71 to 23.85°N. The shifts in patterns likely are due to measures restricting groundwater use along a section of THSR near the point of maximum subsidence (Figs S7 and S8).
Land subsidence can cause angular deflection (AD) of the THSR owing to vertical displacements of the individual structural support columns of the elevated railway. A large AD could weaken the foundational support of the railway and result in serious operational problems if not adequately mitigated 4 . Along the THSR section in Yunlin County (Fig. 2a), heavy groundwater pumping and heterogeneity of sediments with varying degrees of compressibility have led to large ADs, which once approached the 1/1000 limit set by the safety code. TP-051 is nearly parallel to THSR near Tuku (see Fig. S7) and altimetry data have already provided valuable observations to monitor ADs here. The planned California High-Speed Rail (CHSR) now faces the same potential risk of large ADs as the TSHR does: it passes through a region that has been experiencing subsidence due to heavy groundwater depletion. As it happens, EN-684 and TP-043 are perpendicular to the planned route of CHSR near Hanford and the observations from EN and TP follow-on missions can provide an effective means to monitor the subsidence as well as the ADs in the future.
EN identifies two regions of major cropland subsidence in the NCP (Fig. 6, also Figs S9-12): Region 1 covers eastern Hebei Province and northern Shandong Province, with an affected area of 67,900 km 2 ; Region 2 covers eastern Henan Province and northern Anhui Province, with an affected area of 78,000 km 2 . In the foothills of the NCP, VDRs derived from EN and TP, and J2 are indeterminate because the altimetry waveforms are corrupted (see Waveform E, Fig. 7). The J1 rates contain large uncertainties in much of the NCP and are not shown here. TP and J2 detect time-varying subsidence elsewhere in the NCP (Fig. 8). For example, in Boxes 1 and 2 (Fig. 8), located in the Hebei and Anhui provinces, the along-track rates increased dramatically from the time of TP to the time of J2. In Box 3 (Fig. 8), near Tianjin, the J2 subsidence rates are less than those from TP, suggesting that groundwater extraction was reduced here due to measures implemented to control subsidence since the start of high-speed rail operations (Fig. S9). In Box 4 (Fig. 8), many VDRs were derived from TP, but only a few from J2. This is attributed to urbanization around this region, where new man-made structures (at the time of J2) contaminated altimeter waveforms and degraded ranging accuracies (see Waveform D, Fig. 7). The altimeter results identify many subsidence-affected areas and changing subsidence patterns in the NCP that are consistent with the temporal and spatial groundwater storage changes detected by GRACE 22 . The results can be used as guides for future, more detailed, precise geodetic and geotechnical measurements.

Discussion
The causes for subsidence are complicated and not well understood in the study regions (SJV, CT and NCP). Groundwater abstraction may be a major contributing source for land subsidence. Near the point with the largest subsidence in SJV, operations in factory and cattle farms, such as those shown in Box 1 and 2 in Fig. 4a, may contribute partially to the steep subsidence observed along the profiles through abstraction of large amounts of groundwater. Along the THSR in CT, heavy groundwater extraction from surrounding factories, such as a tire factory (Fig. 4b), may be responsible for land subsidence. However, more studies are needed to investigate the roles of the farms and factories on the land subsidence detected in this paper.
Long-term, time-varying land subsidence over croplands from 1992 to present has been determined by TP, J1, J2 and EN altimeters using a dedicated processing method for land altimeter data. The result in the NCP suggests that a mission with a long repeat period and relatively small cross-track spacing (< 50 km) like EN can measure cropland subsidence over a large area. The application of altimetry to subsidence monitoring can be extended globally to other croplands, especially in areas where terrestrial-geodetic subsidence monitoring methods may be lacking. The GNSS and leveling-derived vertical displacement rates in the SJV and CT are used to assess the accuracies of satellite altimetry, which are better than the one cm yr −1 level. Current altimetry technologies can detect subsidence at a 1-km along-track spatial resolution (c.f. Methods for details), and likely will improve as the technologies evolve. Combined use of leveling, GNSS, InSAR and radar altimetry methods for measuring, mapping and monitoring surface deformation and subsidence in particular, would optimally leverage the strengths of each of the different techniques to best address the scientific needs.

Methods
A satellite-borne radar altimeter measures the two-way time delay of the radar pulse which is used to compute the distance between the satellite and the Earth's surface, and after computing the precise orbit of the satellite, the altitude of the surface can be computed. Repeat-period altimetry can be used to measure changes in altitude (height) or height changes of the surface. The ground tracks are re-visited by the satellite at regular time intervals of weeks, and within about ± 1 km of the reference track at the equator. Figure 7a shows the altimeter waveforms from the Pacific Ocean off the coast of California eastward to the SJV and the Sierra Nevada. The same evolution of waveforms is also demonstrated for the NCP (Fig. 7b). The Brown waveforms over the ocean (Waveform A, Fig. 7) and coast (Waveform B, Fig. 7) may result in precise range measurements for marine applications, but waveforms on land may be too contaminated to recover precise ranges for detecting height changes.
The waveform over flat terrain with crops (Waveform C, Fig. 7) or with snowpack, e.g., land regions near Hudson Bay, Canada 29-31 can be retracked to yield a precise range change, after correcting surface gradients using the collinear track analysis. The surface roughness of croplands diffuses radar pulses, similar to wind-induced, small-scale waves in oceans. The surface of a fallow, flat and building-free cropland is similar to a calm lake surface-the altimeter waveform is specular with a steep leading edge compared to other types of waveform. In summer, the cropland surface roughens as crops grow, resulting in a waveform with a less steep leading edge. Similar to a significant wave height at sea, terrain undulation creates a ranging bias and modulates the footprint size of a pulse-limited radar. For example, a terrain undulation of 1 m at a length scale of several km leads to a pulse-limited footprint radius of about 1 km for EN and slightly larger for the TP series altimeters 28,29 (pulse width = 3.125 ns). However, unlike significant wave height at sea, terrain undulation at a given location can remain constant for the satellite mission period, and the undulation-induced ranging bias will be eliminated when heights from repeat cycles are differenced.
Here we adopted a modified, enhanced version of the repeat-track method 32-34 for the space and time reduction of raw height measurements from TP, J1, J2 and EN. For EN, we also estimated parameters that account for the effects of backscatter and the slopes of the leading and trailing edges 32,35 . For all missions, we chose to compute VDRs for bins spaced at about 1 km along satellite ground tracks. A bin is a circular region with a given radius, centered at a location along the mean ground track of all repeat cycles in a satellite mission. Over a flat cropland with a moderate undulation of about 1 m, the radius of an altimeter footprint is about 1 km. Within a given bin, all the raw height measurements (the sampling rate is 18 Hz for EN, and 10 Hz for TP, J1 and J2) from the repeat cycles were least-squares fitted to the following space and time function 3 : where b l , and τ are the means of b, l and τ from all repeat cycles in the bin, and C b ,C l and C t are the parameters that adjust these effects. For EN, up to 12 parameters could be estimated (see Table S2). We used a robust least-squares estimator, which is resilient against outliers, to estimate the parameters in Eq. (1) as follows. The first-round estimated parameters used all height measurements. Then, the residuals of the measurements and the standard deviations of the residuals were computed. We then used the three-σ outlier threshold to remove anomalous heights: if a residual v i j is three times larger than its standard deviation, the corresponding measurement H i j was removed, and the next-round parameter estimation was carried out. The final parameters were estimated from the measurements that pass the three-σ testing. Also, if the standard deviation of the residuals in a bin were larger than 5 m, the vertical rate in the bin was not computed. Then, the space and time corrected H i j was computed as (the time correction accounts for the b, l, τ effects and is for EN data only) are the mean, VDR and vertical acceleration, and e and f account for the seasonal variation of height. Again, we used the robust least-squares estimator to determine the parameters in Eq. (5). Note that the terms  h e f , , (1) are just the initial estimates for VDR and annual variations; the final VDR is  h and the amplitudes of the annual variations are e and f. In Sec. S3 of Supplementary Information (SI), we selected an optimal waveform retracker to recover precise ranges from TP, J1, J2 and EN from waveforms like Waveform C in Fig. 7 for the SJV and NCP. Note that backscatter effects are not considered for TP, J1 and J2 in this study because of a lack of the necessary data.
The repeat-track method also diminishes the anisotropy roughness effect due to non-circular radar polarization, which causes different radar echoes for descending and ascending passes for the EN altimeter system. In fact, such height measurements do not repeat the exact same locations (maximum lateral offset: 1 km) and the measurements are affected by a number of time-varying factors. As a summary, the following conditions should be met to obtain precise heights on land: (a) sufficiently flat terrain within the altimeter footprint that may be covered with short vegetation, such as crops, (b) low percentage of buildings in the footprint area, (c) use of a proper waveform retracker, and (d) proper accounting of radar backscatter effects using an empirical threshold value. Parts of the croplands in the SJV, CT and the NCP meet conditions (a,b). If the proportion of cropland within a "bin" of selected and processed data is greater than 50%, the accuracy of the resulting altimeter-derived VDR is generally better than 1 cm yr −1 (at 1-km resolution), based on the assessments using leveling data (CT , Table S3) and crossover differences (CT and NCP, Table S4).
Typically, there will be inter-mission altimeter range biases between two altimeter missions, and such biases depend on surface attributes (ocean, land, ice and river, etc). The combination of inter-mission altimeter bias and terrain-induced height difference between reference tracks will lead to an apparent discontinuity between the two time series of height changes from two satellite missions at the crossovers of two ground tracks (Fig. S6). Because the VDRs in this study are computed for individual satellites and are the time derivatives of heights, the inter-mission discontinuities will not affect individual rates. In this study, we compute the cumulative subsidence, S, at a given bin from multiple satellite missions, such as TP, J1 and J2, using