Subsidence of sinkholes in Wink, Texas from 2007 to 2011 detected by time-series InSAR analysis

Abstract

West Texas’ Permian Basin, where the Wink sinkholes are located, is underlain by evaporite rocks that have been exposed to the dissolution due to natural processes and/or anthropogenic activities. We used the time-series interferometric synthetic aperture radar technique to process 16 ALOS PALSAR images from 01/01/2007 to 02/27/2011 for inspecting ground stability. Our results clearly show that two major sinkholes (Wink Sinks 1, 2), formed in 1980 and 2002, are both still subsiding, but the maximum subsidence for the 4-year period (2007–2011) occurred over an area about 1 km northeast of Wink Sink 2. The peak subsidence rate reached ∼40 cm/year during 2007–2010 and rose to ∼50 cm/year after 2010 when the area was hit by a record drought. Continuous monitoring of the subsidence in the vicinity of the Wink sinkholes is required for preventing and mitigating catastrophic outcomes of long-term developing geohazards to the area’s oil production facilities, infrastructure, and human safety.

Keywords:

• Wink sinkholes
• InSAR
• time series InSAR

1. Introduction

A sinkhole is a depression or a hole caused by the collapse of land surface (Gutiérrez et al. 2014). Such phenomenon is generally found in geological formations known as karst topography where soluble rocks (evaporites, carbonate rocks) dissolve as a consequence of chemical process and collapse (e.g. Beck and Pearson 1995; Johnson and Neal 2003; Waltham et al. 2005). Karst is often prone to the ground subsidence induced by the contact with freshwater from aquifer systems. Sinkhole formation can be related to dissolution of the water-soluble rocks, the development and expansion of underground cavities, consecutive roof failures, and crop-out on the surface. According to Gutiérrez et al. (2014), natural sinkholes can be classified into two types: “solution sinkhole” and “subsidence sinkhole”. A solution sinkhole is attributed to the surface dissolution, which creates shafts or conduits and can lead to ground instability. In this case, karst rock is exposed with almost no mantle soil. On the other hand, a subsidence sinkhole is created by ground surface subsidence, which may be attributed, for example, to chemical dissolution of caprock and bedrock underground, internal erosion of rock, subsurface dissolution and downward gravitational movement of the overlying material. It commonly shows one of the phenomena: collapse, sagging and/or suffusion (Gutiérrez et al. 2014).

Sinkhole hazards can be found in karst areas worldwide (e.g. Beck and Pearson 1995; Johnson and Neal 2003; Waltham et al. 2005). Well-known areas with sinkholes are located along the dead sea coastline (Baer et al. 2002; Abelson et al. 2003; Wust-Bloch and Joswig 2006; Nof et al. 2013), northeastern Spain (Galve et al. 2009; Castañeda et al. 2009; Gutiérrez et al. 2011; Castañeda et al. 2011; Yerro et al. 2014; Carbonel et al. 2015; Galve et al. 2015), the Netherlands (Chang and Hanssen 2014) and North America (Paine et al. 2012; Vaccari et al. 2013; Rucker et al. 2013; Jones and Blom 2014; Kim et al. 2016; Kim and Lu 2018).

Sinkhole development can be induced by natural processes and/or anthropogenic activities in a karst environment. According to Waltham et al. (2005), the vast majority of newly formed sinkholes are human-made. The anthropogenic factors primarily include (i) incremental water input into ground (cover and bedrock), (ii) aquifer exploitation and mining dewatering resulting in water level decline, and (iii) static and dynamic loadings such as water impoundment, vegetation removal or ground freeze and thawing. Sinkholes can cause economic losses, directly and/or indirectly. For example, sinkhole subsidence can severely destroy the integrity of infrastructure such railways, roads, canals, pipeline systems, buildings, and even nuclear power stations (e.g. Kuniansky et al. 2016). A critical step in sinkhole hazard analysis must therefore include an infrastructure inventory investigation and risk assessment.

