Survey and visualization of land subsistence caused by groundwater depletion using Sentinel-1A IW TOPS Interferometry

The increasing demand for groundwater for agricultural, industrial and urban use is growing in many countries (Dehghani et al., 2013). For decades, continuous harvesting of groundwater aquifers has led to significant declines in plains. Iran is a country with wide spatial and temporal limitations regarding water resources, especially in the central arid areas. Reports suggest that in some areas of Iran the subsidence rate has exceeded 30 cm per year (Amighpey and Arabi, 2016).

More than 25 percent of groundwater aquifers of Iran are in critical condition due to over-exploitation (Bamler and Hartl 1998). Although the high correlation between land subsidence, the reduction of groundwater level and the change in the mechanical properties of the subsurface layers have been known for years and many attempts have been made to fully understand this phenomenon, but no comprehensive and precise model of prediction of land subsistence is available for Iran. By reducing the level of groundwater, an increase in effective stress caused by decrease in pore water pressure will cause subsidence. Typically, reduction of the pore pressure and increase of the effective stress occurs over long time periods; therefore, reduction of the piezometric level and the associated subsidence show up with a delay (Dehghani et al., 2013). The subsidence is mainly due to the loss of the soil surface fluid that is mainly happening in unconsolidated or semi-consolidated sediments adjacent to the sand layers (Salehi et al., 2014). Under such conditions, non-elastic congestion is caused by an increase in the effective stress in the soil while the formation of the soil particles is reduced and the thickness of the vertical layer decreases. The uneven drop in the groundwater level and the heterogeneity of the aquifer texture cause anomalous subsidence, and consequently, lines and gaps in land surface, especially agricultural lands, are observed. The cracks of the earth spread silently and cancerously that lead to more abnormal conditions by altering the aquatic status of the area such as the direction and speed of groundwater flow (Qu et al., 2017).

The high subsistence rate destructs structures such as buildings, roadways, bridges, transfer pipes and earth dams (Hoffmann et al., 2001).  For many years in Iran, where the phenomenon of subsidence was not so commonplace, researchers had realized this phenomenon by observing the outflow of water from the agricultural water wells.  In recent years, extensive efforts have been made to measure changes in the earth's crust using techniques such as precise alignment. During the past two decades, two advanced GPS technology and radar imaging have been used to increase the accuracy of surface deformation measurement using satellite-powered platforms. In the meantime, GPS systems with two features include low horizontal spatial resolution and high-resolution time have provided useful tool for examining small changes in land surface displacement. On the other hand, the coverage and high spatial resolution has transformed satellite radar images as one of the strongest tools for researchers in the field of crust changes such as earthquake cycle modeling, geothermal activities, landslide zoning (Dai et al., 2016), vertical land subsidence measurement due to underground water exploitation, modeling of tectonic movements of faults (Qu et al., 2017), and studying of active salt domes (Abdolmaleki et al., 2014) and volcanoes. Interferometry using radar satellite imagery with the artificial window or synthetic aperture radar (SAR) provides a precise method based on the use of at least two radar images of an area and is capable of measuring the vertical and horizontal displacement with remarkable precision in wide areas with 0.1 to 1 mm accuracy in different intervals (Tung et al., 2016).

The InSAR method is possible for estimating large-scale vertical displacements. The basis of this approach is that if two radar images of a region are taken at two times, the interferogram to illustrate the combination of these two images will show slight changes in the earth's surface on a mm scale. When the surface of the earth in a region is captured twice with the same geometry, if the radar-earth distance is the same in both images, the backscatter phases will remain the same, which means that no altitude difference has been made. However, if the level in the earth's crust is displaced (horizontal or vertical), the radar-earth distance changes and the phase of the second image will be shifted. This phase shift signifies the displacement of the surface. The corresponding motion rate is half the length. In order to measure the vertical displacement of the surface of the earth, the Interferometric synthetic-aperture radar (InSAR) has been considered as one of the non-geodetic methods in recent years (Wang et al., 2018).

