Geophys. J. Int. (2009) 178, 47–56 doi: 10.1111/j.1365-246X.2009.04135.

GJI Geodesy, potential field and applied geophysics

InSAR monitoring of progressive land subsidence in Neyshabour,
northeast Iran

Maryam Dehghani,1∗ Mohammad Javad Valadan Zoej,3 Iman Entezam,3

Ali Mansourian2 and Sassan Saatchi4
1 Facultyof Geomatics and Geodesy, K.N.Toosi University of Technology (KNTU), Tehran, Iran. E-mail: [email protected]
2 Facultyof Geomatics and Geodesy, K.N.Toosi University of Technology (KNTU), Tehran, Iran
3 Engineering Geology Group, Geological Survey of Iran (GSI), Tehran, Iran
4 UCLA Center for Tropical Research, Los Angeles, USA

Accepted 2009 January 31. Received 2009 January 31; in original form 2008 March 15

SUMMARY
The area of Neyshabour, a small historical city located in Northeast Iran, is subject to land
subsidence. To monitor the temporal evolution of the subsidence, the small baseline subset

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015

(SBAS) algorithm is used for interferometric synthetic aperture radar (SAR) time-series anal-
ysis. To limit the spatial and temporal decorrelation phenomena, the interferograms produced
from the raw ENVISAT ASAR data are characterized by small spatial and temporal baselines.
Accordingly, four independent SAR acquisition data sets separated by large spatial and tem-
poral baselines are used in the time-series analysis. To link the separate data sets, a smoothing
constraint that minimizes the curvature of the subsidence temporal evolution is added to the
least-squares method. The optimum smoothing factor estimated in the smoothed time-series
analysis reduces the atmospheric noise, unwrapping and orbital errors whereas it preserves the
non-linear seasonal deformation features. The time-series results show an incremental lower-
ing of the ground surface, accompanied by small seasonal effects. The mean LOS deformation
velocity map computed from the time-series analysis demonstrates a considerable subsidence
rate of up to 19 cm yr–1 . Comparison between the InSAR time-series and continuous GPS
measurements verifies the accuracy of the obtained results.
Moreover, the quantitative integration of the InSAR-derived displacement measurements
with observations of the hydraulic head fluctuations causing these displacements yields infor-
mation about the compressibility and storage properties of the aquifer system.
Key words: Time series analysis; Radar interferometry.

sediments are compacted and land subsidence occurs (Hoffmann

1 I N T RO D U C T I O N
et al. 2003). Environmental consequences of land subsidence mainly
Land-surface subsidence caused by overextraction of groundwater include damage to engineered structures (such as buildings, road-
has been recognized as a major problem with environmental con- ways, pipelines and well casings), earth fissures and surface runoff
sequences in many areas (e.g. Longfield 1932; Tolman & Poland (Hoffmann et al. 2003).
1940; Poland & Davis 1969; Galloway et al. 1998, 1999; Abidin During the past decade, the Neyshabour plain, located in the
et al. 2001). Excessive groundwater pumping results in an increas- Kavir Markazi basin in northeast Iran (Fig. 1), has experienced
ing effective stress within the sediments of the aquifer system significant land subsidence, associated with the compaction of the
(Terzaghi 1925). Hence, the aquifer system has undergone some aquifer system.
degree of deformation in response to the changes in stress. The This plain is bounded on the south by the eastern part of the
seasonal cycle of discharge and recharge from unconsolidated het- Alborz Mountains, called the Binaloud, which trends NW–SE. The
erogeneous aquifer systems causes compaction and expansion, re- Neyshabour plain is mostly covered by young Quaternary terraces.
spectively (Riley 1969; Poland & Ireland 1988; Heywood 1997; The area is bordered by thrust and overthrust faults to the north.
Amelung et al. 1999; Hoffmann et al. 2001; Lu & Danskin 2001). Annual rainfall in the area is estimated to range from 250 to
When water is removed from fine-grained, highly compressible 400 mm yr–1 , occurring mainly during spring and winter seasons
sediments, such as clay and silt interbeds in an aquifer, these (Ghaemi et al. 2002).
Water level measurements from piezometer installations have
been made monthly, by the Water Management Organization, since
∗ IEEE Student Member 1997. Excessive groundwater extraction from pumping wells to


C 2009 The Authors 47
Journal compilation 
C 2009 RAS
48 M. Dehghani et al.

Figure 1. (a) Location of the Neyshabour plain in Iran, (b) LANDSAT 7 ETM+ colour composite image (R:7, G:4, B:2) of the Neyshabour plain.

provide water for the growing cities and cultivated lands has low- compaction at long timescales. However, InSAR time-series anal-

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015