Sinkhole-induced deformation and collapse have been observed by the following methods: microgravity (Paine et al. 2012), Light Detecting And Ranging (Filin and Baruch 2010; Filin et al. 2011), spaceborne InSAR (interferometric synthetic aperture radar) (Chang and Hanssen 2014; Kim et al. 2016; Kim and Lu 2018), field surveys (Bruno et al. 2008; Margiotta et al. 2012), Uninhabited Aerial Vehicle SAR (Jones and Blom 2014), and geological mapping. The detailed subsurface information can also be directly obtained through spelaeological exploration and geophysical field surveys through trenching, probing, and drilling. Incorporating all types of data and field geological investigation methods provides higher resolution sinkhole deformation maps and enables the characterization of the sinkhole-related deformation cause(s) and expected progression. Although considerable research has been increasingly devoted to the investigation of sinkhole activities in recent decades, there is still much room for improvement to further understand the mechanism of sinkholes, in particular, through continuous monitoring of developing sinkholes.

In the Permian basin of west Texas, two large collapse sinkholes emerged in Wink County in 1980 and 2002, dubbed Wink Sink 1 and Wink Sink 2, respectively (Figure 1). They are believed to still be actively developing, likely attributed to a deep supply-water well, since fresh or under-saturated saline water flowing through improperly cased wellbores and/or rock fractures may contribute to salt dissolution (Johnson et al. 2003). Paine et al. (2012) processed microgravimetry, GPS, and InSAR observations from a January to July 2007 advanced land observation satellite (ALOS)-1 interferogram and concluded that the maximum vertical subsidence was approximated as 30 cm/year. Kim et al. (2016) observed the development of Wink sinkholes using Sentinel-1A of May–August 2015 and found that the areas nearby the current sinkholes were also unstable and experiencing a large deformation. The subsidence east of Wink Sink 2 could be more than 40 cm/year during 2015–2017 based on Sentinel-1A/B InSAR imagery (Kim and Lu 2018). Due to the coarse resolution of Sentinel-1A/B and the rapid subsidence, the peak subsidence from Sentinel-1A/B could be even greater (Kim et al. 2016; Kim and Lu 2018). Recently, using high-resolution TerraSAR-X imagery, Kim et al. (2019) calculated the peak subsidence over the area about 1 km east of Wink Sink 2 and found it to be about 60 cm/year during 2015 and 2016. However, there has been no study regarding temporal development of the two sinkholes and surrounding regions from 2007 to 2011. Thus, the major motivation of this paper was to monitor the time-varying sinkhole development, which can be important to assess the state of deformation of these two sinkholes and mitigate potential sinkhole hazards nearby. More specifically, we have processed the time series of 16 HH-polarized L-band ALOS-PALSAR (Phased Array type L-band SAR) images taken between January 2007 and February 2011 and confirmed that the two sinkholes are still subsiding and the area to the east of Wink Sink 2 was already experiencing subsidence at an alarming rate during 2007–2011.

Figure 1. Study area of Wink sinkholes (black box) near Wink and Kermit, Texas that is covered by ALOS PALSAR data. Wink Sink 1 developed in 1980 is located ∼1.5 km north of Wink Sink 2, which outcropped in 2002. Background image is from Sentinel-2. Blue (N4616103) and yellow (USW00023040) triangles represent groundwater well and weather station, respectively. Source: Author

2. Data and methods

A synthetic aperture radar (SAR) system can obtain high-resolution images by emitting radar waves to targets from a spaceborne, airborne, or ground-based platform. InSAR technique combines two SAR images of the same imaging geometry from the same area (e.g. Rosen et al. 2000; Hanssen 2001). The phase differences between the images, or the interferogram, highlights even slight changes in the surface over time. Multiple processing steps are necessary to properly align the images and correct for signal noises and potential errors. The interferometric phase after the topography correction is proportional to a subtle terrain deformation, and contaminated by atmosphere changes between the two acquisitions, topography error, position uncertainty and other noises: (1) $${\varphi }_{}=-\frac{4\pi }{\lambda }\frac{B\text{cos}\left(\theta -\alpha \right)}{r_{}\text{sin}{\theta }_{}}{h}_{\text{err}}+{\varphi }_{\text{def}}+{\varphi }_{\text{atm}}+{\varphi }_{\text{noise}},$$(1)