Dehghani et al., (2009) studied Mashhad Plain in the period 2003-2005 using InSAR and Envisat images and showed that the subsidence of the plain was 23 cm per year in most cases. The researchers also investigated the piezometric wells of Mashhad Plain and showed that in many wells, the piezometric surface loss is correlated with the Earth's surface; however, in some wells, despite the piezometric loss, the subsidence signal was not observed. Motagh et al., (2017) in a study of Rafsanjan Plain, used InSAR time series analysis obtained by the exploitation of Envisat, ALOS and Sentinel-1 (S1) SAR data archives to assess land subsidence in an area of 1000 km². Their results showed that in the recent years, the subsidence rate of the plain had been changed from 5 cm/year to more than 30 cm/year. This event follows mainly from the extreme exploitation of groundwater resources and also is partly influenced by Quaternary faults in the area.

Dai et al., (2016) in their interactive survey, reviewed Sentinel images of 2014 to 2015 in Mexico City and found that the annual subsidence rate is more than 40 cm per year. Ghazifard et al., (2017) studied Damghan Province between 2005-2011 and evaluated annual subsidence rate at 7 cm. They also found a significant relationship between the occurrence of subsidence and the degradation of groundwater piezometric levels and the reduction of fine-grained sediment thickness. Zhu et al., (2015) estimated the subsidence rate to be more than 52 mm per year in a study conducted on the North Beijing Plain. They did not find a meaningful relationship between urban development and subsidence rates. Hey-man et al., (2017) used Radarsat, Sentinel and ALOS images to evaluate subsurface underground subsidence, and concluded that if access to up-to-date imagery is possible for all three sensors, Sentinel images have more capability and accuracy than the other imagery.

Strozzi et al., (2017) studied the consequences of subsidence associated with construction of a 57-kilometer-long Gotthard tunnel in Switzerland and concluded that even without accurate information, Sentinel-1A imagery could be used to map the high-risk subsidence areas.

The Sentinel-1A images are available through the European Space Agency beginning from 14th January 2014 with a set mission of 7 years. This satellite has Sun-synchronous orbit with an orbital height of 693 km, and provides imaging in the C band with a wavelength of 55.5 mm with 12-day repeat cycle. The most important product of Sentinel-1A is the SLC data with a spatial resolution of 20 × 5 m and is able to show the displacement with precision of millimeters in interferometric studies. In the present research, Sentinel-1A data is used to survey and zonate the subsidence of Garmsar Plain, southeast of Tehran.

Materials and Methods

The study area

Garmsar Plain with an area of 320 km² is located between the geographical coordinates of 52˚15'N to 52˚35'N latitudes and 35˚05' E to 35˚17' E longitudes on the southern margin of central Alborz and in the west of Semnan province 90 km southeast of Tehran and in the north of Iran's central desert (Fig. 1). The average height of the plain from the sea level is 875 meters. Hablehrood River after leaving the south of the Alborz Mountains, forms the Garmsar plain at the foot of this mountain and in the north of the desert. Garmsar is part of Iran's central plateau and has hot and dry summers and cool and dry winters. Garmsar alluvial plain has always been affected by water scarcity due to its dry climate. However, it is one of  the most important plains in the country in which traditional agriculture has been modernized in recent years. Due to its proximity to Tehran, Garmsar has a special economic and industrial landscape. Over the past two decades, rapid industrial and agricultural development in the region has been accompanied by rapid population growth and, consequently, an increase in water demand, leading to an increase in the withdrawal of groundwater resources annual dropping groundwater aquifers level by 1.3 m (Azareh et al., 2014).

Figure 1. The location of the study area

Research Methodology

A radar signal in synthetic aperture radars (SAR) consists of two parts: amplitude and phase. The phase is a fraction of a complete cycle from the sinusoidal wave and amplitude is the strength of the return signal., In the microwave wavelength range that radar waves are also part thereof, the dielectric coefficient determines the conductivity and reflectivity of the waves. The dielectric coefficient of the material in the dry state varies from 3 to 8, while the dielectric coefficient of water is about 80. This property of radar waves has made them useful in natural resource studies. The phase in the radar image is essentially determined by the distance between the sensor antenna and the ground fault. Although the phase of a single SAR image gives us useful information, but the phase difference of a couple coherent and co-registered images in two different times provides valuable information for researchers. In radar interferometry, the phase difference between the two images is determined from a region with time and the specified baseline. In Figure 2 and Equations 1 to 3 the radial picture pair phase calculation method and the relationship between the phase differences between the two parameters have been shown. In Equations 1 to 3 and 5, the value of λ or wavelength in the Sentinel-1A is 55.6 mm.

           (2)

