Land Subsidence in the Texas Coastal Bend: Locations, Rates, Triggers, and Consequences

: Land subsidence and sea level rise are well-known, ongoing problems that are negatively impacting the entire Texas coast. Although ground-based monitoring techniques using long-term global navigation satellite systems (GNSS) records provide accurate subsidence rates, they are labor intensive, expensive, time-consuming, and spatially limited. In this study, interferometric synthetic aperture radar (InSAR) data and techniques were used to map the locations and quantify rates of land subsidence in the Texas Coastal Bend region during the period from October 2016 to July 2019. InSAR-derived land subsidence rates were then validated and calibrated against GNSS-derived rates. The factors controlling the observed land subsidence rates and locations were investigated. The consequences of spatial variability in land subsidence rates in Coastal Bend were also examined. The results indicated that: (1) land subsidence rates in the Texas Coastal Bend exhibited spatial variability, (2) InSAR-derived land subsidence rates were consistent with GNSS-derived deformation rates, (3) land subsidence in the Texas Coastal Bend could be attributed mainly to hydrocarbon and groundwater extraction as well as vertical movements along growth faults, and (4) land subsidence increased both ﬂood frequency and severity in the Texas Coastal Bend. Our results provide valuable information regarding not only land deformation rates in the Texas Coastal Bend region, but also the effectiveness of interferometric techniques for other coastal rural areas around the globe.


Introduction
Currently, more than 30% of the world's population lives in coastal areas, and 50% are expected to do so by 2030 [1]. However, most coastal areas are witnessing land surface deformation [2][3][4]. Surface deformation, represented by subsidence and uplift, is a change in the elevation of the Earth's surface in response to localized geologic changes, such as sediment compaction, fluid withdrawal, and fault slipping. Subsidence is the "sinking" or decrease in the elevation of the Earth's surface in relation to surrounding areas, while uplift is a "rise" or increase in the elevation of the Earth's surface relative to the surrounding areas.
Land subsidence makes coastal communities more vulnerable to natural forces, such as flooding, hurricanes, and sea level rise [5][6][7]. Hurricanes and flooding are both associated with a loss of life, livestock, crops, and natural habitat; property damage; and contamination of surface and groundwater resources [8,9]. Sea level rise is usually associated with inundation of wetlands and deltas, enhanced coastal erosion, increased vulnerability of coastal environments to storms, and seawater intrusion/pollution in coastal aquifers and surface water [10,11].
Land subsidence and sea level rise are ongoing phenomena that negatively impact the state of Texas [12]. More than 24% of Texans live in coastal areas. Texas coastal communities are at an increased risk compared to other communities because of land subsidence [13]. Texas coastal regions are also extremely important to the economy of Texas and the entire United States; they provide almost 50% of the national gross domestic product [14].
In Texas, land subsidence is attributed mainly to sediment compaction, growth faulting, and fluid (groundwater and hydrocarbon) withdrawal [15][16][17]. Sediment compaction is a function of sediment types, grain size, porosity, and water content. Highly compressible deltaic sediments (e.g., clays) were reported to subside, erode, break up, and become submerged more rapidly compared to other sediment types (e.g., sands) [18]. Generally, sediment compaction can be classified as natural (e.g., natural overburden loading, glacial isostatic adjustment, downward tectonic movements) or anthropogenic (e.g., fluid withdrawal). The average natural subsidence rate along the 600-km Texas coastline was estimated at 1.4 mm/year from the 1970s to the 2010s [19]. Growth faults in coastal areas are known to cause subsidence because they are usually represented by normal faults with a hanging wall on the coastal side of the fault plane [20]. There is strong evidence that the extraction of fluids from subsurface sediments can accelerate or reinitiate movement along these faults [21]. In addition, fluid extraction is an important controlling factor in land subsidence. When a system is depleted, the effective stress on the pore skeleton changes, which leads to significant compaction in unconsolidated sediments, expressed as a surface elevation change [22,23]. Land deformation caused by fluid pumping varies based on the thickness of the material that is deforming, the magnitude of the change in hydraulic head, and the specific skeletal storage that relates to the head history of the compacting sediment [24]. Within an elastic porous medium, water level changes and land subsidence are expected to be linearly related by compressibility [25].
Although ground-based monitoring techniques, such as global navigation satellite systems (GNSS) stations, provide accurate land subsidence rates, these techniques are labor intensive, expensive, time-consuming, and spatially limited (in other words, they are point measurements) [17,26]. In addition, interruptions in GNSS data acquisition usually result in an exaggerated estimate of land subsidence rates [27]. The use of interferometric synthetic aperture radar (InSAR) provides a significant solution to these limitations. InSAR is a remote-sensing technique that uses multi-temporal radar images for the detection of elevation changes. It provides a better, more homogenous spatial coverage as well as a more complete synoptic view of land surface deformation rates.
InSAR generates interference patterns when there are differences in the phases of radar waves between two or more radar images acquired over the same location at different times. The resulting interferograms represent changes in distance between the ground surface and the radar instrument. These distance changes are used to generate maps of land surface deformation over the areas of interest by determining whether the target moved toward the satellite (uplift) or away from the satellite (subsidence) between acquisitions [27][28][29][30][31].
InSAR has been successfully used to map fault slip from earthquakes [32], mine subsidence [33], aquifer compaction from pumping [34,35], landslides [36], and land subsidence in urban areas [15,16]. Several studies have used InSAR data and techniques for land deformation mapping in areas surrounding Houston, Texas [12,16,37,38]. However, none of these studies have quantified the land deformation rates in the Coastal Bend region of south Texas (Figure 1).
In this study, InSAR data and techniques were used to map the locations, quantify the rates, investigate the controlling factors, and examine the consequences of land subsidence in the Coastal Bend region in Texas ( Figure 1). This research provides valuable information regarding not only land deformation rates in the Texas Coastal Bend region, but also the effectiveness of InSAR techniques for other coastal rural areas around the world that lack significant GNSS coverage and pose significant challenges to InSAR applicability.
Remote Sens. 2022, 14,192 3 of 24 information regarding not only land deformation rates in the Texas Coastal Bend region, but also the effectiveness of InSAR techniques for other coastal rural areas around the world that lack significant GNSS coverage and pose significant challenges to InSAR applicability.