ered aquifer hydraulic heads. Moreover, the sustained decrease of ysis has been recently developed as an appropriate technique to
groundwater levels indicates that the aquifer is not sufficiently quantitatively characterize the spatial and short-term behaviour of
recharged. The Neyshabour plain is directly affected by this type of land subsidence (Berardino et al. 2002).
subsidence, caused by aquifer-system compaction associated with The main purpose of this paper is to show the potential of In-
the decline in hydraulic head. The water level decline in the aquifer SAR to study the mechanics of the Neyshabour aquifer system, as
system, which is composed of highly compressible fine-grained in- well as the spatial heterogeneity of the aquifer system structure, and
terbeds, has led to subsidence rates of several centimetres per year, monitor the temporal and spatial behaviour of compaction of the
resulting in damage to well casings and structures. Land subsi- Neyshabour basin. InSAR time-series analysis is used to measure
dence extending to Neyshabour City would be a cause for concern, the compaction of the aquifer system. To mitigate the atmospheric
since this city is historically important. Several buildings more than artefacts, orbital errors and noise associated with unwrapping, a
1000 yr old are located in Neyshabour, making it important to mit- smoothing constraint that minimizes the curvature of the temporal
igate the negative effects of land subsidence. evolution of subsidence is incorporated into the time-series anal-
To date, no comprehensive study on subsidence in Neyshabour ysis. The stress–strain relationships are then derived using InSAR
and its aquifer system has been carried out. Due to the importance measurements and groundwater information from the piezometric
of Neyshabour as a historical city, monitoring the land subsidence is wells to understand the mechanics of the compaction of the aquifer
required. Subsidence in Neyshabour has previously been monitored system as a function of water level decline. The quantitative in-
using precise levelling observations. A precise levelling survey car- tegration of the InSAR data and groundwater observations yields
ried out from 1992 to 2002 recognized a maximum subsidence rate information about the storage properties of the aquifer system.
of 15 cm yr–1 (Amighpey et al. 2006). Furthermore, a continu- This paper is organized as follows. Section 2 explains the InSAR
ous GPS station (NISH) has registered a cumulative subsidence of time-series analysis algorithm. The time-series results are presented
23 cm from 2005 to 2007. Precise levelling observations and contin- in Section 3. Section 4 is devoted to the quantitative comparison
uous GPS measurements measure the deformation at a few sparse of the InSAR results and groundwater information. Concluding
station locations. Measurements of aquifer compaction resulting remarks are presented in Section 5.
from the hydraulic head decline can reveal information about the
aquifer compressibility and storage characteristics.
The lack of detailed spatial geodetic information has limited the
2 R A D A R I N T E R F E R O M E T RY
study of the spatial heterogeneity of the aquifer system. Among
T I M E - S E R I E S A N A LY S I S
various ground- and space-based techniques available, the Interfer-
ometric SAR technique (InSAR) has shown the ability to provide To produce the deformation time-series, 14 Level0 ENVISAT
precise measurements of land surface deformation over large ar- ASAR images, spanning 2004 and 2005, are used. The data are
eas and at high spatial resolution (Massonnet et al. 1993; Peltzer converted into single look complex (SLC) format using ESA or-
et al. 1998; Fruneau & Sarti 2000; Tesauro et al. 2000; Crosetto bital data and GAMMA software (Wegmuller et al. 1995). The data
et al. 2002; Chatterjee et al. 2006; Motagh et al. 2006). A number set collection with respect to the temporal and spatial baseline is
of studies have used the InSAR technique to detect and monitor illustrated in Fig. 2. Since the study area is covered by vegetation,
land subsidence (Chang et al. 2004; Carbognin et al. 2004; Mo- the temporal decorrelation increases significantly with time (Zebker
tagh et al. 2006). High-resolution InSAR measurements offer great & Villasenor 1992). The temporal baseline greater than 4 months
potential for studying the spatial heterogeneity of aquifer system results in very low coherence in the cultivated parts of the study
properties and monitoring the land subsidence (Hoffmann et al. area, whereas the coherence is well preserved in non-vegetated
2001). Lu & Danskin (2001) and Bawden et al. (2001) defined areas. Therefore, the generated interferograms are characterized
the structure of groundwater basins using InSAR data. Traditional by small temporal and spatial baselines to reduce the temporal
InSAR techniques have been widely used to characterize aquifer and baseline decorrelation, respectively. Fourteen interferograms,


C 2009 The Authors, GJI, 178, 47–56
Journal compilation 
C 2009 RAS
 Land subsidence in Neyshabour, northeast Iran 49

artefacts and noise generated by the unwrapping process rather

than real uplift, as they appear randomly in the interferograms. The
signal of the subsidence can be easily recognized in the different
interferograms due to the constant geographic locations. Other sig-
nals corresponding to error sources, such as atmospheric artefacts,
do not have fixed positions over a long period of time (Massonnet
& Feigl 1998).
The next step is to generate the deformation time-series, using the
interferograms. The main idea in time-series analysis is to invert the
interferograms to obtain the deformation at each date of acquisition
using the least-squares method. Suppose D = [d 1 , d 2 , d 3 . . .] is the
vector of the deformation on each acquisition date to be estimated
and I = [i 12 , i 23 , i 34 . . .] is the observation vector containing the
Figure 2. Available ENVISAT ASAR images for the Neyshabour plain. range changes obtained from the interferograms. These vectors can
Solid lines represent the generated interferograms.
be related together by eq. (1):
representing four independent data sets separated by large base- AD = I, (1)
lines, are generated using the GAMMA software.
To remove topographic effects, a digital elevation model from the where A is the design matrix. In the inversion solution, the deforma-
shuttle radar topography mission (SRTM; Farr & Kobrick 2000), tion of the first date was set to zero (d 1 = 0). As mentioned before,
with a spatial resolution of 90 m, is used. The orbital errors due to there are four independent data sets, separated by large spatial and
the tilts and offsets remaining in the interferograms are removed by temporal baselines. To easily link these independent data sets and

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015

