- Email: [email protected]
Studying land subsidence in Yazd province, Iran, by integration of InSAR and levelling measurements
Author’s Accepted Manuscript Studying land subsidence in Yazd province, Iran, by integration of InSAR and levelling measurements Masoome Amighpey, Siyavash Arabi www.elsevier.com/locate/rsase
PII: DOI: Reference:
S2352-9385(16)30020-9 http://dx.doi.org/10.1016/j.rsase.2016.04.001 RSASE23
To appear in: Remote Sensing Applications: Society and Environment Received date: 10 August 2015 Revised date: 11 April 2016 Accepted date: 12 April 2016 Cite this article as: Masoome Amighpey and Siyavash Arabi, Studying land subsidence in Yazd province, Iran, by integration of InSAR and levelling measurements, Remote Sensing Applications: Society and Environment, http://dx.doi.org/10.1016/j.rsase.2016.04.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Studying land subsidence in Yazd province, Iran, by integration of InSAR and levelling measurements Masoome Amighpeya, Siyavash Arabib a
National Cartographic Center of Iran, Meraj Ave., Azadi Squ., Tehran, Iran, [email protected], Tel.: +98 21 66071001
b
National Cartographic Center of Iran, Meraj Ave., Azadi Squ., Tehran, Iran, [email protected]
Abstract The Yazd-Ardakan plain, located in central Iran, is an arid region with a low rainfall where over-extraction of groundwater for agriculture and industrial purposes has caused serious subsidence in the area. In this study, the small baseline subset (SBAS) algorithm was applied to perform an interferometric synthetic aperture radar (InSAR) time series analysis for the Yazd-Ardakan plain. The result is a space-time deformation product that can be exploited to view not only the smoothly varying long-term surface motion, but also its time-varying patterns. The mean rate of subsidence computed from the repeated first order precise leveling network of Iran was accurately matched by the mean rate of the InSAR time series. Time series analysis of our InSAR mapping suggests that subsidence occurs within a northwest– southeast elliptic-shaped bowl, with peak amplitude of about 120 mm/yr for the 2003–2006 time period. Analysis of piezometric records suggests that subsidence likely results from extensive over-drafting of the aquifer system in the region, which has caused about 16 m of water table decline since 1974. Comparing the subsidence map produced by InSAR and water level change contour allowed the estimation of skeletal storage coefficients, providing a basis for future work that could characterize the storage properties of an aquifer system with a high degree of spatial resolution. The storage coefficient contour map generated in this study simulated the pattern and orientation of the geological structure of the area. 1 Introduction Iran is a country with wide spatial and temporal limitations regarding water resources, especially in the central arid areas. This problem has become more conspicuous in recent decades due to economic and urban development. Economic and population growth have deteriorated available water resources in terms of both quantity and quality aspect. The significant amount of water obtained from groundwater resources is such that land subsidence has become an issue in some areas (Amighpey et al., 2006, Motagh et al., 2007, Amighpey et al., 2008, Motagh et al., 2008, Dehghani et al., 2009). Incidentally, many of these subsidence regions are located in the precise leveling routes of Iran’s National Cartographic Center (NCC) network, which were observed twice in the last four decades and then assessed to determine potential height changes. In fact, these assessments were the first to detect the rate and location of this subsidence (Amighpey et al., 2006). Yazd province, located in central Iran, is one of these areas in which assessments of repeated precise leveling observations showed an extreme degree of subsidence, especially in the Yazd-Ardakan flat area (Amighpey et 1
al., 2006). The Yazd-Ardakan desert, with a length of 60 km and width of 15 km, is located between Yazd and Ardakan cities. Over the past two decades, subsidence and rupture of the earth have led to significant damage in Yazd province. This phenomenon has caused roads to rupture, water systems to break and structures to crack. In fact, one of the main considerations of construction work in Yazd is the numerous instances of earth rupture between Yazd and Ardakan. This phenomenon is the direct result of subsidence. Being located beside the central mountains, far from sea, adjacent to the desert and in a low rainfall region, Yazd has a climate which mostly resembles dry desert climates. Little rain along with high water evaporation, relatively low dampness, considerable heat and great temperature changes are among the factors making this province one of the driest parts of Iran, yet water needs for consumption, agriculture and industry are met mainly by groundwater resources. This has brought about over-extraction of groundwater in recent decades, which is, in turn, the main reason of subsidence in this area. Time series analysis of satellite radar interferometry is a useful method for studying both spatial and temporal patterns of surface deformation and has been extensively used to study subsidence areas (Berardino et al., 2002; Usai, 2003; Schmidt and Bu¨rgmann, 2003). In the present study, the InSAR deformation time series, based on the small baseline subset (SBAS) algorithm (Berardino et al., 2002), was computed for the Yazd-Ardakan desert, using data from the Envisat satellite acquired from the European Space Agency (ESA). We found that the mean rate of InSAR time series analysis agreed well with re-leveling results later used to measure subsidence. The availability of such a joint temporal/spatial deformation data set, when combined with independent hydrologic and geologic assessments, can provide more information on the storage properties of the aquifer system. A better understanding of the areal variability of the aquifer system response to stress can improve the effectiveness of groundwater management schemes and the identification of zones with a high potential of fissure formation, which is valuable information for urban planning. In this study, a simple temporal model for the subsidence was obtained, one capable of also depicting its spatial extent. The spatial and temporal behavior of the subsidence was correlated to water-level fluctuations to characterize aquifer-system response patterns and to further refine measurements of the aquifer storage properties. Analyzing the temporal behavior of subsidence and groundwater levels in the region made it possible to distinguish high-risk areas, where high water over-drafting has caused linear subsidence, from moderate-risk ones, where water over-drafting is partially compensated for in the region’s rainy season. This could prove to be vital information for the management of water pumping in the area so as to prevent further damage. 2 Subsidence in the Yazd-Ardakan plain from precise leveling measurements Geodetic observations of local land subsidence for the Yazd-Ardakan region were first documented using precise leveling data collected by the National Cartographic Center of Iran (NCC) during 1993-2001. The pink color points in Fig. 1 (a) show the location of leveling benchmarks extracted from first-order precise leveling lines along the Naeen-Yazd Highway. A number of stations were destroyed between 1993-2001.
2
Figure1. (a) SRTM DEM of Yazd subsidence area. Green points denote the first order leveling routes and pink points denote parts of these routes that show subsidence. (b) Results of first-order precise leveling surveys in the Yazd-Ardakan region between 1993-2001. Error bars represent the standard deviations in cumulative random error from the reference benchmark (BM15).
Analysis was done on the development of vertical deformation along the Naeen-Yazd Highway between the 1993 and 2001 surveys relative to a datum at BM15 and BM17, near Nogonbad and outside of the subsidence area. Figure 1(b) depicts the height changes along these routes. For the 8-year period spanning1993–2001, the maximum cumulative degree of subsidence recorded by a benchmark is as high as 95.00 cm 0.92cm at a point located between Shamsi and Hojjatabad (BM62). Ardakan (BM50) and Yazd city (BM83) also subsided by 47.00 cm 0.79 cm and 22.00 cm 1.1cm during the same 8-year time span, respectively. Precise leveling accuracy was computed by means of the Vignal formula (Bomford, 1971). Data pertaining to temperature and air pressure are not available for this route, which could have been used to remove the refraction effect, but the fact that there was no correlation between vertical displacements and topography makes it evident that elevation-dependent systematic error has had no impact on the leveling data. InSAR time series analysis is a better resource for studying subsidence versus two existing epochs over a long interval of leveling observation. This temporal/spatial map can be particularly useful in providing good input for modeling subsidence in the consolidation theory. In the following, the results of an investigation of the spatial and temporal behavior of Yazd-Ardakan subsidence patterns using surface displacement maps generated by InSAR are presented. 3 InSAR measurements of subsidence and water level change assessments The concept of combining InSAR information from a large number of SAR Images, which can facilitate the solution of a deformation time series, is a relatively new one (Ferretti et al., 2000; Berardino et al., 2002; Usai, 2003; Mora et al., 2002; Lanari et al., 2004). Using a multi-temporal technique not only allows the estimation of discrete deformation events, but also provides an avenue for a temporal study provided that the whole time series is available. To measure ground subsidence in the Yazd-Ardakan region, we selected 10 Envisat SAR images (see Table 1) in descending orbits. Employing the so-called two-pass technique (Zebker et al. 1994), a set of 43 differential interferograms was generated covering a time period between March 26, 2003 and December 20, 2006. Interferometric processing was done using the public domain InSAR processor 3
DORIS (Kampes & Usai 1999). In this method, the interferograms are corrected for the phase signature due to orbital separation using precise DEOS satellite orbits, and for the topography using a 3-arcsecond SRTM digital elevation model (more information available at http://srtm.usgs.gov). The remaining orbital error was then removed by fitting a plane in the far field area where the leveling data showed no deformation and subtracting this plane from the interferograms. Areas with low correlation were filtered out from the interferograms. The availability of both time and space data makes it possible to effectively identify and remove atmospheric artifacts in the results via stacking and the use of the time dependent characteristic of said atmospheric artifacts. Using multiple independent interferograms helps find atmospheric artifacts where fringes existed in special interferograms only, as this effect depends on the acquisition time. It is also worth noting that precise leveling was used to test the accuracy of the computed deformation map. In the end, only 19 interferograms which enjoyed a high correlation coefficient were chosen for time series analysis (Table 2). Figure 2 shows the spatial vs. temporal baseline of these interferograms. Figure 3(a) shows the interferogram generated from images acquired on September 17, 2003 and September 1, 2004. As the interferograms represent the deformation in the time interval between the acquisition dates of master and slave images, least square adjustment was used to invert deformation in different time intervals (interferograms) to a time series of deformation indicating displacement at each date of acquisition. If D is the displacement vector of the acquisition date of the images, which is unknown, and I is the displacement vector between the slave and master images acquisition interval (interferograms), we have the following equation: AD = I (1) Where A is the coefficient matrix. Applying least square adjustment, we have: D ( AT PA) 1 AT PI (2) CD ( AT PA) 1
(3)
Where P is the inverse of the observation variance matrix, while C D is the variance matrix of the unknowns. Table1. Details of the used images in this study Sensor Track Pass Swath Date Orbit ASAR/IM 20 Descending I2 2003-03-26 5585 ASAR/IM 20 Descending I2 2003-09-17 8090 ASAR/IM 20 Descending I2 2003-11-26 9092 ASAR/IM 20 Descending I2 2004-06-23 12098 ASAR/IM 20 Descending I2 2004-09-01 13100 ASAR/IM 20 Descending I2 2004-12-15 14603 ASAR/IM 20 Descending I2 2005-03-30 16106 ASAR/IM 20 Descending I2 2005-10-26 19112 ASAR/IM 20 Descending I2 2006-02-08 20615 ASAR/IM 20 Descending I2 2006-12-20 25124 Table2. Details of the used interferogram in this study Interferogram Interferogram name Temporal baseline Baseline length (m) number (day) 1 20030326_20031126 245 118 2 20030326_20030917 175 434 3 20030326_20040623 455 153 4 20030326_20040901 525 403 Mission Envisat Envisat Envisat Envisat Envisat Envisat Envisat Envisat Envisat Envisat
4
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
20030326_20041215 20030917_20040901 20030917_20051026 20031126_20040623 20031126_20040901 20031126_20041215 20031126_20061220 20040623_20041215 20040623_20061220 20040901_20051026 20040901_20061220 20041215_20060208 20041215_20061220 20050330_20060208 20051026_20061220
630 350 770 210 280 385 1120 175 910 420 840 420 735 315 420
160 32 127 35 521 42 271 7 306 159 250 427 314 147 409
Figure2. Spatial vs. temporal baseline of the 19 interferograms used in this study.
To compute the InSAR time-series, 19 interferograms were inverted using the least square adjustment method, weighting each data set by the inverse of its far field variance (Funning et al., 2005). 9 deformation maps were obtained into March 26, 2003, the image acquisition date prior to which no deformation is supposed (the first image acquisition date used as the reference epoch). Figure 3(b) shows the mean rate of the subsidence in this area obtained from stacked interferograms. Looking at figure 3(a), the interferogram belonging to the period between September 17, 2003 and September 1, 2004, 3 fringe patterns are evident. The subsidence pattern can be divided into 3 parts: the largest rate occurring in an extensive spatial distribution between Meybod and Zarch shaped like a bowl and directed northwest– southeast along the axis of Yazd-Ardakan road, the second to the west of Ardakan, and the third in the south of Yazd city, capital of Yazd province.
5
(a)
(b)
Figure3. (a) Interferogram generated from images acquired on September 17, 2003 and September 1, 2004, (b) Mean rate of subsidence obtained from stacked interferograms (m/yr). Contour lines represent the water level change rate in m/yr, derived from piezometric wells (black triangles) in the 2003-2007 time interval.
The highest rate of subsidence in the Meybod-Zarch region can be found between Shamsi and Hojjatabad (120 mm/yr). Figure 4(a) shows the time series of the maximum subsidence area. As is apparent, the subsidence demonstrates a linear behavior with the passage of time. R-Square and RMSE (Root Mean Squared Error) statistics were used to evaluate the degree of linear fit (cf. Mood et al., 1974). Rsquare is the square of the correlation between the response values and the predicted response values and is used to measure how successful the fit is in explaining variation the data. . It is also called the square of the multiple correlation coefficients and the coefficient of multiple determinations. RMSE is known as the fit standard error. The R-square of the linear fitness in the time series is 0.97, meaning that the fit explains 97% of the total variation in the data. Also of note is that the RMSE of linear fit is 10 mm. Figure 4(a) also depicts the fitted regression line.
6
(a)
(b)
(c) Figure 4. InSAR time series of subsidence in: (a) maximum subsidence area between Shamsi and Hojjatabad, (b) Ardakan, (c) Yazd. Squares are vertical displacement at the image acquisition time. Error bars represent the estimated accuracy of vertical displacement. The blue line is the fitted regression line of vertical displacement.
As evidenced by the data, the highest rate of subsidence occurred in agricultural areas (where groundwater extraction is common) and caused the land to rupture and buildings to fracture. Figure 5 illustrates an instance of the impact of subsidence in the maximum rate area (Hojjatabad). The earth fissure orientation is perpendicular to the northwest–southeast direction of the subsidence.
Figure 5. Field observation of the subsidence effect on the maximum subsidence area in Hojjatabad.
Another bowl-shaped pattern of subsidence was detected in Ardakan. Figure 4(b) depicts the relevant time series and linear fitted model. The assessment of fitting linear regression reveals linear subsidence with a maximum rate of 100 mm. The Rsquare value of the linear fitness to Ardakan time series is 0.98, meaning that the fit explains 98% of the total variation in the data about the average. The RMSE of linear fit is 7 mm. This subsidence area is located in an agricultural land and industrial region. Yazd, the capital of the Yazd province, showed distinct fringes with a lower rate compared with Ardakan and Shamsi. These fringes appeared near the railway in Yazd city, around an industrial estate. Results of the assessment of the degree of fit yielded an R-square value equal to 0.5747 and RMSE of 20 mm, rejecting a simple linear 7
pattern and emphasizing seasonal effects )figure 4 (c)). Figure 6 shows the time series of deformation in the city. Yazd’s time series oscillations show instances of strong seasonal uplift in the area in the spring of each year due to the effects of annual artificial recharge (e.g. figure 6 (f) and 6 (h)). The highest subsidence rate is about 65 mm/yr.
Figure 6. Time series of elevation deformation in Yazd city.
To assess the accuracy of InSAR measurements, the subsidence rate acquired from leveling observations was projected on the satellite line of sight and then compared with the average velocity from the InSAR time series. Although the two leveling observations belong to 1993 and 2001, whereas the radar images were acquired in the 2003-2007 time interval, (in different time periods), comparing the mean subsidence rate obtained from these two distinct sources yields a good degree of correspondence in leveling benchmarks (figure 7). The existing RMS between them is 15 mm. This correspondence, together with the accuracy of the InSAR results, proves the existence of a constant subsidence rate over these years. The existence of such a constant subsidence rate indicates a lack of water pumping management in recent years.
Figure 7. Comparison between mean subsidence rate computed from leveling projected on the satellite line of sight and InSAR in the path of leveling.
8
To better understand the cause of land subsidence in Yazd-Ardakan plain, variations in groundwater levels in the piezometric wells in the region were studied. Groundwater level measurements have been regularly made in the region every month since 1973. The contour lines in Figure 3(b) show the behavior of the mean rate of water level decline in this zone and around the subsiding area. It is evident that a continuous lowering of the hydraulic head below the ground level has been occurring in the region over recent decades. For the time interval of recorded piezometric measurements from 1973 to 2007, the mean water table in this area declined about 15 m as a result of intensive groundwater pumping and deficient natural recharge. In principle, the over-extraction of groundwater leads to a reduction in water pressure within the aquifer system which is reflected by an increase in the effective stress and compaction of sediments (Terzaghi 1925). The relationship between these two processes depends on many factors, including geotechnical parameters such as sediment porosity, permeability and compressibility as well as the spatial distribution of the thickness of compressible layers within the aquifer system and the history and amounts of pressure drawdown caused by groundwater decline (Freeze & Cherry 1979). Some external processes such as tide, recharge, or regional tectonics might also affect the deformation of an aquifer system (Amelung et al. 1999; Moreau et al. 2006). As evidenced by Figure 3(b), although there is a clear correlation between the rate of water level reduction and the subsidence, the areas with the highest degree of aquifersystem compaction are offset from the zones of maximum water-level decline in all 3 subsidence areas: the Meybod-Zarch region, Ardakan and Yazd. Such offset between the areas of greatest subsidence and most significant water-level decline has previously been reported in the Las Vegas metropolitan area (Bell et al. 2008). This might be due to the higher spatial resolution of the earth surface deformation map produced by InSAR in contrast with the water level changes contours provided by the sparse piezometric wells in the area, indicating that achieving more compact piezometric wells information can potentially remove this offset. On the other hand, these offsets are likely attributable to the differing thickness and hydraulic properties of the aquifers and aquitards underlying these areas. 4 Discussion In the past two decades, land subsidence and ground failure in the Yazd-Ardakan flat area has caused significant damage to existing roads and structures. This phenomenon is a direct result of the aforementioned subsidence. In fact, earth rupture has affected numerous areas in the region, including Meybod, Ardakan, Rastagh, Ashkezar, Zarch and Yazd city. Earth fissure in the region is perpendicular to the northwest–southeast direction of the subsidence. Geologically speaking, the subsidence area is located in a clay flat and cultivated land surrounded by gypsiferous marl bed to the right and young traces and gravel fencing to the left. An assessment of the region’s fault structure reveals that there is evidence of existing undetermined yet inferred faults in the gypsiferous marl bed to the right of the subsidence area. The total length of this fault is about 24 km, beginning from the east of Meybod and ending 6 km north of Zarch. Also, there is another undetermined fault to the west of Ashkezar, which surrounds the right hand of the end of the subsidence zone. Nevertheless, further research needs to be conducted to derive the effect of geological structures on the aquifer system in Yazd province. The storage coefficient of an aquifer system is a parameter that contains the responses of the aquifer and fine-grained interbeds to variations in the hydraulic head. It is a 9
critical hydraulic parameter that strongly influences the non-steady flow of groundwater and is important to ground water resource evaluations. For investigation of the skeletal storage coefficient of the aquifer system, the water level variations ( h ) supposed to represent the effective stress and the InSAR deformation timeseries ( b ) showing the aquifer compaction are used. The inverse slope of the best fitting line to the stress-strain diagram is a rough estimate of the skeletal storage coefficient of the aquifer system (Hoffmann et al. 2001). b Here, the mean rate of subsidence ( ) acquired by InSAR indicates the compaction t of the aquifer system, while the mean rate of the water level variations in piezometric h wells at the same time ( ) denotes the effective stress. The storage coefficient of t the existing piezometric wells ( S ) can be computed as follows: b )4) b S t h h t Fig. 8 illustrates the computed storage coefficients for the existing piezometric wells in the region (green dots) depicted by contour lines. Significant negative values could be observed in the north and west of Yazd city. Here, despite a decline in water levels, the InSAR deformation map shows land surface uplift. This might be the result of low coherence between the Envisat images of these regions, reducing InSAR accuracy in these areas. More accurate methods, such as PSInSAR, are advised for further research in this area. On the other hand, this phenomenon be the result of seasonal effects in this region, as discussed in section 3, which have not been reflected well in comparison with the mean rate of earth surface deformation and mean rate of water level variations. The results show that storage coefficient values in the aquitard section are in the range of 1.19 10 4 in Sharifabad (west of Ardakan) to 0.114 in Roknabad (south of Meybod). The estimated storage coefficient of the maximum subsidence area is 0.1 .
10
Figure 8. Contour map of the storage coefficient of the Yazd-Ardakan aquifer system. Green dots are the piezometric wells.
Comparing the storage coefficient contour map with a geological map of the area showed that the general pattern and orientation of the storage coefficient map (figure 8) is compatible with the different geological structure in the area. The subsidence rate in the Yazd-Ardakan desert, as computed by InSAR in this study, is significantly higher than other subsidence areas, such as Antelope Valley, California (Galloway et al. 1998) and Las Vegas Valley, Nevada (Bell et al., 2008). Also, the storage parameter computed in this study is nearly ten times greater than those of Las Vegas and Antelope Valley (Galloway et al. 1998, Bell et al., 2008). The other critical subsidence area in Iran the storage parameter of which has previously been estimated by InSAR (Dehghani et al. 2009) is Neyshabour. The maximum storage parameter of Neyshabour is about three times greater than that of our case study, whereas the maximum rate of water level decline is about half of the YazdArdakan ones. Thus, it can be said that the higher storage parameter of Neyshabour explains the higher subsidence rate there. Unless improved groundwater management plans are implemented, it is anticipated that continued pressure on the groundwater resources in Yazd province, with rates of extraction far beyond those of natural recharge, might cause serious damages to urban, agricultural and industrial areas in the region. The maximum annual subsidence rate of 120 mm/yr inferred based on InSAR measurements in this study is one of the highest rates of deformation ever recorded for groundwater basins subject to groundwater development. Land subsidence and earth fissures resulting from overexploitation of water are serious geologic hazards the impacts of which will increase as more water is withdrawn from aquifers in the Yazd-Ardakan region without adequate replenishment. 5 Conclusions Differential interferometric SAR techniques were employed to investigate the temporal and spatial behavior of surface deformations in Yazd province, central Iran. A good temporal sampling rate of the monitoring done in the area was obtained, 11
which also enjoyed a high degree of spatial coverage over the area. The correspondence between precise leveling measurements of the subsidence and InSAR results confirm the accuracy of the InSAR deformation map of the region. Based on the subsidence map generated, 3 distinct patterns of fringes were found. A NW-SE elliptical-shaped bowl shape was apparent, with a maximum rate of 120 mm/yr between Meybod and Zarch, 100 mm/yr to the west of Ardakan, and 65 mm/yr to the south of Yazd city, the capital of the province for the 2003–2006 time period. These regions are located on cultivated land, clay flat and sand dunes. Time series analysis revealed that the first two areas suffer from subsidence with linear behavior, whereas the third has a periodic behavior. This means that water over-drafting without compensation in Meybod- Zarch and Ardakan are evident, whereas there is fluctuated land subsidence in Yazd city due to a lower water extraction rate that can partially be compensated for in the rainy season. An assessment of the groundwater level in the region revealed that the over-extraction of groundwater for agricultural and industrial purposes is the common reason of subsidence in all 3 areas. A water level decline of about 15 m has been recorded in the region during the previous three decades. This study evidenced the importance of taking the thickness and hydraulic properties of aquifer systems in water pumping managing into consideration, as the areas where the highest amounts of pumping were recorded did not necessarily suffer from maximum subsidence. The results of this study made it possible to divide the damaged area to high-risk regions (the area between Meybod and Zarch and Ardakan city where water over-drafting has caused significant linear subsidence and requires serious attention so as to prevent further damage) and a middle-risk area (Yazd city, where water over-drafting is partially compensated for in the rainy season. Measures need to be taken for the prevention of subsidence growth in Yazd city. Stress changes in the Yazd-Ardakan aquifer system were determined by means of piezometric well observations, while InSAR measurements were used to estimate the storage coefficient, an important parameter for the management of groundwater resources, in the aquifer system. The generated storage coefficient contour map simulated the pattern and orientation of the geological structure of the area. The spatial-temporal displacement field produced over the course of this study can be a useful resource for modeling the subsidence and determining the permitted amount of water extraction in the area, thus empowering water-pumping management in the area for the prevention of future damage. Acknowledgement SAR data have been made available by ESA, though the CAT-1 project proposal. The authors wish to thank National cartographic center (NCC) of Iran for funding this study. References Amelung, F., Galloway, D. L., Bell, J. W., Zebker, H. A. & Laczniak, R. J., 1999. Sensing the ups and downs of Las Vegas–InSAR reveals structural control of land subsidence and aquifer-system deformation: Geology, 27,483–486. Amighpey, M., Arabi, S., Talebi, A. & Djamour, 2006. Elevation changes of the precise leveling tracks in the Iran leveling network, Scientific report published in National Cartographic Center (NCC) of Iran, Tehran, Iran.
12
Amighpey, M., Mousavi. Z., Nankali, H., Arabi, S., Sedighi, M., Hosseini, S., 2008. Studying subsidence in Iran with leveling and permanent GPS observations, Geophysical Research Abstracts, Vol. 11, EGU2009-8189, EGU General Assembly 2008, Vienna, Austria. Bell, J. W., Amelung, F., Ferretti, A., Bianchi, M., Novali, F., 2008. Permanent scatterer InSAR reveals seasonal and long-term aquifer-system response to groundwater pumping and artificial recharge, Water Resour. Res., 44, W02407, doi:10.1029/2007WR006152. Berardino, P., Fornaro, G., Lanari, R., Sansosti, E., 2002. A new algorithm for surface deformation monitoring based on small baseline differential sar interferograms, IEEE Trans. Geosci. Remote Sens., vol. 40, no. 11,pp. 2375–2383. Bomford, G., 1971. Geodesy (third edition), Oxford University Press, England. Dehghani, M., M. J. Valadan Zoej, I. Entezam, S. Saatchi, 2009. InSAR Monitoring of Progressive Land Subsidence in Neyshabour, Northeast Iran, Geophysical Journal International, 178 (1): 47-56. Ferretti, A., Prati, C., & Rocca, F., 2000. Non-linear subsidence rate estimation using permanent scatterers in differential SAR interferometry, IEEE Trans. Geosci. Remote Sens., 38, 5. Freeze, R.A. & Cherry, J.A., 1979. Groundwater, Englewood Cliffs, N.J.,PrenticeHall, c1979. UCD Phys Sci GB1003.2.F73. Funning, G. J., Parsons, B., Wright, T. J., Jackson, J. A., Fielding, E. J., 2005. Surface displacements and source parameters of the 2003 Bam (Iran) earthquake from Envisat advanced synthetic aperture radar imagery, J. Geophys. Res., 110, B09406, doi:10.1029/2004JB003338. Galloway, D.L., Hudnut, K.W., Ingebritsen, S.E., Phillips, S.P., Peltzer, G.,Rogez, F. & Rosen, P.A., 1998. Detection of aquifer system compaction and land subsidence using interferometric synthetic aperture radar, Antelope valley, Mojave Desert, alifornia, Water Resour. Res., 34, 2573–2585. Hoffmann, J., Zebker, H.A., Galloway, D.L., Amelung, F., 2001. Seasonal subsidence and rebound in Las Vegas Valley, Nevada, observed by synthetic aperture radar interferometry, Water Resour. Res., 37(6), 1551– 1566. Kampes, B. & Usai, S., 1999. Doris: The Delft Object-oriented Radar Interferometric software. 2nd International Symposium on operationalization of Remote Sensing, ITC, Enschede, the Netherlands. 4 pages (on CDROM). Lanari, R., Lundgren, P., Manzo, M. & Casu, F., 2004. Satellite radar interferometry time series analysis of surface deformation for Los Angeles, California, Geophys. Res. Lett., 31, L23613, doi:10.1029/2004GL021294. Mood, A., Graybill, F., Boes, D., 1974. Introduction to the Theory of Statistics. third 13
edition, McGraw-Hill, New York. Mora, O. R. Lanari, J. J. Mallorquí, P. Berardino, and E. Sansosti, 2002. A new algorithm for monitoring localized deformation phenomena based on small baseline differential SAR interferograms,” in Proc. IGARSS, Toronto, ON, Canada, pp. 1237– 1239. Moreau, F., Dauteuil, O., Bour, O. & Gavrilenko, P., 2006. GPS measurements of ground deformation induced by water level variations into a granitic aquifer (French Brittany), Terra Nova, 18, 50–54. Motagh, M., Djamour, Y., Walter, T. R., Wetzel, H. U., Zschau, J., Arabi, S., 2007. Land subsidence in Mashhad Valley, northeast Iran; results from InSAR, levelling and GPS, Geophys. J. Int., 168(2), 518–526. Motagh, M., Walter, T. R., Sharifi, M.A., Fielding, E., Schenk, A., Anderssohn, J., Zschau1, J., 2008. Land subsidence in Iran caused by widespread water reservoir overexploitation, Geophys. Res. Lett., VOL. 35, L16403. Schmidt, D.A. & B¨urgmann, R., 2003. Time dependent land uplift and subsidence in the Santa Clara valley, California, from a large InSAR data set, J. geophys. Res., 108, doi:10.1029/2002JB002267. Terzaghi, K., 1925. Principles of soil mechanics, IV, settlement and consolidation of clay, Eng. News Rec., 95(3), 874–878. Usai, S., 2003. A least squares database approach for SAR interferometric data, IEEE Trans. Geosci. Remote Sens., 41(4), 753−760. Zebker, H.A., Rosen, P.A., Goldstein, R.M., Gabriel, A. & Werner, C.L., 1994. On the derivation of coseismic displacement fields using differential radar interferometry: the Landers earthquake, J. geophys. Res., 99, 19 617–19 634.
Highlights We apply the small baseline algorithm to produce time series of deformatin in the Yazd-Ardakan plain. Mean rate of subsidence computed from repeated first order precise levelling network of Iran is accurately matched by the mean rate of InSAR time series. Peak amplitude of the subsidence is estimated about 120 mm/yr for the 2003– 2006 time period. Applying the subsidence map produced by InSAR and water level change contour, skeletal storage coefficient is estimated. The storage coefficient contour map generated in this study simulated the pattern and orientation of the geological structure of the area.
14
PII: DOI: Reference:
S2352-9385(16)30020-9 http://dx.doi.org/10.1016/j.rsase.2016.04.001 RSASE23
To appear in: Remote Sensing Applications: Society and Environment Received date: 10 August 2015 Revised date: 11 April 2016 Accepted date: 12 April 2016 Cite this article as: Masoome Amighpey and Siyavash Arabi, Studying land subsidence in Yazd province, Iran, by integration of InSAR and levelling measurements, Remote Sensing Applications: Society and Environment, http://dx.doi.org/10.1016/j.rsase.2016.04.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Studying land subsidence in Yazd province, Iran, by integration of InSAR and levelling measurements Masoome Amighpeya, Siyavash Arabib a
National Cartographic Center of Iran, Meraj Ave., Azadi Squ., Tehran, Iran, [email protected], Tel.: +98 21 66071001
b
National Cartographic Center of Iran, Meraj Ave., Azadi Squ., Tehran, Iran, [email protected]
Abstract The Yazd-Ardakan plain, located in central Iran, is an arid region with a low rainfall where over-extraction of groundwater for agriculture and industrial purposes has caused serious subsidence in the area. In this study, the small baseline subset (SBAS) algorithm was applied to perform an interferometric synthetic aperture radar (InSAR) time series analysis for the Yazd-Ardakan plain. The result is a space-time deformation product that can be exploited to view not only the smoothly varying long-term surface motion, but also its time-varying patterns. The mean rate of subsidence computed from the repeated first order precise leveling network of Iran was accurately matched by the mean rate of the InSAR time series. Time series analysis of our InSAR mapping suggests that subsidence occurs within a northwest– southeast elliptic-shaped bowl, with peak amplitude of about 120 mm/yr for the 2003–2006 time period. Analysis of piezometric records suggests that subsidence likely results from extensive over-drafting of the aquifer system in the region, which has caused about 16 m of water table decline since 1974. Comparing the subsidence map produced by InSAR and water level change contour allowed the estimation of skeletal storage coefficients, providing a basis for future work that could characterize the storage properties of an aquifer system with a high degree of spatial resolution. The storage coefficient contour map generated in this study simulated the pattern and orientation of the geological structure of the area. 1 Introduction Iran is a country with wide spatial and temporal limitations regarding water resources, especially in the central arid areas. This problem has become more conspicuous in recent decades due to economic and urban development. Economic and population growth have deteriorated available water resources in terms of both quantity and quality aspect. The significant amount of water obtained from groundwater resources is such that land subsidence has become an issue in some areas (Amighpey et al., 2006, Motagh et al., 2007, Amighpey et al., 2008, Motagh et al., 2008, Dehghani et al., 2009). Incidentally, many of these subsidence regions are located in the precise leveling routes of Iran’s National Cartographic Center (NCC) network, which were observed twice in the last four decades and then assessed to determine potential height changes. In fact, these assessments were the first to detect the rate and location of this subsidence (Amighpey et al., 2006). Yazd province, located in central Iran, is one of these areas in which assessments of repeated precise leveling observations showed an extreme degree of subsidence, especially in the Yazd-Ardakan flat area (Amighpey et 1
al., 2006). The Yazd-Ardakan desert, with a length of 60 km and width of 15 km, is located between Yazd and Ardakan cities. Over the past two decades, subsidence and rupture of the earth have led to significant damage in Yazd province. This phenomenon has caused roads to rupture, water systems to break and structures to crack. In fact, one of the main considerations of construction work in Yazd is the numerous instances of earth rupture between Yazd and Ardakan. This phenomenon is the direct result of subsidence. Being located beside the central mountains, far from sea, adjacent to the desert and in a low rainfall region, Yazd has a climate which mostly resembles dry desert climates. Little rain along with high water evaporation, relatively low dampness, considerable heat and great temperature changes are among the factors making this province one of the driest parts of Iran, yet water needs for consumption, agriculture and industry are met mainly by groundwater resources. This has brought about over-extraction of groundwater in recent decades, which is, in turn, the main reason of subsidence in this area. Time series analysis of satellite radar interferometry is a useful method for studying both spatial and temporal patterns of surface deformation and has been extensively used to study subsidence areas (Berardino et al., 2002; Usai, 2003; Schmidt and Bu¨rgmann, 2003). In the present study, the InSAR deformation time series, based on the small baseline subset (SBAS) algorithm (Berardino et al., 2002), was computed for the Yazd-Ardakan desert, using data from the Envisat satellite acquired from the European Space Agency (ESA). We found that the mean rate of InSAR time series analysis agreed well with re-leveling results later used to measure subsidence. The availability of such a joint temporal/spatial deformation data set, when combined with independent hydrologic and geologic assessments, can provide more information on the storage properties of the aquifer system. A better understanding of the areal variability of the aquifer system response to stress can improve the effectiveness of groundwater management schemes and the identification of zones with a high potential of fissure formation, which is valuable information for urban planning. In this study, a simple temporal model for the subsidence was obtained, one capable of also depicting its spatial extent. The spatial and temporal behavior of the subsidence was correlated to water-level fluctuations to characterize aquifer-system response patterns and to further refine measurements of the aquifer storage properties. Analyzing the temporal behavior of subsidence and groundwater levels in the region made it possible to distinguish high-risk areas, where high water over-drafting has caused linear subsidence, from moderate-risk ones, where water over-drafting is partially compensated for in the region’s rainy season. This could prove to be vital information for the management of water pumping in the area so as to prevent further damage. 2 Subsidence in the Yazd-Ardakan plain from precise leveling measurements Geodetic observations of local land subsidence for the Yazd-Ardakan region were first documented using precise leveling data collected by the National Cartographic Center of Iran (NCC) during 1993-2001. The pink color points in Fig. 1 (a) show the location of leveling benchmarks extracted from first-order precise leveling lines along the Naeen-Yazd Highway. A number of stations were destroyed between 1993-2001.
2
Figure1. (a) SRTM DEM of Yazd subsidence area. Green points denote the first order leveling routes and pink points denote parts of these routes that show subsidence. (b) Results of first-order precise leveling surveys in the Yazd-Ardakan region between 1993-2001. Error bars represent the standard deviations in cumulative random error from the reference benchmark (BM15).
Analysis was done on the development of vertical deformation along the Naeen-Yazd Highway between the 1993 and 2001 surveys relative to a datum at BM15 and BM17, near Nogonbad and outside of the subsidence area. Figure 1(b) depicts the height changes along these routes. For the 8-year period spanning1993–2001, the maximum cumulative degree of subsidence recorded by a benchmark is as high as 95.00 cm 0.92cm at a point located between Shamsi and Hojjatabad (BM62). Ardakan (BM50) and Yazd city (BM83) also subsided by 47.00 cm 0.79 cm and 22.00 cm 1.1cm during the same 8-year time span, respectively. Precise leveling accuracy was computed by means of the Vignal formula (Bomford, 1971). Data pertaining to temperature and air pressure are not available for this route, which could have been used to remove the refraction effect, but the fact that there was no correlation between vertical displacements and topography makes it evident that elevation-dependent systematic error has had no impact on the leveling data. InSAR time series analysis is a better resource for studying subsidence versus two existing epochs over a long interval of leveling observation. This temporal/spatial map can be particularly useful in providing good input for modeling subsidence in the consolidation theory. In the following, the results of an investigation of the spatial and temporal behavior of Yazd-Ardakan subsidence patterns using surface displacement maps generated by InSAR are presented. 3 InSAR measurements of subsidence and water level change assessments The concept of combining InSAR information from a large number of SAR Images, which can facilitate the solution of a deformation time series, is a relatively new one (Ferretti et al., 2000; Berardino et al., 2002; Usai, 2003; Mora et al., 2002; Lanari et al., 2004). Using a multi-temporal technique not only allows the estimation of discrete deformation events, but also provides an avenue for a temporal study provided that the whole time series is available. To measure ground subsidence in the Yazd-Ardakan region, we selected 10 Envisat SAR images (see Table 1) in descending orbits. Employing the so-called two-pass technique (Zebker et al. 1994), a set of 43 differential interferograms was generated covering a time period between March 26, 2003 and December 20, 2006. Interferometric processing was done using the public domain InSAR processor 3
DORIS (Kampes & Usai 1999). In this method, the interferograms are corrected for the phase signature due to orbital separation using precise DEOS satellite orbits, and for the topography using a 3-arcsecond SRTM digital elevation model (more information available at http://srtm.usgs.gov). The remaining orbital error was then removed by fitting a plane in the far field area where the leveling data showed no deformation and subtracting this plane from the interferograms. Areas with low correlation were filtered out from the interferograms. The availability of both time and space data makes it possible to effectively identify and remove atmospheric artifacts in the results via stacking and the use of the time dependent characteristic of said atmospheric artifacts. Using multiple independent interferograms helps find atmospheric artifacts where fringes existed in special interferograms only, as this effect depends on the acquisition time. It is also worth noting that precise leveling was used to test the accuracy of the computed deformation map. In the end, only 19 interferograms which enjoyed a high correlation coefficient were chosen for time series analysis (Table 2). Figure 2 shows the spatial vs. temporal baseline of these interferograms. Figure 3(a) shows the interferogram generated from images acquired on September 17, 2003 and September 1, 2004. As the interferograms represent the deformation in the time interval between the acquisition dates of master and slave images, least square adjustment was used to invert deformation in different time intervals (interferograms) to a time series of deformation indicating displacement at each date of acquisition. If D is the displacement vector of the acquisition date of the images, which is unknown, and I is the displacement vector between the slave and master images acquisition interval (interferograms), we have the following equation: AD = I (1) Where A is the coefficient matrix. Applying least square adjustment, we have: D ( AT PA) 1 AT PI (2) CD ( AT PA) 1
(3)
Where P is the inverse of the observation variance matrix, while C D is the variance matrix of the unknowns. Table1. Details of the used images in this study Sensor Track Pass Swath Date Orbit ASAR/IM 20 Descending I2 2003-03-26 5585 ASAR/IM 20 Descending I2 2003-09-17 8090 ASAR/IM 20 Descending I2 2003-11-26 9092 ASAR/IM 20 Descending I2 2004-06-23 12098 ASAR/IM 20 Descending I2 2004-09-01 13100 ASAR/IM 20 Descending I2 2004-12-15 14603 ASAR/IM 20 Descending I2 2005-03-30 16106 ASAR/IM 20 Descending I2 2005-10-26 19112 ASAR/IM 20 Descending I2 2006-02-08 20615 ASAR/IM 20 Descending I2 2006-12-20 25124 Table2. Details of the used interferogram in this study Interferogram Interferogram name Temporal baseline Baseline length (m) number (day) 1 20030326_20031126 245 118 2 20030326_20030917 175 434 3 20030326_20040623 455 153 4 20030326_20040901 525 403 Mission Envisat Envisat Envisat Envisat Envisat Envisat Envisat Envisat Envisat Envisat
4
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
20030326_20041215 20030917_20040901 20030917_20051026 20031126_20040623 20031126_20040901 20031126_20041215 20031126_20061220 20040623_20041215 20040623_20061220 20040901_20051026 20040901_20061220 20041215_20060208 20041215_20061220 20050330_20060208 20051026_20061220
630 350 770 210 280 385 1120 175 910 420 840 420 735 315 420
160 32 127 35 521 42 271 7 306 159 250 427 314 147 409
Figure2. Spatial vs. temporal baseline of the 19 interferograms used in this study.
To compute the InSAR time-series, 19 interferograms were inverted using the least square adjustment method, weighting each data set by the inverse of its far field variance (Funning et al., 2005). 9 deformation maps were obtained into March 26, 2003, the image acquisition date prior to which no deformation is supposed (the first image acquisition date used as the reference epoch). Figure 3(b) shows the mean rate of the subsidence in this area obtained from stacked interferograms. Looking at figure 3(a), the interferogram belonging to the period between September 17, 2003 and September 1, 2004, 3 fringe patterns are evident. The subsidence pattern can be divided into 3 parts: the largest rate occurring in an extensive spatial distribution between Meybod and Zarch shaped like a bowl and directed northwest– southeast along the axis of Yazd-Ardakan road, the second to the west of Ardakan, and the third in the south of Yazd city, capital of Yazd province.
5
(a)
(b)
Figure3. (a) Interferogram generated from images acquired on September 17, 2003 and September 1, 2004, (b) Mean rate of subsidence obtained from stacked interferograms (m/yr). Contour lines represent the water level change rate in m/yr, derived from piezometric wells (black triangles) in the 2003-2007 time interval.
The highest rate of subsidence in the Meybod-Zarch region can be found between Shamsi and Hojjatabad (120 mm/yr). Figure 4(a) shows the time series of the maximum subsidence area. As is apparent, the subsidence demonstrates a linear behavior with the passage of time. R-Square and RMSE (Root Mean Squared Error) statistics were used to evaluate the degree of linear fit (cf. Mood et al., 1974). Rsquare is the square of the correlation between the response values and the predicted response values and is used to measure how successful the fit is in explaining variation the data. . It is also called the square of the multiple correlation coefficients and the coefficient of multiple determinations. RMSE is known as the fit standard error. The R-square of the linear fitness in the time series is 0.97, meaning that the fit explains 97% of the total variation in the data. Also of note is that the RMSE of linear fit is 10 mm. Figure 4(a) also depicts the fitted regression line.
6
(a)
(b)
(c) Figure 4. InSAR time series of subsidence in: (a) maximum subsidence area between Shamsi and Hojjatabad, (b) Ardakan, (c) Yazd. Squares are vertical displacement at the image acquisition time. Error bars represent the estimated accuracy of vertical displacement. The blue line is the fitted regression line of vertical displacement.
As evidenced by the data, the highest rate of subsidence occurred in agricultural areas (where groundwater extraction is common) and caused the land to rupture and buildings to fracture. Figure 5 illustrates an instance of the impact of subsidence in the maximum rate area (Hojjatabad). The earth fissure orientation is perpendicular to the northwest–southeast direction of the subsidence.
Figure 5. Field observation of the subsidence effect on the maximum subsidence area in Hojjatabad.
Another bowl-shaped pattern of subsidence was detected in Ardakan. Figure 4(b) depicts the relevant time series and linear fitted model. The assessment of fitting linear regression reveals linear subsidence with a maximum rate of 100 mm. The Rsquare value of the linear fitness to Ardakan time series is 0.98, meaning that the fit explains 98% of the total variation in the data about the average. The RMSE of linear fit is 7 mm. This subsidence area is located in an agricultural land and industrial region. Yazd, the capital of the Yazd province, showed distinct fringes with a lower rate compared with Ardakan and Shamsi. These fringes appeared near the railway in Yazd city, around an industrial estate. Results of the assessment of the degree of fit yielded an R-square value equal to 0.5747 and RMSE of 20 mm, rejecting a simple linear 7
pattern and emphasizing seasonal effects )figure 4 (c)). Figure 6 shows the time series of deformation in the city. Yazd’s time series oscillations show instances of strong seasonal uplift in the area in the spring of each year due to the effects of annual artificial recharge (e.g. figure 6 (f) and 6 (h)). The highest subsidence rate is about 65 mm/yr.
Figure 6. Time series of elevation deformation in Yazd city.
To assess the accuracy of InSAR measurements, the subsidence rate acquired from leveling observations was projected on the satellite line of sight and then compared with the average velocity from the InSAR time series. Although the two leveling observations belong to 1993 and 2001, whereas the radar images were acquired in the 2003-2007 time interval, (in different time periods), comparing the mean subsidence rate obtained from these two distinct sources yields a good degree of correspondence in leveling benchmarks (figure 7). The existing RMS between them is 15 mm. This correspondence, together with the accuracy of the InSAR results, proves the existence of a constant subsidence rate over these years. The existence of such a constant subsidence rate indicates a lack of water pumping management in recent years.
Figure 7. Comparison between mean subsidence rate computed from leveling projected on the satellite line of sight and InSAR in the path of leveling.
8
To better understand the cause of land subsidence in Yazd-Ardakan plain, variations in groundwater levels in the piezometric wells in the region were studied. Groundwater level measurements have been regularly made in the region every month since 1973. The contour lines in Figure 3(b) show the behavior of the mean rate of water level decline in this zone and around the subsiding area. It is evident that a continuous lowering of the hydraulic head below the ground level has been occurring in the region over recent decades. For the time interval of recorded piezometric measurements from 1973 to 2007, the mean water table in this area declined about 15 m as a result of intensive groundwater pumping and deficient natural recharge. In principle, the over-extraction of groundwater leads to a reduction in water pressure within the aquifer system which is reflected by an increase in the effective stress and compaction of sediments (Terzaghi 1925). The relationship between these two processes depends on many factors, including geotechnical parameters such as sediment porosity, permeability and compressibility as well as the spatial distribution of the thickness of compressible layers within the aquifer system and the history and amounts of pressure drawdown caused by groundwater decline (Freeze & Cherry 1979). Some external processes such as tide, recharge, or regional tectonics might also affect the deformation of an aquifer system (Amelung et al. 1999; Moreau et al. 2006). As evidenced by Figure 3(b), although there is a clear correlation between the rate of water level reduction and the subsidence, the areas with the highest degree of aquifersystem compaction are offset from the zones of maximum water-level decline in all 3 subsidence areas: the Meybod-Zarch region, Ardakan and Yazd. Such offset between the areas of greatest subsidence and most significant water-level decline has previously been reported in the Las Vegas metropolitan area (Bell et al. 2008). This might be due to the higher spatial resolution of the earth surface deformation map produced by InSAR in contrast with the water level changes contours provided by the sparse piezometric wells in the area, indicating that achieving more compact piezometric wells information can potentially remove this offset. On the other hand, these offsets are likely attributable to the differing thickness and hydraulic properties of the aquifers and aquitards underlying these areas. 4 Discussion In the past two decades, land subsidence and ground failure in the Yazd-Ardakan flat area has caused significant damage to existing roads and structures. This phenomenon is a direct result of the aforementioned subsidence. In fact, earth rupture has affected numerous areas in the region, including Meybod, Ardakan, Rastagh, Ashkezar, Zarch and Yazd city. Earth fissure in the region is perpendicular to the northwest–southeast direction of the subsidence. Geologically speaking, the subsidence area is located in a clay flat and cultivated land surrounded by gypsiferous marl bed to the right and young traces and gravel fencing to the left. An assessment of the region’s fault structure reveals that there is evidence of existing undetermined yet inferred faults in the gypsiferous marl bed to the right of the subsidence area. The total length of this fault is about 24 km, beginning from the east of Meybod and ending 6 km north of Zarch. Also, there is another undetermined fault to the west of Ashkezar, which surrounds the right hand of the end of the subsidence zone. Nevertheless, further research needs to be conducted to derive the effect of geological structures on the aquifer system in Yazd province. The storage coefficient of an aquifer system is a parameter that contains the responses of the aquifer and fine-grained interbeds to variations in the hydraulic head. It is a 9
critical hydraulic parameter that strongly influences the non-steady flow of groundwater and is important to ground water resource evaluations. For investigation of the skeletal storage coefficient of the aquifer system, the water level variations ( h ) supposed to represent the effective stress and the InSAR deformation timeseries ( b ) showing the aquifer compaction are used. The inverse slope of the best fitting line to the stress-strain diagram is a rough estimate of the skeletal storage coefficient of the aquifer system (Hoffmann et al. 2001). b Here, the mean rate of subsidence ( ) acquired by InSAR indicates the compaction t of the aquifer system, while the mean rate of the water level variations in piezometric h wells at the same time ( ) denotes the effective stress. The storage coefficient of t the existing piezometric wells ( S ) can be computed as follows: b )4) b S t h h t Fig. 8 illustrates the computed storage coefficients for the existing piezometric wells in the region (green dots) depicted by contour lines. Significant negative values could be observed in the north and west of Yazd city. Here, despite a decline in water levels, the InSAR deformation map shows land surface uplift. This might be the result of low coherence between the Envisat images of these regions, reducing InSAR accuracy in these areas. More accurate methods, such as PSInSAR, are advised for further research in this area. On the other hand, this phenomenon be the result of seasonal effects in this region, as discussed in section 3, which have not been reflected well in comparison with the mean rate of earth surface deformation and mean rate of water level variations. The results show that storage coefficient values in the aquitard section are in the range of 1.19 10 4 in Sharifabad (west of Ardakan) to 0.114 in Roknabad (south of Meybod). The estimated storage coefficient of the maximum subsidence area is 0.1 .
10
Figure 8. Contour map of the storage coefficient of the Yazd-Ardakan aquifer system. Green dots are the piezometric wells.
Comparing the storage coefficient contour map with a geological map of the area showed that the general pattern and orientation of the storage coefficient map (figure 8) is compatible with the different geological structure in the area. The subsidence rate in the Yazd-Ardakan desert, as computed by InSAR in this study, is significantly higher than other subsidence areas, such as Antelope Valley, California (Galloway et al. 1998) and Las Vegas Valley, Nevada (Bell et al., 2008). Also, the storage parameter computed in this study is nearly ten times greater than those of Las Vegas and Antelope Valley (Galloway et al. 1998, Bell et al., 2008). The other critical subsidence area in Iran the storage parameter of which has previously been estimated by InSAR (Dehghani et al. 2009) is Neyshabour. The maximum storage parameter of Neyshabour is about three times greater than that of our case study, whereas the maximum rate of water level decline is about half of the YazdArdakan ones. Thus, it can be said that the higher storage parameter of Neyshabour explains the higher subsidence rate there. Unless improved groundwater management plans are implemented, it is anticipated that continued pressure on the groundwater resources in Yazd province, with rates of extraction far beyond those of natural recharge, might cause serious damages to urban, agricultural and industrial areas in the region. The maximum annual subsidence rate of 120 mm/yr inferred based on InSAR measurements in this study is one of the highest rates of deformation ever recorded for groundwater basins subject to groundwater development. Land subsidence and earth fissures resulting from overexploitation of water are serious geologic hazards the impacts of which will increase as more water is withdrawn from aquifers in the Yazd-Ardakan region without adequate replenishment. 5 Conclusions Differential interferometric SAR techniques were employed to investigate the temporal and spatial behavior of surface deformations in Yazd province, central Iran. A good temporal sampling rate of the monitoring done in the area was obtained, 11
which also enjoyed a high degree of spatial coverage over the area. The correspondence between precise leveling measurements of the subsidence and InSAR results confirm the accuracy of the InSAR deformation map of the region. Based on the subsidence map generated, 3 distinct patterns of fringes were found. A NW-SE elliptical-shaped bowl shape was apparent, with a maximum rate of 120 mm/yr between Meybod and Zarch, 100 mm/yr to the west of Ardakan, and 65 mm/yr to the south of Yazd city, the capital of the province for the 2003–2006 time period. These regions are located on cultivated land, clay flat and sand dunes. Time series analysis revealed that the first two areas suffer from subsidence with linear behavior, whereas the third has a periodic behavior. This means that water over-drafting without compensation in Meybod- Zarch and Ardakan are evident, whereas there is fluctuated land subsidence in Yazd city due to a lower water extraction rate that can partially be compensated for in the rainy season. An assessment of the groundwater level in the region revealed that the over-extraction of groundwater for agricultural and industrial purposes is the common reason of subsidence in all 3 areas. A water level decline of about 15 m has been recorded in the region during the previous three decades. This study evidenced the importance of taking the thickness and hydraulic properties of aquifer systems in water pumping managing into consideration, as the areas where the highest amounts of pumping were recorded did not necessarily suffer from maximum subsidence. The results of this study made it possible to divide the damaged area to high-risk regions (the area between Meybod and Zarch and Ardakan city where water over-drafting has caused significant linear subsidence and requires serious attention so as to prevent further damage) and a middle-risk area (Yazd city, where water over-drafting is partially compensated for in the rainy season. Measures need to be taken for the prevention of subsidence growth in Yazd city. Stress changes in the Yazd-Ardakan aquifer system were determined by means of piezometric well observations, while InSAR measurements were used to estimate the storage coefficient, an important parameter for the management of groundwater resources, in the aquifer system. The generated storage coefficient contour map simulated the pattern and orientation of the geological structure of the area. The spatial-temporal displacement field produced over the course of this study can be a useful resource for modeling the subsidence and determining the permitted amount of water extraction in the area, thus empowering water-pumping management in the area for the prevention of future damage. Acknowledgement SAR data have been made available by ESA, though the CAT-1 project proposal. The authors wish to thank National cartographic center (NCC) of Iran for funding this study. References Amelung, F., Galloway, D. L., Bell, J. W., Zebker, H. A. & Laczniak, R. J., 1999. Sensing the ups and downs of Las Vegas–InSAR reveals structural control of land subsidence and aquifer-system deformation: Geology, 27,483–486. Amighpey, M., Arabi, S., Talebi, A. & Djamour, 2006. Elevation changes of the precise leveling tracks in the Iran leveling network, Scientific report published in National Cartographic Center (NCC) of Iran, Tehran, Iran.
12
Amighpey, M., Mousavi. Z., Nankali, H., Arabi, S., Sedighi, M., Hosseini, S., 2008. Studying subsidence in Iran with leveling and permanent GPS observations, Geophysical Research Abstracts, Vol. 11, EGU2009-8189, EGU General Assembly 2008, Vienna, Austria. Bell, J. W., Amelung, F., Ferretti, A., Bianchi, M., Novali, F., 2008. Permanent scatterer InSAR reveals seasonal and long-term aquifer-system response to groundwater pumping and artificial recharge, Water Resour. Res., 44, W02407, doi:10.1029/2007WR006152. Berardino, P., Fornaro, G., Lanari, R., Sansosti, E., 2002. A new algorithm for surface deformation monitoring based on small baseline differential sar interferograms, IEEE Trans. Geosci. Remote Sens., vol. 40, no. 11,pp. 2375–2383. Bomford, G., 1971. Geodesy (third edition), Oxford University Press, England. Dehghani, M., M. J. Valadan Zoej, I. Entezam, S. Saatchi, 2009. InSAR Monitoring of Progressive Land Subsidence in Neyshabour, Northeast Iran, Geophysical Journal International, 178 (1): 47-56. Ferretti, A., Prati, C., & Rocca, F., 2000. Non-linear subsidence rate estimation using permanent scatterers in differential SAR interferometry, IEEE Trans. Geosci. Remote Sens., 38, 5. Freeze, R.A. & Cherry, J.A., 1979. Groundwater, Englewood Cliffs, N.J.,PrenticeHall, c1979. UCD Phys Sci GB1003.2.F73. Funning, G. J., Parsons, B., Wright, T. J., Jackson, J. A., Fielding, E. J., 2005. Surface displacements and source parameters of the 2003 Bam (Iran) earthquake from Envisat advanced synthetic aperture radar imagery, J. Geophys. Res., 110, B09406, doi:10.1029/2004JB003338. Galloway, D.L., Hudnut, K.W., Ingebritsen, S.E., Phillips, S.P., Peltzer, G.,Rogez, F. & Rosen, P.A., 1998. Detection of aquifer system compaction and land subsidence using interferometric synthetic aperture radar, Antelope valley, Mojave Desert, alifornia, Water Resour. Res., 34, 2573–2585. Hoffmann, J., Zebker, H.A., Galloway, D.L., Amelung, F., 2001. Seasonal subsidence and rebound in Las Vegas Valley, Nevada, observed by synthetic aperture radar interferometry, Water Resour. Res., 37(6), 1551– 1566. Kampes, B. & Usai, S., 1999. Doris: The Delft Object-oriented Radar Interferometric software. 2nd International Symposium on operationalization of Remote Sensing, ITC, Enschede, the Netherlands. 4 pages (on CDROM). Lanari, R., Lundgren, P., Manzo, M. & Casu, F., 2004. Satellite radar interferometry time series analysis of surface deformation for Los Angeles, California, Geophys. Res. Lett., 31, L23613, doi:10.1029/2004GL021294. Mood, A., Graybill, F., Boes, D., 1974. Introduction to the Theory of Statistics. third 13
edition, McGraw-Hill, New York. Mora, O. R. Lanari, J. J. Mallorquí, P. Berardino, and E. Sansosti, 2002. A new algorithm for monitoring localized deformation phenomena based on small baseline differential SAR interferograms,” in Proc. IGARSS, Toronto, ON, Canada, pp. 1237– 1239. Moreau, F., Dauteuil, O., Bour, O. & Gavrilenko, P., 2006. GPS measurements of ground deformation induced by water level variations into a granitic aquifer (French Brittany), Terra Nova, 18, 50–54. Motagh, M., Djamour, Y., Walter, T. R., Wetzel, H. U., Zschau, J., Arabi, S., 2007. Land subsidence in Mashhad Valley, northeast Iran; results from InSAR, levelling and GPS, Geophys. J. Int., 168(2), 518–526. Motagh, M., Walter, T. R., Sharifi, M.A., Fielding, E., Schenk, A., Anderssohn, J., Zschau1, J., 2008. Land subsidence in Iran caused by widespread water reservoir overexploitation, Geophys. Res. Lett., VOL. 35, L16403. Schmidt, D.A. & B¨urgmann, R., 2003. Time dependent land uplift and subsidence in the Santa Clara valley, California, from a large InSAR data set, J. geophys. Res., 108, doi:10.1029/2002JB002267. Terzaghi, K., 1925. Principles of soil mechanics, IV, settlement and consolidation of clay, Eng. News Rec., 95(3), 874–878. Usai, S., 2003. A least squares database approach for SAR interferometric data, IEEE Trans. Geosci. Remote Sens., 41(4), 753−760. Zebker, H.A., Rosen, P.A., Goldstein, R.M., Gabriel, A. & Werner, C.L., 1994. On the derivation of coseismic displacement fields using differential radar interferometry: the Landers earthquake, J. geophys. Res., 99, 19 617–19 634.
Highlights We apply the small baseline algorithm to produce time series of deformatin in the Yazd-Ardakan plain. Mean rate of subsidence computed from repeated first order precise levelling network of Iran is accurately matched by the mean rate of InSAR time series. Peak amplitude of the subsidence is estimated about 120 mm/yr for the 2003– 2006 time period. Applying the subsidence map produced by InSAR and water level change contour, skeletal storage coefficient is estimated. The storage coefficient contour map generated in this study simulated the pattern and orientation of the geological structure of the area.
14