Study Area
Much of the Coastal Bend region of Texas ( Figure 1) is comprised of rural countryside and small towns. Therefore, conventional methods of tracking land deformation, such as GNSS (green dots; Figure 1), do not provide a spatially adequate dataset over the entire region, because the GNSS stations are concentrated only in the larger towns. Acquiring additional field GNSS data depends mainly on resource availability and site access. Land cover in the study area consists primarily of agriculture and rangelands. The agriculture and urban areas have been steadily growing in area, while the forest, rangeland, and wetland portions are all steadily shrinking [39].
The Coastal Bend region is situated on a passive depositional margin with geologic formations consisting primarily of alternating sand, silt, and clay layers and ranging in age from the Miocene to Holocene [40] (Figure 2a). The sedimentary sections are estimated to be 4,500-5,500 m thick near the Gulf of Mexico coastline. However, these sedimentary

Study Area
Much of the Coastal Bend region of Texas ( Figure 1) is comprised of rural countryside and small towns. Therefore, conventional methods of tracking land deformation, such as GNSS (green dots; Figure 1), do not provide a spatially adequate dataset over the entire region, because the GNSS stations are concentrated only in the larger towns. Acquiring additional field GNSS data depends mainly on resource availability and site access. Land cover in the study area consists primarily of agriculture and rangelands. The agriculture and urban areas have been steadily growing in area, while the forest, rangeland, and wetland portions are all steadily shrinking [39].
The Coastal Bend region is situated on a passive depositional margin with geologic formations consisting primarily of alternating sand, silt, and clay layers and ranging in age from the Miocene to Holocene [40] (Figure 2a). The sedimentary sections are estimated to be 4500-5500 m thick near the Gulf of Mexico coastline. However, these sedimentary sections are heavily faulted. The faulting tends to be regional, and varies in complexity from system to system [41] (Figure 2a,b). Growth faults, parallel to the coastline, are common within the Texas coast and are attributed to the loading of the unconsolidated sediments in this region [42,43]. These faults are part of the Cenozoic Wilcox and Frio fault trends, which are composed of large displacements, dominantly down to the basin growth faults [44].
sections are heavily faulted. The faulting tends to be regional, and varies in complexity from system to system [41] (Figures 2a and 2b). Growth faults, parallel to the coastline are common within the Texas coast and are attributed to the loading of the unconsolidated sediments in this region [42,43]. These faults are part of the Cenozoic Wilcox and Frio faul trends, which are composed of large displacements, dominantly down to the basin growth faults [44]. (a) Surface lithology within the Texas Coastal Bend region; also shown is the spatial distribution of growth faults (orange lines). (b) Growth fault density (in km/km 2 ) within the study area. (c) Spatial distribution of sand percentage, extracted from geophysical logs (black dots), within the GCA system. Circles 1-7 are areas that were identified to have witnessed a significant land subsidence during the investigated period (October 2016-July 2019).
The hydrogeology of the study area is dominated by the Gulf Coast aquifer (GCA system, which runs parallel to the Gulf of Mexico coastline from the Louisiana border in the northeast to the border with Mexico in the southwest (Figure 1, inset). The GCA sys tem consists of three main aquifer units: the Jasper, Evangeline, and Chicot aquifers. These aquifers thicken and dip toward the Gulf of Mexico [45]. Groundwater is typically uncon fined or semi-confined in this system. The specific yields from GCA are quite low The hydrogeology of the study area is dominated by the Gulf Coast aquifer (GCA) system, which runs parallel to the Gulf of Mexico coastline from the Louisiana border in the northeast to the border with Mexico in the southwest (Figure 1, inset). The GCA system consists of three main aquifer units: the Jasper, Evangeline, and Chicot aquifers. These aquifers thicken and dip toward the Gulf of Mexico [45]. Groundwater is typically Remote Sens. 2022, 14, 192 5 of 24 unconfined or semi-confined in this system. The specific yields from GCA are quite low compared to other unconfined sedimentary aquifer systems, because interbedded silt/clay lenses confine the aquifers in smaller areas [46]. The large volume of water pumped from the aquifers has caused subsidence, particularly in the Houston area [12,16]. Groundwater extraction within the GCA is estimated to be 518 × 10 9 m 3 per year. Most of the pumping in the area is used for irrigation [46]. Figure 3a shows the spatial distribution of the average annual groundwater extraction rates in each county within the Coastal Bend region between January 2000 and January 2015, as retrieved from the Texas Water Development Board (TWDB) website (https://www.twdb.texas.gov accessed on 1 January 2020). Higher groundwater extraction activities were observed in Victoria, Hidalgo, and Atascosa counties ( Figure 3a). compared to other unconfined sedimentary aquifer systems, because interbedded silt/clay lenses confine the aquifers in smaller areas [46]. The large volume of water pumped from the aquifers has caused subsidence, particularly in the Houston area [12,16]. Groundwater extraction within the GCA is estimated to be 518 × 10 9 m 3 per year. Most of the pumping in the area is used for irrigation [46]. Figure 3a shows the spatial distribution of the average annual groundwater extraction rates in each county within the Coastal Bend region between January 2000 and January 2015, as retrieved from the Texas Water Development Board (TWDB) website (https://www.twdb.texas.gov accessed on 1 January 2020). Higher groundwater extraction activities were observed in Victoria, Hidalgo, and Atascosa counties ( Figure 3a). Oil and gas extraction practices are significant within the study area. High oil and gas extraction rates, along with the fluid injection operations that accompany these activities, have been shown to cause localized land deformation within oil and gas fields [47]. Figure 4a shows the average annual oil/gas extraction rates from the Texas Railroad Commission (TRC) website (https://www.rrc.state.tx.us/ accessed on 1 January 2020) for each county between January 2000 and January 2015. Hidalgo, Cameron, and Karnes counties experienced the most hydrocarbon extraction activities (Figure 4a). Oil and gas extraction practices are significant within the study area. High oil and gas extraction rates, along with the fluid injection operations that accompany these activities, have been shown to cause localized land deformation within oil and gas fields [47]. Figure 4a shows the average annual oil/gas extraction rates from the Texas Railroad Commission (TRC) website (https://www.rrc.state.tx.us/ accessed on 1 January 2020) for each county between January 2000 and January 2015. Hidalgo, Cameron, and Karnes counties experienced the most hydrocarbon extraction activities (Figure 4a).

Data
We used synthetic aperture radar (SAR) data to map areas experiencing land subsidence and quantify the displacement rates. The SAR data was acquired by the Sentinel-1 mission (C-band with a 5.6 cm wavelength) operated by the European Space Agency (ESA) under the Copernicus program. Sentinel-1 was launched in April 2014 to provide SAR imaging for all global landmasses and coastal zones with high spatial (5 m range × 20 m azimuth) and temporal (6-12 days) resolutions [48]. The Sentinel-1 data used in this study were downloaded from the Alaskan Satellite Facility Distributed Active Archive Center (https://asf.alaska.edu/ accessed on 1 August 2019). A total of 68 single look complex Sentinel-1 scenes spanning from October 2016 to July 2019 were used in this study (Tables A1 and A2). The Sentinel-1 scenes were acquired in the interferometric wide imaging mode (scene width: 250 km) along track 107, frames 83 and 88, in the ascending direction (Table 1; Figure 1). Note that it is preferable to combine both the descending and ascending datasets to extract three-dimensional (3D) displacement rates [49]. However, the descending direction did not provide adequate coverage of our study area.

Data
We used synthetic aperture radar (SAR) data to map areas experiencing land subsidence and quantify the displacement rates. The SAR data was acquired by the Sentinel-1 mission (C-band with a 5.6 cm wavelength) operated by the European Space Agency (ESA) under the Copernicus program. Sentinel-1 was launched in April 2014 to provide SAR imaging for all global landmasses and coastal zones with high spatial (5 m range × 20 m azimuth) and temporal (6-12 days) resolutions [48]. The Sentinel-1 data used in this study were downloaded from the Alaskan Satellite Facility Distributed Active Archive Center (https://asf.alaska.edu/ accessed on 1 August 2019). A total of 68 single look complex Sentinel-1 scenes spanning from October 2016 to July 2019 were used in this study (Tables A1 and A2; Appendix A). The Sentinel-1 scenes were acquired in the interferometric wide imaging mode (scene width: 250 km) along track 107, frames 83 and 88, in the ascending direction (Table 1; Figure 1). Note that it is preferable to combine both the descending and ascending datasets to extract three-dimensional (3D) displacement rates [49]. However, the descending direction did not provide adequate coverage of our study area. To validate and calibrate the results of the land deformation estimates derived from the InSAR technique, we used GNSS displacement records retrieved over the study area (Table A4; Appendix A). The datasets were acquired from the Nevada Geodetic Laboratory (NGL) (http://geodesy.unr.edu accessed on 1 June 2020) [50]. The NGL processes GNSS data and provides daily precise point positioning (PPP) solutions as well as velocity products [50]. They employ the Median Interannual Difference Adjusted for Skewness (MIDAS) method in calculating the site velocities and uncertainties [51]. This method employs the median of slopes of a great number of selected displacements 1 year apart to represent the overall site velocity. The uncertainties in the generated velocities were represented by the standard deviation of these slopes. The reference frame for the NGL GNSS data was the International GNSS Service Reference Frame 2014 (IGS14) [52]. The entire time span (average: approximately 10 years; range: approximately 3-24 years) for 13 GNSS sites (Figure 1), spatially distributed throughout the study area, was used since longer records provide more accurate velocity estimates [53].
The generated land subsidence rates were corresponded with the areas that experienced flooding following two major flood events, Hurricane Harvey (August 2017) and Hurricane Hanna (July 2020), in order to observe potential codependent relationships between the two processes. SAR data have also been extensively used for flood mapping under widely varying weather conditions [38,54]. The flood mapping approach relies on using multiple pre-event reference SAR images that do not exhibit flooding in contrast with an image acquired after the flood event. The images acquired before the flood event provide a baseline backscatter intensity value for the study area, while the post-event image is selected to compare against this baseline. Flooded areas show a significantly lower backscatter response as water attenuates electromagnetic waves. The ratio between the pre-and post-event images can then be used to classify as flooded or unflooded areas [55,56]. For flood mapping in this study, we used the ascending track 107 frame 88 (Table 2) for the Hurricane Harvey flood event. For Hurricane Hanna, we used the descending track 41 frame 503, because it provided better temporal coverage of the storm event (Tables 2 and A3; Appendix A). Several processes that directly or indirectly control land deformation were examined in a geographic information system (GIS) environment to investigate the factors controlling the observed land deformation rates. These include geologic structures (e.g., growth faults), subsurface lithology, and fluid (groundwater and hydrocarbon) extraction rates.
Growth fault data were digitized from earlier studies, focusing on the Wilcox and Frio growth fault trends, published by the Bureau of Economic Geology at the University of Texas at Austin [57]. In addition, localized growth faults in the Corpus Christi and Nueces Bay areas were also digitized from Hammes et al. [58] (Figure 2). ArcGIS tools were used to calculate fault density, shown in Figure 2b. In this study, the growth faults were classified as low (0-0.1 km/km 2 ), medium (0.1-0.15 km/km 2 ), and high (>0.15 km/km 2 ) density ( Figure 2b).
A continuous distribution of the sand percentage within the GCA members was constructed from point measurements taken at geophysical logs ( Figure 2c) and published by the TWDB [59]. Only logs with 70% coverage across the entire aquifer thickness were used in these calculations. The sand percentage for each geologic unit was calculated by summing the total sand amount across the thickness of the geologic unit and dividing by the total thickness for which the lithology was characterized. These calculations were normalized to the known lithology; for example, if the unit was 100 m thick, the lithology was only known for 85 m, and the sand thickness was 75 m, then the value would be 88% rather than 75% [60]. The spatial distribution of the sand percentage within the study areas is shown in Figure 2c.

Methods
This research involved three major steps to address the objectives of the study outlined above. The first step was to map the spatial distribution of land subsidence occurrences across the study area and quantify their rates. Second, relevant datasets (e.g., fluid extraction rates, geologic structures, subsurface lithology) that directly or indirectly control land subsidence processes were analyzed to investigate the factors that drive the observed subsidence rates and locations. Finally, the generated land subsidence rates were correlated, in a GIS environment, with areas that experienced flooding following major flood events induced by severe weather (Hurricane Harvey and Hurricane Hanna) to examine the significance and implication of land subsidence processes in enhancing the spatial distribution of floodwater and potentially altering the gradient of flow.

Mapping Land Subsidence Rates and Locations
Persistent scatterer interferometry (PSI) was used in this study to quantify land subsidence rates in the Texas Coastal Bend region. PSI uses a series of SAR images to detect coherent persistent/permanent scatterers (PS) points in the region of interest. A PS is a pixel or point that remains stable throughout the entire period of the SAR acquisition. A primary image is chosen from the group of images based on favorable geometry, high coherence, and minimal atmospheric interference. Generally, urban areas are best suited for PS points because many fixed reflectors (e.g., buildings and utility poles) are available. Each image in the rest of the SAR image stack, called a secondary image, is geometrically aligned with the primary image. Interferograms are then generated from each pair.
The main goal of PSI analysis is to decipher the coherent radar signals from incoherent noise by using the PS points as a reference [61]. PS candidates are initially selected based on a high amplitude dispersion threshold (0.4); however, temporal coherence is the main parameter that is used to retain PS that exhibit phase stability [62]. Since our study area comprises significant non-urban areas, the method selected in this study [62] ensures the selection of PS points outside the urban areas and their deformation rates are calculated, providing a comprehensive overview of the regional deformation process.
The interferometric phase observation for each point is generated according to this equation: where ϕ int is the interferometric phase, ϕ topo is the topographic phase, ϕ de f o is the deformation phase, ϕ orb is the deterministic flat earth phase and the residual phase signal due to orbit in determination, ϕ atm is the atmospheric phase, ϕ scat is the phase due to a temporal and spatial change in the scatter characteristics of the earth surface between the observation times, and ϕ noise is the phase degradation factors caused by thermal noise, coregistration noise, and interpolation noise [63,64].
InSAR uses the phase difference between different acquisitions to determine the displacement of each PS in the line of sight (LOS). By identifying and removing other terms, the displacement is calculated using this equation: where ∆ϕ represents the change in phase (the difference in ϕ int related to different acquisition times or ϕ int2 − ϕ int1 ), λ represents wavelength (Sentinel-1: 5.6 cm), and d r represents the relative displacement [63]. Assuming negligible horizontal motion, the vertical motion component was then extracted from the InSAR-derived LOS measurements using the following equation [64]: where d v represents vertical motion and θ is the incidence angle. The incidence angle of the pixels in our study area ranged from 30.06 • to 45.24 • . InSAR processing was completed using the open-source software packages provided by the ESA, the Sentinel Application Platform (SNAP) [65], and the Stanford method for persistent scatterers (StaMPS; [28]). The PSI processing was conducted in two separate workflows: single primary differential interferometric synthetic aperture radar (DInSAR) processing using SNAP (Figure 5a), and PSI analysis using StaMPS (Figure 5b).
The primary scene was selected using the stack overview tool in SNAP. The primary image was then split into a sub-swath using the S-1 tops split function to make it compatible as a snap2stamps input. The snap2stamps package provides a set of Python scripts that call tools from SNAP to perform automatic interferogram stacking that is compatible with StaMPS PSI analysis [65]. Snap2stamps is implemented by running the following steps ( Figure 5a):

1.
Secondary image preparation: The available SAR data within the project folder is separated into individual folders by acquisition date.

2.
Secondary splitting: The secondary images are split into the correct area of interest and orbital corrections are applied. 3.
Co-registration and interferogram generation: The SAR data is co-registered, producing a set of interferograms with the topographic phase removed, and separate files are prepared for the execution of StaMPS analysis. 4.
StaMPS export: The interferograms are converted into binary files that are compatible as StaMPS inputs.
We then ran the StaMPS process chain from steps 1-7 as described in the StaMPS user manual [66] to generate land subsidence rates (Figure 6a). These steps included noise calculation, PS selection and weeding, phase correction and unwrapping, look angle error estimation, and atmospheric filtering (Figure 5b). We used a GNSS station (txae; Figure 1) with stable (0 mm/yr; Table A4; Appendix A) deformation rate as a reference point (in StaMPS processing) against which all PS points in the study area were referenced and compared. For the removal of the atmospheric phase screen (APS), we employed the Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) and the Generic Atmospheric Correction Online Service for InSAR (GACOS) methods for corrections [67][68][69]. Given the size of our study area, it is possible that there were remnants of the atmospheric effects even after applying the TRAIN and GACOS corrections. We quantified this effect as the difference between the InSAR-derived deformation rates and the closest GNSS-derived deformation rates for each of the seven investigated regions (polygons 1-7; Figure 6a). and orbital corrections are applied.
3. Co-registration and interferogram generation: The SAR data is co-registered, producing a set of interferograms with the topographic phase removed, and separate files are prepared for the execution of StaMPS analysis.
4. StaMPS export: The interferograms are converted into binary files that are compatible as StaMPS inputs.  We then ran the StaMPS process chain from steps 1-7 as described in the StaMPS user manual [67] to generate land subsidence rates (Figure 6a). These steps included noise calculation, PS selection and weeding, phase correction and unwrapping, look angle error estimation, and atmospheric filtering (Figure 5b). We used a GNSS station (txae; Figure 1) with stable (0 mm/yr; Table A4) deformation rate as a reference point (in StaMPS processing) against which all PS points in the study area were referenced and compared. For the removal of the atmospheric phase screen (APS), we employed the Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) and the Generic Atmospheric Correction Online Service for InSAR (GACOS) methods for corrections [68][69][70]. Given the size of our study area, it is possible that there were remnants of the atmospheric effects even after applying the TRAIN and GACOS corrections. We quantified this effect as the difference between the InSAR-derived deformation rates and the closest GNSS-derived deformation rates for each of the seven investigated regions (polygons 1-7; Figure 6a).

InSAR-GNSS Validation
The linear vertical velocities were estimated at the 13 GNSS locations using daily PPP data provided by NGL (Table A4). A least-square regression was used to fit a trend to the daily GNSS data. These GNSS-derived velocities were then compared to the InSAR-derived velocities. Correlation coefficients (R) were calculated for the two datasets (Figure

InSAR-GNSS Validation
The linear vertical velocities were estimated at the 13 GNSS locations using daily PPP data provided by NGL (Table A4; Appendix A). A least-square regression was used to fit a trend to the daily GNSS data. These GNSS-derived velocities were then compared to the InSAR-derived velocities. Correlation coefficients (R) were calculated for the two datasets (Figure 6b). We used the average velocity of InSAR points within a 500 m radius of the GNSS stations. The uncertainties in GNSS-and InSAR-derived velocities (Figure 6b) were calculated as follows: (1) horizontal error bars (x-axis; InSAR) represent the standard deviation of all measurements located within a circle of 500 m radius that is centered at the location of the GNSS station; (2) vertical error bars (y-axis; GNSS) represent the 95% confidence interval (CI) calculated using the empirical equation provided by Wang [53]: where T is the year range of the time series in years. Table A4 (Appendix A) lists T, rates and CIs for all GNSS stations used in this study. It is worth mentioning that the horizontal error bar was overestimated, since Sentinel-1 can detect spatial variations in deformation rates within a 500 m radius; however, the 500 m was found to be the minimum radius to generate a statistically significant number of deformation points that could be used to calculate the standard deviation.

Mapping Flooded Areas
Senitinel-1 images were used for flood mapping in our study area. We mapped the flooded areas using the SARscape Analytics package within the environment for visualizing images (ENVI) software (Figure 5c). The SAR flood mapping routine was used [70]. To map inundated areas after Hurricane Harvey, we used four pre-event images as a control background for the unflooded area. We then picked a post-event image that was as close as possible to the major event. Using the SAR flood mapping tool, we then picked inundated areas (resolution: 25 m) by selecting a ratio between pre-and post-event coherence to determine which areas were flooded by the storm events (Figure 7). The generated flood maps were validated, by visual inspection, against the available visible images and reported rainfall rates during each event.
where T is the year range of the time series in years. Table A4 lists T, rates and CIs for all GNSS stations used in this study. It is worth mentioning that the horizontal error bar was overestimated, since Sentinel-1 can detect spatial variations in deformation rates within a 500 m radius; however, the 500 m was found to be the minimum radius to generate a statistically significant number of deformation points that could be used to calculate the standard deviation.

Mapping Flooded Areas
Senitinel-1 images were used for flood mapping in our study area. We mapped the flooded areas using the SARscape Analytics package within the environment for visualizing images (ENVI) software (Figure 5c). The SAR flood mapping routine was used [71]. To map inundated areas after Hurricane Harvey, we used four pre-event images as a control background for the unflooded area. We then picked a post-event image that was as close as possible to the major event. Using the SAR flood mapping tool, we then picked inundated areas (resolution: 25 m) by selecting a ratio between pre-and post-event coherence to determine which areas were flooded by the storm events (Figure 7). The generated flood maps were validated, by visual inspection, against the available visible images and reported rainfall rates during each event. .

Land Subsidence Rates in the Texas Coastal Bend
The InSAR-derived vertical ground deformation rates extracted from Sentinel-1 data during the period from October 2016 to July 2019 over the Texas Coastal Bend region are shown in Figure 6a. The deformation in the study area was spatially heterogenous with localized areas of subsidence as high as −27 mm/yr. In addition, some localized areas showed a land uplift greater than 10 mm/yr. Figure 6c shows a histogram of land deformation rates over the Texas Coastal Bend region. This histogram indicates that the average InSAR-derived deformation rate was estimated at −1.59 ± 4.07 mm/year (average ± standard deviation).
The InSAR rates were then compared to those extracted from 13 GNSS stations located within the study area. Figure 6b shows a scatterplot of InSAR-derived versus GNSS-derived ground deformation rates along with a 1:1 line. Note that a complete 1:1 correspondence between InSAR-and GNSS-derived deformation rates was not expected, since we averaged InSAR points within a 500 m radius of the GNSS stations. Figure 6b shows that some of the GNSS rates were consistent with the deformation rates obtained from InSAR analysis. The R between the two datasets was estimated to be 0.24. The mean difference between the GNSS and InSAR-derived velocities at all 13 stations was estimated to be −0.36 ± 3.28 mm/yr (average ± standard deviation). The difference between the InSAR-and GNSS-derived deformation rates generally increased with distance from the reference GNSS station (txae). For example, a difference of −0.14 mm/yr was found at the txfe station (54 km from txae), whereas a difference of −6.04 mm/yr was found at the txvc station (162 km from txae). This is consistent with other findings from similar studies [71]. Despite the fact that the InSAR deformation rates were located within the order of variability of the GNSS-derived rates, there are several potential factors that contribute to the inconsistency between the two deformation rates. These include differences in sampling rate, uncompensated atmospheric artefacts, orbital errors, unwrapping errors over low coherence areas, and/or inherent errors associated with GNSS observations [72][73][74].
Seven regions were observed to experience significant (>3.0 mm/yr) land subsidence rates along the Texas Coastal Bend region during the investigated period: Victoria, George West, Refugio, Falfurrias, Karnes City, McAllen, and Nueces Bay/Corpus Christi (Table 3; polygons 1-7; Figure 6a). Each of these regions contained a minimum of 200 PS points; they extended outward in a circular shape, from the highest subsidence value at each location (polygons 1-7; Figure 6a). The reported land subsidence rate for each region (Figure 6a) was calculated as the average of all InSAR rates within that region ( Table 3). Histograms of InSAR-derived land subsidence rates over these regions are shown in Figure A1; Appendix A. Errors in land subsidence rates in each region (σ InSAR ) were calculated by using this equation: where SD InSAR is the standard deviation of all InSAR rates within each region (Table 3), and ∆ is the deviation of the InSAR rates from the closest GNSS station rate (txvc for Victoria; txgw for George West; txgo for Refugio; txff for Falfurrias; txbe for Karnes City; txpr for McAllen; and txcc for Nueces; Table A4; Appendix A). The InSAR uncertainties generated from this approach (Equation (5)) were higher compared to that (~1.5 mm/yr) derived from Equation (4) using almost 3 years of InSAR data (October 2016-July 2019). Surprisingly, no significant correlation was observed between the spatial distribution of these seven regions and the sand percentage (Figure 2c), since a wide range (20-100%) was observed. The values for the majority of these regions, other than George West (polygon 2) and Falfurrias (polygon 4), were found to be within 50-70% (Figures 2c and 6a). George West experienced a higher average percentage (60-80%), whereas Falfurrias experienced a lower (20-40%) percentage. This might be true only for Polygon 4 (e.g., lower sand percentage and higher subsidence); however, we do not expect subsided areas to have sand as subsurface lithology. Land subsidence in the Texas Coastal Bend is believed to be attributed mainly to fluid extraction and vertical movements along growth faults (Figure 8). Land subsidence rates in polygons 1 (Victoria) and 2 (George West) were attributed mainly to oil/gas extraction, in polygons 5 (Karnes) and 6 (McAllen) to groundwater extraction, and in polygons 3 (Refugio), 4 (Falfurrias), and 7 (Nueces) to high fault density. Below, we explain in detail the factors that control the observed land subsidence rates in each of the seven (Figure 6a) regions in the Texas Coastal Bend. 33 Surprisingly, no significant correlation was observed between the spatial distribution of these seven regions and the sand percentage (Figure 2c), since a wide range (20-100% was observed. The values for the majority of these regions, other than George West (pol ygon 2) and Falfurrias (polygon 4), were found to be within 50-70% (Figures 2c and 6a) George West experienced a higher average percentage (60-80%), whereas Falfurrias expe rienced a lower (20-40%) percentage. This might be true only for Polygon 4 (e.g., lower sand percentage and higher subsidence); however, we do not expect subsided areas to have sand as subsurface lithology.
Land subsidence in the Texas Coastal Bend is believed to be attributed mainly to fluid extraction and vertical movements along growth faults (Figure 8). Land subsidence rates in polygons 1 (Victoria) and 2 (George West) were attributed mainly to oil/gas extraction in polygons 5 (Karnes) and 6 (McAllen) to groundwater extraction, and in polygons 3 (Refugio), 4 (Falfurrias), and 7 (Nueces) to high fault density. Below, we explain in detai the factors that control the observed land subsidence rates in each of the seven (Figure 6a regions in the Texas Coastal Bend.  Figure 6a). Fluid (groundwater and hydrocarbon) extraction is reported as the difference between the extraction rates during the study period (2016-2019) and the historical period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). Land subsidence rates in polygons 1 and 2 were attributed mainly to oil/gas extraction, in polygons 5 and 6 to groundwater extraction, and in polygons 3, 4, and 7 to high fault density. Generally, land subsidence increased with fluid extraction and fault density.

Factors Controlling Observed Subsidence Rates in the Texas Coastal Bend
The highest subsidence rates were observed in the city of Victoria, with an average rate of -7.55 ± 5.51 mm/yr (polygon 1; Figure 6a). Localized subsidence rates as high as 15 mm/yr were also observed south of Victoria. Increased subsidence rates could be at  Figure 6a). Fluid (groundwater and hydrocarbon) extraction is reported as the difference between the extraction rates during the study period (2016-2019) and the historical period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). Land subsidence rates in polygons 1 and 2 were attributed mainly to oil/gas extraction, in polygons 5 and 6 to groundwater extraction, and in polygons 3, 4, and 7 to high fault density. Generally, land subsidence increased with fluid extraction and fault density.

Factors Controlling Observed Subsidence Rates in the Texas Coastal Bend
The highest subsidence rates were observed in the city of Victoria, with an average rate of -7.55 ± 5.51 mm/yr (polygon 1; Figure 6a). Localized subsidence rates as high as -15 mm/yr were also observed south of Victoria. Increased subsidence rates could be attributed to enhanced oil/gas extraction activities, as well as vertical movements along growth faults. This assumption is supported by the observed increase in oil/gas extraction rates (431 × 10 3 BBL/yr; Figure 4b) during the investigated period compared to the long-term average annual extraction rates (oil/gas: 664 × 10 3 BBL/yr; Figure 4a). In addition, this area has a high fault density >0.15 km/km 2 (Figure 2b), which could also accelerate the subsidence rates. Victoria witnessed a decline in groundwater extraction rates (−1.7 × 10 6 m 3 /yr; Figure 3b) during the investigated period compared to the historical rates (18.3 × 10 6 m 3 /yr [high]; Figure 3a).
Live Oak County near the town of George West also experienced land subsidence rates of −9.02 ± 5.38 mm/yr (polygon 2, Figure 6a). We believe the subsidence in this area could be attributed to a combination of both growth faulting and enhanced oil/gas extraction activities. Live Oak County reported a very large increase in oil/gas extraction (7.65 × 10 6 BBL/yr; Figure 4b) compared to its high historical rate (2.98 × 10 6 BBL/yr; Figure 4a), while there was a slight decrease in groundwater extraction (−78.2 × 10 3 m 3 /yr; Figure 3b) compared to the medium historical average of 5.2 × 10 6 m 3 /yr (Figure 3a). In addition, there was a high density of growth faults observed in the area surrounding George West (>0.15 km/km 2 , Figure 2b) where the increased subsidence rates were observed.
The town of Refugio witnessed land subsidence with an average rate of −8.67 ± 5.06 mm/yr (polygon 3, Figure 6a). We believe vertical movement along growth faults is the main contributing factor for the observed subsidence in this area, because the extraction rates for both oil/gas (−1.25 × 10 6 BBL/yr; Figure 4b) and groundwater (−251 × 10 3 m 3 /yr Figure 3b) decreased compared to their respective historic annual averages (oil/gas: 3.84 × 10 6 BBL/yr (high), Figure 4a; groundwater: 2.84 × 10 6 m 3 /yr (medium), Figure 3a). However, Refugio County had a high fault density (>0.15 km/km 2 ) overlapping the area of increased subsidence ( Figure 2b). Movement along these faults might be activated due to the ongoing fluid extraction activities.
The city of Falfurrias in Jim Wells County (polygon 4, Figure 6a) experienced a land subsidence rate of −9.48 ± 3.99 mm/yr. We believe this area is primarily controlled by growth faulting; the extraction rates in the area of both oil/gas and groundwater decreased (oil/gas: −4.8 × 10 3 BBL/yr: Figure 4b; groundwater: −2.0 × 10 6 m 3 /yr; Figure 3b) in comparison to the historical averages (oil/gas: 138 × 10 3 BBL/yr (low), Figure 4a; groundwater: 6.05 × 10 6 m 3 /yr (low), Figure 3a). Although this area had a low to medium amount of fault density (0-0.15 km/km 2 ; Figure 2b), there were two prominent growth faults in the area that sandwich the subsided region, which could potentially explain the observed subsidence.
Karnes City experienced an averaged land subsidence rate of -5.24 ± 4.67 mm/yr (polygon 5, Figure 6a). Fluid extraction is the primary controlling factor in the subsidence; Karnes County reported a large increase in both oil/gas (78.04 × 10 6 BBL/yr; Figure 4b) and groundwater (9.29 × 10 6 m 3 /yr; Figure 3b) extraction rates during the study period as compared to the historical rates (oil/gas: 18.23 × 10 6 BBL/yr (high), Figure 4a; 6.85 × 10 6 m 3 /yr (medium), Figure 3a). The fault density in the area was 0 km/km 2 ( Figure 2b); thus, we do not believe that growth faults were a controlling factor in this area.
Hidalgo County experienced land subsidence (−2.64 ± 2.92 mm/yr) centered around the McAllen metro area, which also includes the cities of Edinburgh and Mission. Localized areas experienced land subsidence rates as high as −13 mm/yr (polygon 6, Figure 6a). Land subsidence in this area is believed to be due to groundwater extraction. Hidalgo County reported an increase in groundwater extraction of 5.27 × 10 6 m 3 /yr during the study period, in conjunction with the high historical average annual rate (14.62 × 10 6 m 3 /yr; Figure 3a). Hidalgo County had a lower fault density value (<0.1 km/km 2 ). Hidalgo County also witnessed a slight decrease in oil/gas extraction during the study period (−6.30 × 10 3 BBL/yr), and historically has had a low amount of hydrocarbon extraction (59.39 × 10 3 BBL/yr).
The Nueces Bay area, especially in the port corridor, experienced land subsidence rates of −3.53 ± 1.92 mm/yr (polygon 7, Figure 6a). This is believed to be due to the presence of growth faulting in the area (fault density: >0.15 km/km 2 ; Figure 2b), which affects both Corpus Christi and Nueces Bays. We do not think fluid extraction is contributing to subsidence here, as Nueces and San Patricio Counties experienced a decrease in oil/gas extraction rates (−62.5 × 10 3 BBL/yr in San Patricio County; −262.5 × 10 3 BBL/yr in Nueces County; Figure 4b) compared to the historical annual rates (452 × 10 3 BBL/yr in San Patricio County; 501 × 10 3 BBL/yr in Nueces County; Figure 4a). San Patricio County even experienced a large decrease in groundwater extraction (−4.5 × 10 6 m 3 /yr Figure 3b), compared to high historical rates of 13.3 × 10 6 m 3 /yr (Figure 3a).

Consequences of Land Subsidence in the Texas Coastal Bend
The InSAR-derived vertical deformation rates were spatially correlated, in a GIS environment, with the spatial distribution of areas affected by recent flood events in the Texas Coastal Bend (Figure 7). We decided to examine areas that were inundated after Hurricane Harvey in August 2017 and Hurricane Hanna in July 2020. These two storms heavily impacted the South Texas coast with heavy rainfall rates [75,76]. Figure 7 shows the rainfall amounts at each rain gauge for each event. Hurricane Harvey was associated with rainfall rates ranging between 93 mm in Orange Grove and 330 mm in Refugio. Hurricane Hanna, however, was associated with rainfall rates between 165 mm in Weslaco and 311 mm north of McAllen. Mapped flooded areas after Harvey and Hanna were estimated at 1,020 and 2,225 km 2 , respectively. The majority of these areas are cultivated cropland. This is expected, because urban areas have better drainage systems than cultivated land. We investigated whether there was a connection between the observed subsidence rates and flooding in the aftermath of storm events within the study area.
We focused on two main areas, Hidalgo County (polygon 6, Figure 7) and Nueces Bay (polygon 7, Figure 7). These two areas were selected for two reasons: (1) they were most affected by flooding events after Hurricane Harvey (polygon 6, Figure 7) and Hurricane Hanna (polygon 7, Figure 7), and (2) they contain the two largest population centers within the Coastal Bend, McAllen (polygon 6) and Corpus Christi (polygon 7).
We compared the average subsidence rates of the non-flooded areas to the average subsidence rates of the flooded areas in each polygon. The subsidence rates within flooded areas were significantly higher than in the non-flooded areas. For example, the average subsidence rate of the non-flooded area in polygon 6 was −2.61 ± 1.91 mm/yr, while the average subsidence rate of the areas flooded in the aftermath of Hurricane Hanna within polygon 6 was −3.25 ± 1.72 mm/yr; the subsidence of non-flooded area was 80% of that of flooded areas. The same applies to polygon 7: the average subsidence rate of the non-flooded areas was estimated at −3.20 ± 2.14 mm/yr, which is only 62% of the average subsidence rate (−5.16 ± 2.42 mm/yr) of the areas flooded in the aftermath of Hurricane Harvey within polygon 7.
We examined the statistical significance of the difference in subsidence rates between flooded and non-flooded areas. A t-test was used to explore whether there was indeed a significant difference between the mean subsidence rates in flooded and non-flooded areas in polygons 6 and 7. For polygon 6, our results showed that the inundated areas in the aftermath of Hurricane Hanna experienced a significantly higher rate of subsidence than the non-flooded areas (t = 7.07, p = 3.14 × 10 −12 ). Similarly, for polygon 7, the areas flooded during Hurricane Harvey experienced a significantly higher rate of subsidence than the non-flooded areas (t = 7.36, p = 2.59 × 10 −11 ). The increased subsidence in these areas is believed to be a contributing factor in the areas becoming severely inundated following heavy rainfall events.
Although we saw a significant connection between flooded and subsided areas within polygon 6 and polygon 7, some regions within the Coastal Bend experienced land subsidence but not floods (e.g., polygons 1, 3, and 4). This could be due to the fact that the Coastal Bend is very diverse in the spatial domain and additional factors could exert controls on areas that are being flooded in the aftermath of heavy rainfall. Other regions (polygons 2 and 5) were away from the storm path. We also saw flooding in areas that were not experiencing severe subsidence in our study area. While we believe that subsidence is one component in causing susceptibility to floods, we also recognize that it is possible for these areas to experience flooding due to other major factors that affect flood susceptibility, such as the actual elevation of the flood plains, local topography, vegetation, the path of the storm, spatial variability in rainfall intensity, and hydrologic conditions preceding the storm event.

Discussion
This research integrated different remote-sensing data and techniques to develop a new land subsidence dataset for the south Texas coast. This study explained how natural (e.g., growth faults) and anthropogenic (e.g., fluid extraction) processes control coastal subsidence. Our study showed significant correlations between growth faulting, fluid extraction, and land subsidence in the Texas Coastal Bend. The average correlation coefficient between the InSAR-derived deformation rates and fluid extraction in Victoria, George West, Karnes City, and McAllen (polygons 1, 2, 5, and 6; Figure 6a) was found to be 0.80. A relatively lower correlation coefficient (0.42) was found between land subsidence and growth fault density in Victoria, George West, Refugio, Falfurrias, and Nueces Bay/Corpus Christi (polygons 1, 2, 3, 4, and 7; Figure 6a). This is congruent with previous InSAR studies that have been conducted in the Houston area and in other parts of the United States; these studies showed that growth faulting and fluid extraction were correlated with observed subsidence rates. For example, Qu et al. [16] reported a land subsidence rate of up to 53 mm/yr in northwest Houston, Texas between the 1990s and 2000s that was mainly related to fluid withdrawal. Bawden et al. [12] found that the northwestern portions of Harris County, Texas experienced rapid subsidence (rate: 15−30 mm/yr) due to groundwater extraction. Miller and Shirzaei [38] reported that between 2007 and 2011, northwest Houston, Texas experienced severe land subsidence (24-30 mm/yr) due to fluid extraction. They also found a significant correlation between areas that were inundated after Hurricane Harvey and observed land subsidence rates.
Unlike the previous studies that restricted PS-InSAR applications for non-vegetated areas because dense vegetation can lead to rapid decorrelation between SAR acquisitions [28,77], our results over the rural Texas Coastal Bend region showed a reasonable correspondence between the PS-derived and GNSS-derived land deformation rates. This research demonstrated the effectiveness of PS-InSAR techniques and their potential for quantifying land deformation rates in other rural and coastal areas around the globe.
This study has a few limitations. Some are related to the temporal coverage of the Sentinel-1 scenes used to map land subsidence rates and locations. Our study period was relatively short (October 2016-July 2019); however, combining these data with datasets from other satellite platforms (e.g., ENVISAT, PALSAR) makes it possible to extend this range to cover the entire 21st century. We believe that the short study period could be an issue in the sense that subsidence is a long-term problem, and we have only captured a snapshot of this with a 4-year study period. This means that the rates presented in this study-while accurate-may not be representative of the historical land subsidence rates.
Another limitation is the lack of bidirectional coverage of Sentinel-1 scenes within the study area. Within our study area, only ascending scenes were available that provided adequate coverage, and no descending scenes were available. Due to this limitation, we assumed that the movement of the PS was purely vertical with no horizontal movement. LOS velocities are the product of movement in three spatial axes (north, east, vertical), although primarily in the vertical direction. If two acquisition directions are available, it would be possible to use trigonometry in conjunction with the incidence angle of each pixel to subtract the north and east movement vectors, leaving the true vertical displacement that is made clear from the slight movements in the other two directions.
It is also challenging to compare the GNSS-and InSAR-derived land deformation rates in this study (R between the two datasets was estimated to be 0.24). GNSS stations are sparse and GNSS rates were calculated from daily measurements, while the InSAR rates reported in this study had a temporal baseline of roughly 1 month. Other factors that contribute to the inconsistency between the two rates include uncompensated atmospheric artefacts, orbital errors, unwrapping errors over low coherence areas, and/or inherent errors associated with GNSS observations. We faced two challenges in flood mapping: the lack of visible images required for accurate validation of flooding results, and the relatively large temporal separation between the storm time and next available SAR image. The high-resolution optical images from NOAA only exist in the Rockport area for Hurricane Harvey, and no imagery is available for Hurricane Hanna. The optical images for Harvey were acquired on 29 August 2017, which was 4 days prior to the available post-event image. We believe that this temporal difference makes the use of this optical image as a validation method less accurate. We also attempted to use optical satellite images from Landsat 8 and Sentinel 2 to verify the results of the flood mapping exercise; however, the images close in temporal baseline to the flood events were cloudy and extremely challenging to use for validation. In addition, there were no SAR images available in the immediate hours or days after the storm events. In the case of Harvey, the first available SAR image was 7 days post-landfall (3 September 2017). Therefore, it is likely that the areas that were shown as inundated did not show the full extent of the flooding that actually occurred within the study area. We believe that, due to the temporal difference in acquisition dates, our map underestimates the spatial extent of areas flooded by Hurricane Harvey, especially in the areas near the coast, where floodwaters would dissipate quickly. Despite this, land subsidence is believed to be a contributing factor in the areas becoming severely flooded following hurricane/storm events.
Our results will enhance the prediction of coastal communities' responses to natural forces (e.g., flooding, hurricanes, and sea level rise) and will eventually facilitate extremeevent mitigation and remediation. This, in turn, will improve the resiliency of Texas coastal communities. Decision makers could use these results to explore sustainable mitigation scenarios for subsided areas while preserving the unique nature of valuable coastal resources.
Mapping land subsidence rates and locations in South Texas significantly improves the current understanding of the factors that control the variations in land deformation in coastal areas by providing high-resolution spatial and temporal datasets that were previously not available on this scale. In addition, most models that are currently used to assess the impacts of sea level rise on Texas coastal resiliency apply land deformation rates that are constant in both space and time. The results of this study could be used to improve the performance of many sea level rise and/or storm impact models in Texas and other parts of the globe.

Conclusions
The objectives of this study were to map land subsidence rates and locations in the Texas Coastal Bend between October 2016 and July 2019, investigate the factors controlling the observed subsidence rates, and explain the consequences of the subsidence. We identified rapid local land subsidence, with spatial variability, in seven Coastal Bend regions: Victoria, George West, Refugio, Falfurrias, Karnes City, McAllen, and Nueces Bay/Corpus Christi. We suggest that the rapid subsidence rates in these regions are the result of both natural and anthropogenic factors, namely the heavy growth faulting dissecting the coastal plains of Texas along with the high fluid (hydrocarbon and groundwater) extraction rates seen in some counties. We were also able to map the areas flooded in the aftermath of two recent tropical cyclones and compare the subsidence rates in the flooded areas to those of the areas not affected by flooding in the two largest population centers in the study area. We found that there was a significant increase in the average subsidence rates within the flooded areas compared to those that were not flooded.
With the projected rise in sea levels coupled with both natural and anthropogenic induced subsidence, Texas Coastal Bend communities face an increased risk of flooding. The study area has already witnessed a significant number of tropical cyclones; however, land subsidence in this region amplifies what is already one of the most destructive forces on Earth-flooding. Our results signal the need to generate a complementary longer land deformation time series using InSAR platforms, such as ENVISAT and PALSAR, to better understand the spatial and temporal variabilities in land subsidence locations and rates, as well as the factors that control the observed subsidence rates. These InSAR platforms could also be used to better map inundated areas, which could help to mitigate the effects of future flood events. This study serves as a warning that subsidence is an ongoing issue in south Texas, and if not monitored closely it could have significant detrimental effects for the region in the future.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
This study was performed based on public-access data. No new data were created in this study. Data sharing is not applicable to this article. 28 8     Figure A1. Histograms and statistics of InSAR-derived land subsidence rates over regions/polygons 1−7 shown in Figure  6a.