where B is the length of the baseline, α is the baseline orientation angle, λ is the radar wavelength, θ is the look angle, r is the slant-range from the target to the reference satellite, h_err is the topography height error, and φ_def, φ_atm and φ_noise are, respectively, the components of deformation, atmospheric and noise in the interferogram phase. InSAR has been widely applied to various studies of geohazards, for example, landslides, earthquakes, volcanoes, and land subsidence (e.g. Massonnet et al. 1993; Amelung et al. 1999; Bawden et al. 2001; Lu and Danskin 2001; Lu et al. 2010; Lu and Dzurisin 2014). Decorrelation in space and time between the two consecutive SAR acquisitions can, however, affect the robustness of the results (Zebker and Villasenor 1992). Artefacts due to atmospheric and orbital errors could significantly degrade the accuracy of measurement (e.g. Ferretti et al. 2001; Li et al. 2005). Multi-interferogram techniques, including SBAS (Small Baseline Subset) InSAR, PSInSAR (Persistent Scatterers InSAR), and SqueeSAR, have been proposed to overcome these problems and retrieve time-series deformation histories (e.g. Ferretti et al. 2001, 2011; Hooper et al. 2004, 2007; Hooper and Zebker 2007; Hooper 2008; Lu and Zhang 2014; Qu et al. 2015).

The SAR data used in this study were provided by the ALOS that was launched by the Japan Aerospace Exploration Agency (JAXA) in January 2006. We utilized the HH-polarized Fine-beam mode PALSAR (Phased Array type L-band SAR) collectively spanning from 1 January 2007 to 27 February 2011 (Figure 2).

Figure 2. InSAR pair configuration of perpendicular and temporal baselines. Each circle expresses a SAR acquisition. Lines represent the generated differential interferograms. Perpendicular baselines of all InSAR pairs are calculated with respect to the reference image (first SAR acquisition).

All SAR images are co-registered to a single master image by computing the offsets in range and azimuth between each of the slave images and the master image and then resampling the slave image to precisely match with the master image. The calculation of the offsets between two SAR images is based on the local spatial correlation for a number of small windows (e.g. 64 pixels by 64 pixels) throughout the image; the optimum range and azimuth offsets that maximize the local correlation can be determined by the cross correlation of the SAR intensity images. The calculated offsets at many locations throughout the images are then used to define the coefficients of the polynomial that allows the resampling of the slave image to the master one. Initially, all interferograms with spatial baseline of less than 2500 m and temporal baselines of 12 months were generated. Inspection of these interferograms indicates the coherence over the peak of the deformation zone is significantly reduced if spatial or temporal baseline increases markedly. This is caused by the nature of the deforming zone: large deformation rate in time and high deformation gradient in space. After quality screening, a total of 18 interferograms with high coherence were used in this analysis (Figure 3). The topographic effect is removed using the 1-arc-second (∼30 m) SRTM (Shuttle Radar Topography Mission) DEM data. Based on all 18 interferograms, the residual topography components and orbital errors (including some large-scale atmospheric effects) were modelled and corrected by considering the relationship of residual phases with perpendicular baselines (residual topography phases) and two-dimensional (2D) phase variation as a polynomial form (orbital errors) (e.g. Ferretti et al. 2001; Hooper et al. 2004; Lu and Dzurisin 2014).

Figure 3. Unwrapped differential interferograms from 2007 to 2011, showing LOS range change (mostly subsidence). Each fringe represents a LOS range change of 11.8 cm. Source: Author

As the study area in West Texas is flat (topography variations are less than 20 m), atmospheric artefacts (caused by turbulent troposphere) in interferograms can be reduced by stacking (e.g. Shirzaei and Bürgmann 2012). We can then estimate the mean deformation phase rate as follows: (2) $${\varphi }_{ph\_\mathit{\text{rate}}}=\frac{{\sum }_{j=1}^N{\varphi }_j}{{\sum }_{j=1}^N\mathrm{\Delta }{t}_j},$$(2)

where N is the total number of the differential interferograms, φ_{ph_rate} is the mean phase rate, Δt_j is the jth time interval, and φ_j is the phase of the jth differential interferogram.