subtracting a plane fitted to the data in the far field, away from the mitigate several error types, including atmospheric artefacts, noise
deformation signal (e.g. Funning et al. 2005). To test and estimate and unwrapping errors, the smoothing constraint is incorporated
the uncertainty in the data, the root mean square (rms) misfit be- into the inversion problem (e.g. Lundgren et al. 2001; Schmidt &
tween the data and the fitted plane in the far field is calculated for Burgmann 2003). Using a finite difference approximation for the
all interferograms. Table 1 presents the uncertainty in the interfero- second-order differential of the time-series as a smoothing con-
grams. straint, eq. (1) is written as
   
According to Table 1, a misfit on the order of a millimetre for all A I
D= , (2)
interferograms is representative of high performance of referencing γ ∂ /∂t
2 2 2
0
the phase in the far field to zero. Fig. 3 presents several unwrapped
interferograms, recording the incremental deformation between two where γ is the smoothing factor. A new smoothing constraint must
acquisition dates. A few unwrapping errors appear clearly on the be added for each additional acquisition date. If the acquisition
interferograms because of low-correlated areas. Adaptive filter is dates are not evenly spaced, which is usually the case, the irregu-
used to decrease the number of residues and avoid unwrapping lar time spacing must be included in the finite difference expres-
errors. sion by applying the minimum curvature concept, that is, constant
The interferometric phase is converted into displacement along velocity.
the line of sight (LOS) direction via multiplication by the correction The smoothing factor, γ , should be determined optimally. The
factor λ/4π . Since the incidence angle of the radar is quite small optimal estimation of the smoothing factor results in the smooth
(23◦ ) and the major part of the deformation resulting from land sub- deformation time-series whereas the non-linear deformation com-
sidence can be assumed to be vertical, the majority of the surface ponents are preserved. The time-series results can become rough by
deformation caused by subsidence can be measured in the computed decreasing the smoothing factor. In this case, the observed fluctu-
interferograms. The negative values in the interferograms represent ations in the results are mostly due to the sources of error in each
subsidence signals. It should be noted that the large positive val- interferogram. If a large smoothing factor is selected, any non-linear
ues observed in all interferograms are mostly due to atmospheric deformation will be damped, and the overweighted smoothing will
force the time-series to be a straight line. If the smoothing is in-
Table 1. Root-mean-square misfit between the data and far-field fitted plane. finite, the best-fit average deformation rate is achieved. Indeed, a
trade-off methodology should be applied to select the most appro-
Interferograms Data uncertainty: Far-field rms (m)
priate smoothing factor regarding reduction of the errors as well as
2004/11/10–2004/05/29 0.00929 preservation of the subtle deformation signal. There are several ap-
2004/05/29–2004/10/16 0.00655 proaches to determine the optimum smoothing factor; one of these
2005/05/14–2005/07/23 0.00737 is to plot the overall rms (misfit) that results from least-squares so-
2005/07/23–2005/11/05 0.00921
lution against different corresponding smoothing factors. The rms
2004/04/24–2004/08/07 0.0078
2004/11/20–2004/12/25 0.00325
is the root mean square of the residuals in least-squares solution as
2004/12/25–2005/01/29 0.0027 follows:

2004/12/25–2005/03/05 0.0057 
 n
2004/12/25–2005/04/09 0.0074  (r̂i )2
2005/03/05–2005/04/09 0.00607  i=1
rms = , (3)
2005/03/05–2005/06/18 0.00686 n
2005/04/09–2005/06/18 0.00515
2004/11/20–2005/01/29 0.00539 where, according to eq. (1), r̂ is defined as
2005/01/29–2005/03/05 0.00190
r̂ = AD̂ − I, (4)


C 2009 The Authors, GJI, 178, 47–56
Journal compilation 
C 2009 RAS
50 M. Dehghani et al.

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015

Figure 3. Unwrapped interferograms of subsidence in Neyshabour superimposed on the LANDSAT ETM+ image. Time intervals covered by the interferograms
are: (a) 2004/05/29 to 2004/10/16; (b) 2004/04/24 to 2004/08/07; (c) 2004/11/20 to 2004/12/25 and (d) 2005/03/05 to 2005/06/18.

3 R E S U LT S
To highlight the major deformation features in the area, the mean
displacement velocity map is extracted using the InSAR time-series
results. The mean displacement velocity map contains the subsi-
dence rate along the LOS direction. The extracted mean displace-
ment velocity map, superimposed on the shaded relief map of the
area, is shown in Fig. 5.

Figure 4. Optimum smoothing factor selected at the elbow of the curve.

Horizontal and vertical axes represent the smoothing factor and least-squares
misfit, respectively.

where D̂ is the estimated phase via the least-squares solution. The

misfit will be minimum with a zero smoothing factor and will
increase with higher values. The resulting plot is a curve with an
elbow. Although the choice of the optimum smoothing factor is
arbitrary, depending on the plot scale, the middle point on the curve
is an estimate for the smoothing factor, which is a good compromise
between removal of noisy fluctuations and preservation of the non-
linear seasonal deformation. Fig. 4 shows the plot of misfit against
smoothing factor and the best smoothing factor. Figure 5. Mean displacement velocity map superimposed on the shaded
Finally, the entire time-series process is repeated using the opti- relief extracted from the SRTM DEM. Yellow to red colours represent sub-
mally selected smoothing factor. sidence.