(3)

Figure 2. Geometry of SAR Interferometry (Bamler & Hartl, 1998)

In general, the word interferometry is used to study the concept of combining a pair of signals in order to obtain information about their difference. The basis of the interferometry method is on the difference between two images of an area. In the process, the wave ranges of two images are multiplied and the formed interferometer will have a phase difference compared to the initial image quantified as follows:

(4)

Where φ∆_orb represents the error introduced by orbital differences and the φ∆_flat is the error due to imprecise Earth's ellipsoidal surface. φ∆_atm  corresponds to the difference in atmospheric propagation times between the two acquisition used to form the interferogram. φ∆_noiseis the difference between the temperature of the satellite
sensor antenna, the volume distribution of the objects and the total effect of the noise component. φ∆_atmis the effect of the orbital component. φ∆_elv is the difference of the phase due to the topography of the region and is an important component of the interferogram related to the spatial basis between the two radar images (line B in Fig. 2). According to Figure 2, the topographic phase can be calculated by Equation 5 as below (Bamler and Hartl 1998):

          (5)

In equation 4 φ∆_dist is the difference in phase due to the displacement and the deformation of the earth, and the unknown part of this equation. In fact, the main objective of interferometry is to extract the amount of phase of displacement by eliminating or minimizing the effect of other components. In satellite radar systems, the waves are polarized. This means that forward and reverse waves have a direction relative to the perpendicular direction of propagation. In polarimetric radar system with two transmitted and received waves in two vertical (V) and horizontal (H) forms, four modes of VV, VH, HV and HH combinations are created. The two combinations of VV and HH are more suitable for studies of subsidence (Mullissa et al., 2017). In this study, two pairs of descending and ascending tracks in single look complex (SLC) format with the highest degree of co-integration with VV polarization were used (Table 1).

Table 1. Details of used SAR images in the present research

Before interferogram formation, we should ensure coherence between the two images. The coherency test of two radar images is one of the best methods for assessing the quality of the interferometry phase (Mullissa et al., 2017). In the first step, two images must be contacted. For this very important step, the orbital parameters of both images, along with the region’s DEM and permanent dispersion points were used. In this step, the newer image was selected as the master and the older was taken as the slave image. Then, using the sampling points, the image was re-sampled relative to the reference image. After correspondence and aligning references of the two images, the phase of the two images was combined (Fig. 5). The combined phase has a high correlation with the topography of the region and the surface displacement of the earth. After the formation of the initial interferogram, using DEM with a resolution of the 30-m, the phase of the topography was removed. Using the Goldstein filter, the noise level of the interferogram decreased and the accuracy increased (Goldstein and Warner 1998). This filter largely modifies the atmospheric effects. In Equation 4, the expression (2πα) is the remainder of the phase difference equation and using the unwrapping fan, this amount is eliminated from the final phase.

• Phase conversion to line of sight (LOS):

Due to the geometry of SAR sensors, the displacement is calculated by interferometry along the line of the satellite view (Tofani et al., 2013). In order to convert the interferometric phase to line of sight (LOS) due to subsidence, the following Equation 6 is used:

      (6)

where d is the amount of displacement along the satellite's view taken in the time interval of the two images per meter, λ is radar wavelength, inc_angle is the angle of inclination of the radar wave, and Ø∆ is the remaining phase of the subsidence resulting from Equation 4. Finally, the distortions in radar imaging geometry are corrected using the Range-Doppler method and converted to geographic coordinates.

• Converting line of sight displacement (LOS) to vertical displacement:

After the displacement in the direction of the satellite view was calculated for each pair of ascending and descending images, the vertical displacement component using geometric combinations of the two view directions (Fig. 3) was extracted using Equation 7 (Dai et al., 2015):

    (7)        

where V_A is displacement in the direction of the ascending circuit and V_D is displacement in the direction of the descending circuit, V_V and V_E are vertical and horizontal components of the displacement.   and  show the descending angle of the ascending and descending wave respectively.

Figure 3. Three-dimensional displacement field geometry using combined

ascending and descending images

Results and Discussion

After selecting two pairs of SAR images (Table 1), a preliminary interferometer was prepared. To ensure greater accuracy of the results, areas with a coherent rate of less than 0.5 were eliminated. These areas were not sufficiently correlated and generally show great variation. Basically, areas with vegetation, snow or with high humidity patterns have a low level of ccoherence and urban, desert and rocky areas have a high level of ccoherence. The coherence of the two pairs of images of the studied area due to desert and the poor vegetation cover was good. The final results reflected that Garmsar Plain has a minimum and  maximum 33 and 37 cm subsidence respectively, and the average annual subsidence of this plain is estimated to be 36 cm (Fig. 4). As shown in Figure 5, the profile of changes of subsidence in ab line indicates that as we move from the southeast to the northwest, the subsidence rates increases. Previously, Azareh et al., (2014) investigated the spatial variations of underground water in the Garmsar Plain and concluded that when we move from northwest to the southeast, the decline in groundwater levels decreases. This is consistent with the pattern obtained from the subsidence of the Garmsar Plain (Fig.4 and fig.5). Other researchers also achieved a similar pattern of connection between land subsidence and groundwater harvesting (Gazifard et al., 2016).

Figure 4. The subsidence rate of Garmsar Plain

Figure 5. The general trend of vertical displacement changes along the ab

line (Northwest to Southeast).

Conclusion

In this study, the Garmsar Plain subsidence was investigated using Sentinel-1A images between 2014 and 2015. Using the DInSAR technique, the average of the plain subsidence was evaluated at 30 centimeters. After zonation of the subsidence in this area, the drop level in groundwater was measured. The observed data from at least 10 piezometric wells in the region over the past 5 years showed a drop of 1 to 1.5 meters for groundwater. The level of groundwater drop in Garmsar Plain was estimated to be 1.5 m / year in the northwest of the plain and 1 m / year in the southeast. Comparison of the final subsidence zonation map of the area and the groundwater level indicated a close correlation between subsidence and groundwater drop.

Haghighatmehr et al., (2011) estimated the subsidence rate of the Hashtgerd Plain using interferometry and spatial information systems to be 56 centimeters per year, and stated that this phenomenon results from the high amount of harvesting underground water and the heterogeneity of input and output of this plain. Data shows substantial decrease in the flow of Hablehrood River during the years 1978 to 2007; therefore, the reduction of input in the Garmsar Plain on the one hand and the increase in groundwater harvesting on the other, has inevitably increased the pressure on this dry plain rising the subsidence rate of Garmsar Plain to 36 centimeters per year. The Garmsar geological structure entails conglomerate, marl layers, masses of gypsum and salts being extensive in the west of the plain. The effect of dissolution on salts and gypsum units exacerbate the severity of the subsidence phenomenon.

Therefore, it can be said that part of the subsidence in this region is due to its geological conditions. The delineation of the amount of subsidence due to geology and groundwater harvest is a matter to be considered in the future research. It is recommended to use this technique in the critical plains of middle east countries at short intervals to determine the behavior and interaction of the plains with regard to management plans. The management of water resources in order to compensate for the deficit of the groundwater reservoir, preventing of over discharge, blocking unauthorized wells, applying supportive and incentive policies, using new irrigation methods based on cropping patterns, installing smart meters and volumetric meters are some of the strategies that can be used to prevent the subsidence of Garmsar Plain and others suffering the same fate.

 Acknowledgements

The present study was supported by Iran National Science Foundation, Vice-Presidency for Science and Technology (grant number: 96016471).  We would like to thank Dr. James W. LaMoreaux and two anonymous reviewers for their very insightful comments, which greatly improved the quality of the original manuscript.