A least-squares method can be applied to invert for the time series of the deformation from a series of time-dependent differential interferogram phases as follows: (3) $$Gm=d,\mathrm{}G=\left[\begin{array}{c}\begin{array}{ccc}1& 0& \begin{array}{ccc}0& 0& \begin{array}{cc}0& \end{array}\end{array}\end{array}\\ \begin{array}{ccc}0& 1& \begin{array}{ccc}1& 0& \begin{array}{cc}0& \cdots \end{array}\end{array}\end{array}\\ \begin{array}{c}\begin{array}{ccc}0& 1& \begin{array}{ccc}1& 1& \begin{array}{cc}0& \end{array}\end{array}\end{array}\\ \begin{array}{ccc}0& 0& \begin{array}{ccc}0& 1& \begin{array}{cc}1& \end{array}\end{array}\end{array}\\ \begin{array}{ccc}& & \begin{array}{ccc}⋮& & \begin{array}{cc}& \ddots \end{array}\end{array}\end{array}\end{array}\end{array}\right],$$(3)

where m= [m₁, m₂,…, m_i,…, m_k−1]^T, m_i is the ith element of the incremental chronological range change between two adjacent SAR acquisitions, and k is the total SAR acquisitions; d= [d₁, d₂,…, d_i,…, d_n]^T, d_i is the observed phase of the ith interferogram, and n is the total number of interferograms (Figure 2); the element of G will be 1 if the interferogram spans the corresponding SAR acquisitions, otherwise 0. After m is obtained through an inversion of the matrix G, the cumulative LOS (line-of-sight) range change (deformation), can be calculated.

3. Results and discussions

The map of the mean deformation rate based on the stacking technique and the average deformation rates along two profiles are shown in Figures 4 and 5, respectively; the time-series deformation maps estimated from our simple SBAS InSAR method are shown in Figure 6.

Figure 4. Mean LOS deformation rates (cm/yr) around Wink sinkholes. Maximum LOS deformation rate reaches ∼30 cm/yr (∼40 cm/year subsidence when it is converted from LOS to vertical direction). Source: Author

Figure 5. Mean LOS deformation rate (cm/yr) along (a) P₁–P₂ and (b) P₃–P₄ in Figure 4.

Figure 6. Cumulative time-series LOS deformation (each fringe cycle represents a LSO range change of 11.8) around Wink sinkholes estimated from 16 ALOS-PALSAR scenes. Subfigures Ti (i = 0, 1, 2, …, 15) correspond to the cumulative LOS deformation images between the first (2007.01.01) and the ith acquisition date. We mark the number and the corresponding date as following (0: 2007.01.01; 1: 2007.07.04; 2: 2008.05.21; 3: 2008.07.06; 4: 2008.11.21; 5: 2009.01.06; 6: 2010.01.09; 7: 2010.02.24; 8: 2010.04.11; 9: 2010.05.27; 10: 2010.07.12; 11: 2010.08.27; 12: 2010.1012; 13: 2010.11.27; 14: 2011.01.12; 15: 2011.02.27). Source: Author

From January 2007 to February 2011 (∼4-year span), the cumulative LOS deformation was approximated as 125 cm over the area located about 1 km northeast of Wink Sink 2 (line D₅ on Figure 7(a)). Because the LOS deformation reflects the ground motion in vertical direction (subsidence) dominantly (Kim et al. 2016; Kim and Lu 2018) and the LOS deformation can be converted to the ground motion in vertical direction by using an incidence angle of 38.74°. The peak subsidence reached to ∼40 cm/year. Subsidence was also occurring in the existing sinkholes. Over Wink Sink 1, the cumulative LOS deformation was more than 30 cm (∼10 cm/year vertical subsidence). The region 500 m north of Wink Sink 2 shows 60 cm cumulative LOS changes (∼20 cm/year subsidence). However, the most troubling place in the Wink area is the region 1 km northeast of Wink Sink 2, experiencing an average subsidence rate of ∼40 cm/year.

Figure 7. (a) Time-series cumulative LOS deformation at four points (D₁, D₂, D₃, D₄). (b) Precipitation (black lines) and groundwater depth (blue diamonds) during ALOS acquisition period (January 2007–February 2011) at weather station (USW00023040) and groundwater well (N4616103), respectively. The precipitation data are from NOAA (National Oceanic and Atmospheric Administration) National Weather Service.

To delineate the temporal deformation patterns at the points of the deforming areas, we selected six subsidence points (Figure 4). The patterns of these deformations from 01/01/2007 to 02/27/2011 are shown in Figure 7. It is obvious that the subsidence was steady without seasonal fluctuations during 2007–2011, which represents an important signature of ground subsidence due to salt dissolution as suggested by Kim and Lu (2018).