C 2009 The Authors, GJI, 178, 47–56
Journal compilation 
C 2009 RAS
 Land subsidence in Neyshabour, northeast Iran 51

Figure 6. Average subsidence maps in (a) winter and (b) summer.

As observed in Fig. 5, the maximum deformation rate is not steady in time, including seasonal effects (e.g. Lanari et al.
19 cm yr–1 . However, the subsidence rates in the discharge and 2004). To study the evolution of deformation, time-series plots of a
recharge seasons are not identical. To investigate the subsidence selection of points located within, along the margin of and outside
rate in the different seasons, average subsidence maps for winter the subsidence area were generated (Figs 8a–c). These points are
2004–2005 and summer 2005 are produced (Fig. 6). This reveals the piezometric wells shown in Fig. 5.

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015

a maximum accelerating subsidence rate of 17 cm yr–1 in summer The deformation time-series of the wells located outside the sub-
and a decelerating deformation rate of 11 cm yr–1 in winter. sidence area, that is, wells 6, 7 and 9, do not show the long-term
Faults and structural lineaments of the area extracted from the subsidence trend, rather seasonal fluctuations (Fig. 8a). However,
1 : 100, 000 geology map and satellite images (LANDSAT 7 ETM+ two other groups of points illustrate the subsidence signals on
and RADARSAT1 images in ScanSAR narrow beam mode) are il- which the seasonal effects are superimposed (Figs 8b and c). Ac-
lustrated in Fig. 5. As observed in this figure, the subsidence pattern cording to the deformation time-series plots of wells 5, 8, 17, 33,
is strictly controlled by the faults and structures surrounding the 14, 16, 31 and 32, relatively minor decelerating subsidence can
area. These faults and structures act as barriers in the aquifer sys- be observed during winter. In the first three acquisition dates of
tem that controls the groundwater flow. Furthermore, they may be the time-series (2004/01/10, 2004/0/24 and 2004/05/29), there is
considered as the boundaries of material zones with different con- decelerating subsidence or, perhaps, uplift related to the recharge
trolling hydrogeological parameters that affect the subsidence, such season. This can be observed in the time-series plots in Figs 8(a)–
as specific storage and hydraulic conductivity. Lineaments A and (c). From the 4th to the 6th acquisition dates in the deforma-
B in Fig. 5 agree well with the boundaries of the subsidence area. tion sequence (2004/08/07, 2004/10/16, and 2004/11/20), accel-
An abrupt change in subsidence rate from 13 to 7 cm yr–1 can be erating subsidence is involved because of the discharge season.
observed across lineament C. Similarly, the subsidence rate changes During this period, subsidence has its maximum rate, as shown in
from 9 to 5 cm yr–1 and from 5 to 0.2 cm yr–1 across lineaments Figs 8(a)–(c). The next recharge season is from the 7th to 11th acqui-
D and E, respectively. More interestingly, the comparison between sition dates (2004/12/25, 2005/01/29, 2005/03/05, 2005/04/09 and
the LANDSAT ETM+ image of the Neyshabour plain and the sub- 2005/05/14). As expected, and as can also be seen in Figs 8(a)–(c),
sidence zone indicates that the subsidence pattern exactly follows decelerated subsidence is observed in this time period. Finally, the
the Neyshabour plain fan. Geophysical and geotechnical investiga- last period (2005/06/18, 2005/07/23 and 2005/11/05) corresponds
tions that are not presently available, may reveal information about to discharge season. The accelerated subsidence is clearly observed
the underlying layers in the aquifer system, making further studies during this period in Figs 8(a)–(c). Involvement of the seasonal dis-
about the roles of faults possible. charge and recovery of the aquifer system in the InSAR time-series
The mean displacement velocity map is extracted under the results is representative of optimum smoothing, not damping out of
assumption of infinite smoothing. Therefore, any non-linear de- the seasonal surface fluctuations.
formation, including atmospheric artefacts, has been damped out. To evaluate the results of the time-series analysis, continuous
To verify the optimum smoothing in the time-series analysis, sev- GPS measurements were used. The National Cartographic Center
eral subsidence maps are generated using the time-series Fig. 7). (NCC) of Iran, as a part of the Iranian Permanent Network for
Figs 7 (a) and (c) show raw interferograms applied in the InSAR Geodynamic (IPGN) programme, has established a continuously
time-series analysis. Deformation sequences obtained by the time- recording GPS station (NISH) near Neyshabour City. The location
series analysis are used to generate the displacement maps spanning of the NISH GPS station is shown by a black square in Fig. 5.
the same time intervals as the interferograms (Figs 7b and d). This station has recorded data continuously since 2005. The InSAR
In the generated subsidence maps, the artefacts generally con- time-series for the location of the NISH GPS station, in addition
sidered to be atmospheric, are significantly decreased, which shows to the continuous GPS measurements, are illustrated in Fig. 9. The
the effect of smoothing. On the other hand, it will be shown below GPS measurements and the InSAR time-series for the short period
that the degree of smoothing does not significantly damp out the of overlap agree only in the trend from 2005.5 to 2006. Because of
seasonal fluctuations. the large scatter of GPS measurements until 2006, it is not possible
The mean velocity map allows for the identification of the long- to compare the seasonal fluctuations of the InSAR-derived results
term average subsidence. However, the main strength of the time- with those of the GPS data. If the period of overlap between both
series analysis is its potential to detect deformation signals that are data sets were longer, their agreements would be more obvious.


