TIME SERIES ANALYSIS OF INSAR DATA TO STUDY LAND SUBSIDENCE INDUCED BY GROUNDWATER LEVEL DECLINE IN SIRJAN PLAIN

Sirjan located in Southwest of Kerman City, Iran, is a city with average annual rainfall of 132mm. Over-exploitation of groundwater used for irrigation of pistachios gardens has caused serious land subsidence in Sirjan Plain. In this research we have used the Interferometric Synthetic Aperture Radar technique to estimate the subsidence rate in Sirjan plain. 18 Interferograms extracted from 12 ENVISAT ASAR images spanning between April 2004 and September 2010, have been studied. The SBAS algorithm was then applied in order to estimate the mean subsidence velocity map. Maximum subsidence rate is estimated as 28cm/yr. Furthermore, groundwater level fluctuations in the study area has been investigated in the piezometer wells located in the study area. Comparing between the results obtained from the interferometry and groundwater level fluctuation maps, shows a strong correlation between head decline in groundwater and land subsidence.


INTRODUCTION
Aquifers are one of the most important water sources for agriculture, industry and drinking.Excessive extraction of groundwater is followed by population growth and industrial activities' development, which will result in subsidence in ground surface (Gumilar, 2015).In fact, subsidence is downward movement of ground surface either gradually or suddenly, both can be result of natural causes or human activities.Groundwater level decline, degradation of subsurface structure strength and karst destruction, are the most general reasons for subsidence occurrence in Iran.Subsidence can disturb the stability of human made structures including historical buildings.The important point is to study and analyse accurately the subsidence, especially in potential places, to inhibit the consequences such as change in groundwater quality, increase in energy consumed for water extraction from aquifer, vulnerability of plain against drought, land subsidence, ecosystem extinction, and gardens drying out.Subsidence phenomena can be studied by hydrogeology methods such as groundwater level, extensometer, and piezometer observation and also by geodetic methods like mapping points using GPS and radar interferometry technique (Abidin, 2008).Radar interferometry technique is an effective and accurate method with high spatial resolution which is used for ground surface phenomena aimed to measure any kind of movement or deformation on ground surface (Kumar, 2011).Considering studies accomplished so far, it is concluded that radar interferometry technique is a reliable method with desired accuracy for estimating subsidence.On the other hand, interferometry failure due to temporal decorrelation can be reduced by Small Baseline Subset (SBAS) InSAR time series analysis method during a long time period (Zhou, 2013).Many researches have been carried out in the field of investigating subsidence using the potential of satellite images and InSAR time series analysis method.Dehghani et al. estimated the subsidence rate up to 19 cm per year for Neyshabour city employing radar interferometry technique using ENVISAT ASAR images (Dehghani, 2009).Khakim et al. used interferograms calculated from 21 images in ascending mode for the time period of Jan 1 st 2007 to Mar 3 rd 2011, and estimated subsidence rate about 45 cm for Bandung Basin, Indonesia (Khakim, 2014) In this paper, interferometry technique is used to study the subsidence of Sirjan Plain.12 images of ENVISAT ASAR as well as 53 piezometric wells spanning between Apr 2004 and Sep 2010 were applied to interpret the subsidence signal.In the next section, InSAR method to obtain the interferograms and SBAS time series algorithm used in this study will be introduced.In section 3, the results and discussion will be presented.Concluding remarks are presented in the last section.

METHODOLOGY
InSAR time series analysis of a significant number of interferograms was to study the spatial and temporal behaviour of the subsidence in Sirjan.Groundwater level fluctuations at piezometric wells were then compared to the results obtained from InSAR.Basic structure of the methodology containing the main processing steps is depicted in Figure 1.

Interferograms generation
12 ENVISAT ASAR radar images of track 260 of ENVISAT acquired from Apr 1 st 2004 to Sep 28 th 2010 are used to generate 18 differential interferograms.The acquisition dates are demonstrated in table 1.In order to mitigate the sptial and temporal decorrelation effect, the interferograms are characterized by small spatial and temporal baselines.The interferograms information are presented in table 2. The interferometric phase which is the phase difference between master and slave images is calculated as: Where decorrelation (Zhou, 2013).In order to estimate the deformation signal, all other components have to be estimated and subtracted from the interferometric phase.
On the other hand, the phase are modulated between   and   , called wrapped phase.Phase unwrapping describes the subsequent step in which the integer phase cycles associated with each relative wrapped phase measurement are estimated.