The peak subsidence during 2007–2011 is around D₅, however, from 2015 to 2016 (Kim et al. 2019), it was slightly southeast of D₅, over the intersection between county roads 201 and 204 (Kim and Lu 2018; Kim et al. 2019). This suggests the peak of subsidence over the area east of Wink Sink 2 is migrating slightly south and east with time. The subsidence near Wink Sink 1 was ∼10 cm/year while it was ∼ 20 cm/year over the area north of Wink Sink 2 during 2007–2011. The subsidence over these regions during 2015–2017 was ∼5 cm/year (Kim et al. 2019). These suggest the subsidence over Wink Sink 1, Wink Sink 2 and the area north of Wink Sink 2 decayed with time, which is consistent with the stabilization of two existing sinkholes in term of their dimensions (Kim et al. 2019). Since the time of each of the two Wink sinkhole collapses, their respective underground cavities have continuously filled with the debris from upper formations (mostly sandstones). The void-filling process, called a suffusion (Waltham et al. 2005), results in the surface subsidence in early stage of sinkhole collapse and then the gradual surface stabilization as a consequence of nearly fully-filled cavity.

When compared with the linear 30 cm/year line (dotdash line in Figure 7(a)), the LOS deformation at the peak of subsiding area (D₅) accelerated from 30 cm/year before 2010 to ∼38.9 cm/year afterwards. Therefore, the peak subsidence accelerated from ∼40 cm/year (30 cm/year in LOS) to 50 cm/year (38.9 cm/year in LOS) during 2010 and 2011 (Figure 7) (note: the peak LOS deformation is generally not affected by horizontal components, because it is often located where east-west movements become nearly zero (Kim and Lu 2018). Therefore, the peak LOS deformation can be converted to the peak subsidence without considering deformation in horizontal direction). As the location of the peak subsidence migrated southward from 2011 to 2016, the rate of peak subsidence also rose from ∼40 cm/year during 2007–2010 (Figure 7(a)), to ∼50 cm/year in 2010 and 2011 (Figure 7(a)), to ∼60 cm/year during 2015 and 2016 (Kim et al. 2019). The acceleration could be due to the lower precipitation in 2009 and 2010 (Figure 7(b)). In fact, the 2011 drought was the worst in the past 30 years in Texas: the precipitation in West Texas was less than a quarter of the 30-year annual average precipitation (Kim et al. 2019). Drought events increase vulnerability to sinkhole development due to the increased overburden stress caused by the lowered groundwater level (Gutiérrez et al. 2014). The increased effective stress and internal erosion in the overlying layers can result in more rapid downward percolation of freshwater in aquifer systems into the underlying salt beds (Linares et al. 2017), finally causing the acceleration of the salt dissolution and surface/subsurface subsidence. Unfortunately, we do not have continuous groundwater-level measurements near Wink the sinkholes, however, annual measurements are available from a nearby well N4616103 (blue triangle in Figure 1; data on Figure 7(b)). Our InSAR results suggest that the subsidence near Wink sinkholes due to the dissolution of salt beds has accelerated and expanded after the 2011 drought (Figure 7).

4. Conclusion

InSAR techniques can indeed provide important information on sinkhole development. We have processed four years of InSAR data collected during 01/01/2007–02/27/2011so as to monitor the deformation of existing sinkholes and surrounding areas near Wink, Texas. The time series of the InSAR deformation clearly shows that the existing sinkholes are continuously subsiding and that conditions that could lead to new ones are developing at an alarming rate. Particularly, the area east of Wink Sink 2, which has shown dramatic subsidence of ∼60 cm/year during 2015 and 2016 (Kim et al. 2019) was already sinking at a rate of 40–50 cm/year during 2007–2011. The subsidence rate ramps up after 2010, as precipitation declined, particularly during the 2011 record drought. We conclude it is critically important to not only continue monitoring the area east of Wink, but also to develop a plan to mitigate potential disasters in case the ground collapses.

Acknowledgements

The ALOS PALSAR datasets were provided by the Japan Aerospace Exploration Agency (JAXA) via the Alaska Satellite Facility (ASF) Distributed Active Archive Center (DAAC).

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work is partially supported by the Natural Science Foundation of China (No. 41674013, 4187402), the China Geological Survey project “The Urban Geo-hazards Investigation in Lvliang mountaineous area” (No. 121201011000150018), and the China Scholarship Council.