C 2009 The Authors, GJI, 178, 47–56
Journal compilation 
C 2009 RAS
52 M. Dehghani et al.

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015

Figure 7. Comparison of the interferograms and displacement maps generated from the time-series results. (a) Interferogram (2005/05/14 to 2005/07/23)
involving atmospheric artefacts, (b) displacement map (2005/05/14 to 2005/07/23) with mitigated atmospheric artefact, (c) interferogram (2004/04/24 to
2004/08/07) involving atmospheric artefacts and (d) displacement map (2004/04/24 to 2004/08/07) with mitigated atmospheric artefact.

The root mean square error (RMSE) between the InSAR time- To study the effective stress changes in the Neyshabour aquifer,
series results and the GPS measurements is estimated to be 0.0104 water level information from piezometric wells located in the subsi-
m. It should be noted that the GPS measurements are projected dence area is analysed. Water level measurements have been made
along the LOS direction, assuming an incidence angle of 23◦ . monthly by the Water Management Organization since 1997.
In addition to the GPS measurements, precise levelling data were Figs 8(d)–(f) present groundwater level fluctuations of the wells,
collected by NCC in two periods, that is, 1992 and 2002 (Amighpey whose corresponding InSAR time-series are plotted in Figs 8(a)–
et al. 2006). The levelling track is depicted in Fig. 5. Fig. 10 shows (c). The water level of piezometric wells has been set to zero at the
the deformation rates along the levelling track derived from the earliest time, that is, 1997, showing the water level changes over
InSAR results, as well as the precise levelling survey. The different time. In this case, any value less than zero corresponds to the water
time intervals covered by the InSAR data (2004 and 2005) and the level decline. The piezometric wells in Figs 8(d)–(f) mostly demon-
levelling survey (1992 and 2002) make the extracted deformation strate the decline in water level caused by excessive groundwater
rates from the two methods different. According to Fig. 10, the extraction. According to Fig. 8(d), the water level of piezometric
subsidence rate has increased during the InSAR time interval (from wells located outside the subsidence area (except for well 6) does
2004 to 2005). However, similar patterns are seen in both profiles. not show a significant decline. The maximum water level decline
of 15 m during 10 yr occurred in well 5 (Fig. 8e), reflecting the
high subsidence rate as well (Fig. 8b). The water level measured
4 DISCUSSION
at well 17 shows fluctuations rather than a falling trend. However,
The decreasing trend observed in the InSAR time-series at the loca- such fluctuations are not reflected in the corresponding deforma-
tions of the pumping wells indicates that the aquifer is compacting at tion time-series (Fig. 8b). This effect can also be observed in well
a constant long-term rate. An aquifer system is typically composed 33. In contrast, the water level has declined continuously in well 6,
of a series of low permeability and highly compressible fine-grained though it is located outside the subsidence area. Hence, it can be
interbeds. The theoretical basis of interbed compaction is based concluded that land subsidence is a function of not only the wa-
on the Terzaghi’s (1925) principle of effective stress. Accordingly, ter head decline but also other factors. The best way to investigate
changes in effective stress will result from changes in the total stress the relationship between groundwater level fluctuations and surface
or changes in pore pressure. If the total stress, which is the geostatic displacements is to map them in a unique plot. It means that the wa-
load of the overlying saturated and unsaturated sediments and tec- ter level variations representing the stress are plotted on the y-axis,
tonic stresses, is assumed to remain constant, a change in effective whereas the InSAR deformation time-series showing the aquifer
stress is accompanied by a change in pore pressure. Groundwa- compaction are plotted on the x-axis. The inverse slope of the best
ter extraction results in decreasing pore pressures and increasing fitting line to the plotted points is a rough estimate of the skeletal
effective stress, which act to compress the interbeds. storage coefficient of the aquifer system (Hoffmann et al. 2001).


C 2009 The Authors, GJI, 178, 47–56
Journal compilation 
C 2009 RAS
 Land subsidence in Neyshabour, northeast Iran 53

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015

Figure 8. (a)–(c) InSAR time-series at piezometric wells located outside the subsidence area (w6, w7 and w9), along the margin of the subsidence area (w5,
w8, w17 and w33) and in the middle of the subsidence area (w14, w16, w31 and w32), respectively; (d)–(f) Water table fluctuations of piezometric wells
located outside the subsidence area (w6, w7 and w9), along the margin of the subsidence area (w5, w8, w17 and w33) and in the middle of the subsidence area
(w14, w16, w31 and w32), respectively.

Figure 9. Comparison of GPS (blue dots) and InSAR (red solid line) time- Figure 10. Comparison of subsidence rates inferred from InSAR and pre-
series at the NISH station. cise levelling.