Time series analysis
Temporal decorrelation, i.e. changes in backscattering behaviour in time, is the most important limiting factor in InSAR technology.Interferograms with large temporal baselines are mainly subject to decorrelation [11].Therefore, it is preferred to generate interferograms characterised by small temporal baselines.Time series analysis using a significant number of interferograms allows for studying the temporal behaviour of the subsidence over large period of time.In order to generate the subsidence time series, the deformation at each acquisition date is estimated.Least squares solution will be applied for interferograms inversion.Suppose that: ' ,, 12 is a vector of differential interferometric phase, , 1, 2, i i     .associatedwith the temporal baseline i .Each interferogram is considered to be the phase difference between master and slave images as follows: is the phase of image j and n is the number of images.Observations equations used in the least squares approach are written as: The first image is set to be a known parameter as the time reference ( 0 1 ) Eq. ( 4) can be written as: AX L  (5) X is a vector containing the unknown values, A is the design matrix and L is the interferometric observations in least-square problem (Dehghani et al., 2009).To decrease residual orbital errors and atmospheric effects and also make the connection between separate subsets in the interfrogram network, a smoothing constraint is added to the least-square problem (Lundgren, 2001).Minimize root mean square (RMS) of the residual in least-square solution corresponding to smoothing factor as the objective function, determines Smoothing factor optimally.Time series analysis is finally carried out using the optimum smoothing factor (Dehghani, 2009).

Groundwater level monitoring at piezometric wells
Piezometric wells are being used to monitor the groundwater level's fluctuations and generate the water level decline map which can be used for subsidence interpretation (Calderhead, 2011).There are 53 piezometric wells in the study area which have been monthly monitored from Oct 2004 to Mar 2005.The results obtained from the InSAR time series are compared to the map demonstrating groundwater level decline.

RESULTS
As mentioned above, 18 differential interferograms were processing.Some of the interferograms, i.e.  4. Before least squares inversion, the remaining orbital error of interferograms must be eliminated by subtracting a plate fitted to the points away from the subsidence area.Moreover, a point is considered as a spatial reference.The mean displacement velocity map obtained from InSAR time series analysis are presented in Figure 5.a.The map highlights the major features of the subsidence.The maximum subsidence rate was estimated as 28 cm/yr.Various sinkholes and fractures were obvious in field observation in Sirjan plain.Figure 6 shows some of these evidences.
A map for water level changes in aquifer was produced over a  and was compared to the InSAR results in order to better interpret the subsidence occurrence.The correlation between the water level decline and aquifer compaction is clearly observed from Figure 5.For example, the area with the highest deformation rate corresponds to the area with the maximum groundwater level decline.This area is shown in a black square in Figure 5.  Water level fluctuations at a selection of piezometric wells whose locations are shown in Figure 5 are presented in Figure 7. Accordingly, a water level decline is observed in Wells 5 and 30 located in the middle of the subsidence area.On the other hand, well 16 which is located outside the subsidence area shows increase in water level.
It should be noticed that other important factors are at work in subsidence occurrence.Aquifer compaction is highly affected from the soil type.If an aquifer is composed of fine-grained sediments, a little water level decline causes a significant subsidence.Aquifer compaction modelling considering other hydrogeological information is considered as future work.
Figure 7. Water level changes of aquifer at piezometric wells located in different parts of the subsidence area.
In order to quantify the relationship between groundwater level fluctuations and ground surface subsidence, we showed both of them in a unique plot for some of piezometric wells in  The values of skeletal storage coefficients for all piezometric wells located in difference part of the subsidence area are calculated and shown in Table 3.

Conclusions
In this study, 12 ENVISAT ASAR images spanning between 2004 and 2010 were used for subsidence monitoring in Sirjan plain.The subsidence is induced by over-exploitation of groundwater for agricultural purposes.SBAS time series analysis was applied to generate the mean subsidence velocity map to highlight the major features of the deformation.Time series analysis results were compared with groundwater levels' information at piezometric wells.The comparison shows the high correlation between water level decline and aquifer compaction.For more accurate discussion, multi-sensor SAR data and modelling aquifer system considering geology and hydrogeology parameters are suggested.
a) Mean subsidence velocity map and b) water level changes in the aquifer system.Piezometric wells 5, 30 and 16 in the middle and outside the subsidence area are depicted.

Figure 6 .
Figure 6.Fractures in the Sirjan plain as consequences of subsidence.
Figure 8.As it is shown, the values of water level spanning between Oct 2004 and Mar 2005 are plotted on the y-axis and InSAR deformation time-series representing ground subsidence are plotted on the x-axis.Values shown on x-axis are over the span of measuring water level depth in piezometric wells, roughly.As shown in the figure, red solid line is a linear regression fitted to the plotted points.The inverse slop of solid line can be considered as an approximate of skeletal storage coefficient at the well location(Hoffmann, 2001).

Figure 8 .
Figure 8. Sample graphs for calculating the skeletal storage coefficient at piezometric wells located in different parts of the subsidence area.

Table 1 .
Acquisition date of the images available

Table 3 .
Values calculated for skeletal storage coefficient at wells location in over the subsidence area.