The storage coefficient of an aquifer system is a parameter that as well as the surface displacement, is not significant. In well 6,
contains the responses of the aquifer and fine-grained interbeds to there is no significant subsidence signal despite the water level de-
variations in hydraulic head. Using this approach, the aquifer system cline. In wells 5, 8 and 17, located at the margin of the subsidence
compaction mechanism can be derived from the InSAR displace- area, ground surface lowering is linearly dependant on water level
ment time-series and contemporaneous measurements of water lev- decline. However, the water level in well 33 remains unchanged,
els in piezometric wells. Fig. 11 illustrates the water level variations though its subsidence rate is considerable. Depending on the hydro-
plotted versus the ground displacements for the piezometric wells. geology, the time-consolidation history and characteristics of the
Since the water level observations are made monthly, ground dis- fine-grained interbeds, including low vertical hydraulic conductiv-
placements are linearly interpolated to the water level observation ity, residual compaction may occur in thicker interbeds, so that even
dates. though there may be seasonal water-level recovery, the thicker in-
As shown in Fig. 11, the ground surface is typically lowered by terbeds may continue to compact. The same effect can be observed
the decline of water level. In wells 7 and 9, the water level decline, in well 31, which is located within the subsidence area. Other wells


C 2009 The Authors, GJI, 178, 47–56
Journal compilation 
C 2009 RAS
54 M. Dehghani et al.

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015

Figure 11. Calculation of the storage coefficient from stress displacement analysis. Data are plotted in a stress–strain diagram. The inverse slope of the best
fitting line (solid line) is the least-squares estimate of the storage coefficient for each piezometric well. ‘S’ in the plots corresponds to the storage coefficient.

located in the middle of the subsidence area (14, 16 and 32) show estimates of the storage coefficients, additional knowledge of the
linear relationships between water level decrease and aquifer com- interbed thicknesses and distributions in the aquifer system is re-
paction. The solid lines in Fig. 11 correspond to the regression line quired. Furthermore, using the hydraulic heads calculated from a
fitted to the stress displacement data. The storage coefficients at groundwater flow model, the InSAR-derived displacements can be
wells’ location estimated by the inverse slop of the regression line applied to generate a map of storage coefficients. This will be con-
are shown in Fig. 11. This parameter can also be determined by sidered in future investigations.
pumping tests or in the laboratory from core samples, which are
generally costly. In the former method, the extracted storage coef-
ficients are representative of only the most permeable fraction of
5 C O N C LU S I O N S
the aquifer system, whereas in the latter, the measurements are not
representative of in situ conditions. However, the values extracted In this case study, the potential of Interferometric SAR (InSAR)
by the relationship between the InSAR-derived displacements and in hydrogeological applications has been shown. In the first part
groundwater levels yield spatially varying estimates of storage co- of the paper, the InSAR approach is used to monitor land subsi-
efficients at the well locations. The estimated values of the slope dence induced by groundwater overexploitation in the Neyshabour
of the best fitting line allow us to predict the amount of subsidence plain. Fourteen ENVISAT ASAR Level0 images are used to gener-
caused by water level decline. The computed storage coefficient ate 14 interferograms characterized by small spatial and temporal
for each piezometric well is shown in Fig. 11. The negative values baselines. The processed interferograms split into four indepen-
are incorrect estimates, which are due to the residual compaction dent data sets separated by large spatial and temporal baselines.
occurring in thicker interbeds. To link these independent data sets and mitigate the atmospheric
Because of the characteristically low vertical hydraulic conduc- artefacts, a smoothed time-series analysis is carried out, using an
tivity of compressible sediments that make up the interbeds of the optimal smoothing factor. The InSAR time-series results agree well
aquifer system, there is a time delay between the water table decline with the trend of GPS measurements during their short period of
and compaction of the interbeds. Factors that control the timing of overlap. A maximum deformation rate of 19 cm yr–1 is estimated
the land subsidence include the vertical diffusivity and the thickness from the mean displacement velocity map extracted from the In-
of the fine-grained interbeds (Hoffmann et al. 2003). Because of the SAR time-series results. The InSAR time-series for the locations
time lag involved in the equilibration of the fine-grained interbeds, of a selection of piezometric wells distributed in different parts of
and consequently, the lower stress changes of the interbeds than is the subsidence zone show the subsidence signal with a long-term
imposed, the estimated values of storage coefficients do not exactly constant rate accompanied by seasonal effects. Groundwater infor-
reflect the elastic compressibility of the aquifer system; however, mation from piezometric wells showed a large decline over 10 yr,
they are an acceptable estimate. In addition to the time delay, other resulting in compaction of the sediments of the aquifer system and
hydrogeological parameters such as hydraulic conductivity, specific land subsidence.
storage, permeability and depth of underlying interbeds will define In the second part of the paper, surface displacements measured
the correlation between hydraulic head changes and compaction by InSAR and water level observations were used to determine
occurring in the fine-grained interbeds. To obtain more accurate the stress–strain relationship at the piezometric wells. Quantitative


C 2009 The Authors, GJI, 178, 47–56
Journal compilation 
C 2009 RAS
 Land subsidence in Neyshabour, northeast Iran 55

integration of groundwater level information and InSAR-derived Carbognin, L., Teatini, P. & Tosi, L., 2004. Eustacy and land subsidence in
displacements at the locations of piezometric wells yield compress- the Venice Lagoon at the beginning of the new millennium, J. Mar. Syst.,
ibility properties of the aquifer system. In most wells, the strain 51, 345–353.
increases with the decrease in water level. However, depending on Chang, C.P., Chang, T.Y., Wang, C.T., Kuo, C.H. & Chen, K.S., 2004. Land-
surface deformation corresponding to seasonal ground-water fluctuation,
the hydrogeology, the time-consolidation history and characteris-
determining by SAR interferometry in the SW Taiwan, Math. Comput.
tics of the fine-grained interbeds, residual compaction may occur in
Sim., 67, 351–359.
thicker interbeds. In this case, even though there may be seasonal Chatterjee, R.S. et al., 2006. Subsidence of Kolkata (Calcutta) City, India
water-level recovery, the thicker interbeds may continue to compact. during the 1990s as observed from space by Differential Synthetic Aper-
The non-linear relation between the InSAR time-series and the wa- ture Radar Interferometry (D-InSAR) technique, Remote Sens. Environ.,
ter level changes in the piezometric wells suggests that various 102, 176–185.
hydrogeological and soil properties, such as hydraulic conductivity, Crosetto, M., Tscherning, C.C., Crippa, B. & Castillo, M., 2002. Sub-
specific storage, permeability and depth of underlying interbeds, sidence monitoring using SAR interferometry: reduction of the atmo-
as well as faults and structures, may affect the subsidence rate and spheric effects using stochastic filtering, Geophysi. Res. Lett., 29(9), 26.1–
pattern. To model the subsidence and aquifer-system compaction, 26.4.
Farr, T. & Kobrick, M., 2000. Shuttle Radar Topography Mission produces
complementary geotechnical studies are required to determine the
a wealth of data, Eos Trans. Am. geophys. Un., 81, 583–585.
hydrogeology properties of geological units. The combination of In-
Fruneau, B. & Sarti, F., 2000. Detection of ground subsidence in the city
SAR data, providing the surface subsidence, and hydrogeological, of Paris using radar interferometry: isolation from atmospheric artifacts
lithologic and geotechnical information to model the compaction of using correlation, Geophys. Res. Lett., 27(24), 3981–3984.
the aquifer-system, is considered as future work. Funning, G.J., Parsons, B., Wright, T.J. & Jackson, J.A., 2005. Surface
Land subsidence due to aquifer compaction is generally con- displacements and source parameters of the 2003 Bam (Iran) earthquake
sidered as an environmental hazard and has several consequences, from Envisat advanced synthetic aperture radar imagery, J. geophys. Res.,

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015

such as damage to buildings, railways and pipelines, and can result 110, B09406, doi:10.1029/2004JB003338.
in earth fissures (Shemshaki et al. 2005). Furthermore, the sustained Galloway, D.L., Hudnut, K.W., Ingebritsen, S.E., Phillips, S.P., Peltzer, G.,
water level decline may be considered as a cause for concern since Rogez, F. & Rosen, P.A., 1998. Detection of aquifer system compaction
and land subsidence using interferometric synthetic aperture radar, An-
the aquifer system is not being replenished sufficiently. Hence, an
telope valley, Mojave Desert, California, Water Resour. Res., 34, 2573–
efficient management of groundwater exploitation is required. Re-
2585.
cently, the Water Management Organization has launched a novel Galloway, D.L., Jones, D.R. & Ingebritsen, S.E., 1999. Land subsidence in
program to limit the drilling of wells and groundwater extraction. the United States, US Geological Survey Circular, No.1182, 175 p.
This may allow the aquifer system to recharge and mitigate the Ghaemi, F., Hosseini, K. & Ghaemi, F., 2002. “ Report of Neyshabour
subsidence rate and its consequences. geology map with the scale of 1:100, 000”, No. 7762, Geological Survey
of Iran.
Heywood, C.E., 1997. Piezometric-extensometric estimations of specific
AC K N OW L E D G M E N T S storage in the Albuquerque Basin, New Mexico, in U.S. Geological Sur-
vey Subsidence Interest Group Conference: Proceedings of the technical
The presented work was supported by the Geological Survey of meeting, U.S. Geological Survey Open-File Report 97–47, pp. 21–26.
Iran (GSI). The authors wish to thank the European Space Agency Hoffmann, J., Galloway, D.L., Zebker, H.A. & Amelung, F., 2001. Sea-
(ESA) and Remote Sensing group of the GSI for providing the sonal subsidence and rebound in Las Vegas Valley, Nevada, observed by
ENIVSAT ASAR data and GAMMA software, respectively. They synthetic aperture radar interferometry, Water Resour. Res., 37(6), 1551–
convey their sincere gratitude to the Water Management Organiza- 1566.
tion of Mashhad for providing the groundwater information used in Hoffmann, J., Leake, S.A., Galloway, D.L. & Wilson, A.M., 2003.
this work. The authors also wish to thank the National Cartographic MODFLOW-2000 ground-water model—user guide to the subsidence
and aquifer-system compaction (SUB) package, U.S. Geological Survey
Center of Iran (NCC) for collecting the continuous GPS measure-
Open-File Report 03-233.
ments and precise levelling data. Suggestions and comments by the Lanari, R., Lundgren, P., Manzo, M. & Casu, F., 2004. Satellite radar inter-
Engineering Geology group of the GSI are appreciated as well. ferometry time series analysis of surface deformation for Los Angeles,
California, Geophys. Rese. Lett., 31, doi:10.1029/2004GL021294.
Longfield, T.E., 1932. The subsidence of London, Ordnance Survey Profes-
sional Papers, no. 14, HMSO, London.
REFERENCES
Lu, Z. & Danskin, W.R., 2001. InSAR analysis of natural recharge to define
Abidin, H.Z. et al., 2001. Land subsidence of Jakarta (Indonesia) and its structure of a ground-water basin, San Bernardino, California, Geophys.
geodetic monitoring system, Nat. Hazards, 23, 365–387. Res. Lett., 28(13), 2661–2664.
Amelung, F., Galloway, D.L., Bell, J.W., Zebker, H.A. & Laczniak, R.J., Lundgren, P., Usai, S., Sansosti, E., Lanari, R., Tesauro, M., Fornaro, G. &
1999. Sensing the ups and downs of Las Vegas—InSAR reveals structural Berardino, P., 2001. Modeling surface deformation observed with SAR
control of land subsidence and aquifer-system deformation, Geology, interferometry at Campi Flegrei Caldera, J. geophys. Res., 106, 19 355–
27(6), 483–486. 19 367.
Amighpey, M., Arabi, S., Talebi, A. & Djamour, 2006. Elevation changes of Massonnet, D. & Feigl, K.L., 1998. Radar interferometry and its ap-
the precise leveling tracks in the Iran leveling network, Scientific report plication to changes in the Earth’s surface, Rev. Geophys., 36, 441–
published in National Cartographic Center (NCC) of Iran, Tehran, Iran. 500.
Berardino, P., Fornaro, G., Lanari, R. & Sansosti, E., 2002. A New algo- Motagh, M., Djamour, Y., Walter, T.R., Wetzel, H.U., Zschau, J. & Arabi, S.,
rithm for surface deformation monitoring based on small baseline differ- 2006. Land subsidence in Mashhad Valley, northeast Iran: results from
ential SAR interferograms, IEEE Trans. Geosci. Remote Sens., 40, 2375– InSAR, levelling and GPS, Geophys. J. Int., 168, doi: 10.1111/j.1365–
2383. 246X.2006.03246.x.
Bawden, G.W., Thatcher, W., Stein, R.S., Hudnut, K.W. & Peltzer, G., 2001. Peltzer, G., Rosen, P., Rogez, F. & Hudnut, K., 1998. Poro-elastic rebound
Tectonic contraction across Los Angeles after removal of groundwater along the Landers 1992 earthquake surface rupture, J. geophys. Res., 103,
pumping effects, Nature, 412, 812–815. 30 131–30 145.


C 2009 The Authors, GJI, 178, 47–56
Journal compilation 
C 2009 RAS
56 M. Dehghani et al.

Poland, J.F. & Davis, G.H., 1969. Land subsidence due to withdrawal of Terzaghi, K., 1925. Principles of soil mechanics, IV: settlement and consol-
fluids, Rev. Eng. Geol., 2, 187–269. idation of clay, Eng. News Rec., 95(3), 874–878.
Poland, J.F. & Ireland, R.L., 1988. Land subsidence in the Santa Clara Tesauro, M., Beradino, P., Lanari, R., Sansoti, E., Fornaro, G. &
Valley, California, as of 1982, U.S. Geological Survey Professional Paper, Franceschetti, G., 2000. Urban subsidence inside the City of Napoli (Italy)
No.497-F, 61 p. observed with synthetic aperture radar interferometry at Campi Flegrei
Riley, F.S., 1969. Analysis of Borehole Extensometer Data from Central caldera, J. geophys. Res., 27, 1961–1964.
California, International Association of Scientific Hydrology Publication, Tolman, C.F. & Poland, J.F., 1940. Ground-water, salt-water infiltration,
Vol.89, pp. 423–431, International Association of Scientific Hydrology. and ground-surface recession in Santa Clara Valley, Santa Clara County,
Schmidt, D.A. & Burgmann, R., 2003, Time-dependent land uplift and California, Am. geophys. Un. Trans., 21, 23–34.
subsidence in the Santa Clara valley, California, from a large interfero- Wegmüller, U., Werner, Ch., Strozzi, T., Wiesmann, A. & Santoro, M., 1995.
metric synthetic aperture radar data set, J. geophys. Res., 108(B9), 2416, GAMMA Remote Sensing Research and Consulting AG, Worbstrasse
doi:10.1029/2002JB002267 225, CH-3073 Gümligen, Switzerland, available at http://www.gamma-
Shemshaki, A., Blourchi, M.J. & Ansari, F., 2005. Preliminary report on rs.ch.
Tehran subsidence, Engineering Geology Scientific report on Geological Zebker, H.A. & Villasenor, J., 1992. Decorrelation in interferometric radar
Survey of Iran, available at www.gsi.ir. echoes, IEEE Trans. Geosci. Remote Sens., 30(5), 950–959.

Downloaded from http://gji.oxfordjournals.org/ by guest on July 7, 2015


C 2009 The Authors, GJI, 178, 47–56
Journal compilation 
C 2009 RAS