Nonlinear Processes in Geophysics Investigation of correlation of the variations in land subsidence ( detected by continuous GPS measurements ) and methodological data in the surrounding areas of Lake Urmia

Lake Urmia, a salt lake in the north-west of Iran, plays a valuable role in the environment, wildlife and economy of Iran and the region, but now faces great challenges for survival. The Lake is in immediate and great danger and is rapidly going to become barren desert. As a result, the increasing demands upon groundwater resources due to expanding metropolitan and agricultural areas are a serious challenge in the surrounding regions of Lake Urmia. The continuous GPS measurements around the lake illustrate significant subsidence rate between 2005 and 2009. The objective of this study was to detect and specify the nonlinear correlation of land subsidence and temperature activities in the region from 2005 to 2009. For this purpose, the cross wavelet transform (XWT) was carried out between the two types of time series, namely vertical components of GPS measurements and daily temperature time series. The significant common patterns are illustrated in the high period bands from 180–218 days band ( ∼ 6–7 months) from September 2007 to February 2009. Consequently, the satellite altimetry data confirmed that the maximum rate of linear trend of water variation in the lake from 2005 to 2009, is associated with time interval from September 2007 to February 2009. This event was detected by XWT as a critical interval to be holding the strong correlation between the land subsidence phenomena and surface temperature. Eventually the analysis can be used for modeling and prediction purposes and probably stave off the damage from subsidence phenomena.


Introduction
The subsidence of the Earth's surface is a phenomenon that occurs in some places in the world, which overuses underground sources of water.Iran has a semi-arid and arid climate and the rate of rainfall is lower than the mean rate in the world.Moreover, it is faced with over-exploitation of groundwater in agricultural areas and extension of the cities and industrial areas.Thus, increased requirement for groundwater, and the high rate of subsidence resulting from overutilization of this valuable resource are likely to become a serious challenge for the future development of the groundwater basins of the north-west of Iran.Consequently, unorganized groundwater extraction and the associated subsidence in the future can also endanger non-agricultural areas, resulting in costly damages to infrastructure due to changing ground levels.While damage to infrastructure can be repaired, the impact of depleted groundwater supplies will be more harmful, with immense social, economic and environmental consequences.Future climate change is expected to put additional stress on ground water resources in Iran (Sarraf et al., 2005).
Due to recent years of progressive dry climate, the water level of Lake Urmia, has been dramatically dropped about 6 m during the last 20 yr (Jalili et al., 2011).This lake, as a closed lake, is located in an area where evaporation is greater than rainfall.It obtains its water from a region with much higher downfall than the area around the lake itself, which is often a depression of some sort.Nevertheless, constructing dams and devoting water sources to agricultural, industrial and domestic uses are the most important factors for water surface slump of this lake (Hassanzadeh et al., 2012).
These types of lakes are sensitive to temperature changes.Further, the temperature is a determining factor in closed lakes, especially depthless lakes, as it controls many phenomena like evaporation, water quality and biological process.In this lake, knowing the spatial pattern of the temperature is as essential as the temporal pattern.However, there is a close relationship between the rapid increase in Earth's average surface temperature and the occurrence of land surface subsidence, which creates a serious environmental problem in urban and agricultural areas (De Bremaecker, 1983;Sneed et al., 2002).
The objective of this study was to detect and specify the non-linear correlation of land subsidence and temperature activities around the lake from 2005 to 2009.First, in order to detect localized and quasi-periodic inconstancy of the vertical components of GPS measurements, they were transformed into time-frequency spaces using the continuous wavelet transform (CWT).In addition, the researchers made use of cross wavelet transform (XWT) and wavelet coherence (WTC) between two types of time series, namely daily temperature time series and vertical components of the continuous GPS measurements, to detect and quantify the temporal and spectral interrelationships between them.Typical results of this analysis could build up a time-frequency map that can be used to analyze the time-frequency correlations between two kinds of time series.

Land subsidence in surrounding areas of Lake Urmia
Prior to the use of GPS observations, subsidence estimates were made using the InSAR observations and first order class levelling networks in the north-west of Lake Urmia.In order to study the temporal behavior of the deformation in high spatial resolution, Sedighi et al. (2009) used InSAR observations consisting of four images acquired by ENVISAT satellite chosen from descending track 92 and ascending track 228 passes from 22 September 2003 to 24 May 2004 (eight months).They reported that during this period the subsidence rate was about 6 cm yr −1 in the north-west of the lake (Gharagheshlagh) which was consistent with the results of the first order class of repeated levelling networks in the north-west of Lake Urmia.K. Moghtased-Azar et al.: Land subsidence in Lake Urmia  However, the GPS measurements have the ability to detect surface displacements with millimeter accuracy.This technique has become an important tool over the last decade for the detection of land subsidence in developed groundwater basins.National Cartographic Center of Iran (NCC) has installed more than 100 permanent stations (since 2004) across Iran to monitor ground surface movement for the geodynamic purposes (Tavakoli, 2007).As illustrated in Fig. 1, there are four permanent GPS stations in the surrounding areas of Lake Urmia.In order to study impermanent behavior of the deformation, the observed motion of each site in vertical direction could be written (Nikolaidis, 2002) where t i for i = 1 • • • N are the daily solution epochs in units of year, and H is the Heaviside step function.The first two terms are the site position a and linear rate b, respectively.Coefficients c and d describe the annual periodic motion, while e and f describe semi-annual motion.The next term corrects for any number n g of offsets, with magnitude g and epochs T g .The last term ν i denotes the measurement errors, assuming that E(ν) = 0. Using weighted least squares solution estimation of the unknown parameters X = [a b c d e f g] T with associated covariance matrix ĈX could be performed.Accordingly, to reveal how the different scales of the nonstationary time series change over the time, the researchers applied the wavelet analysis for performing time-scale decompositions of the deformation signals.
3 Threat of land subsidence in time-frequency space

The continuous wavelet transform (CWT)
There are two classes of the wavelet transforms: discrete wavelet transform (DWT) and the continuous wavelet transform (CWT).The former is a compact representation of the data and is particularly useful for noise reduction and data compression, whereas the latter is a suitable mathematical tool for the analysis of localized intermittent oscillations in a non-stationary time series as well as for feature extraction purposes (Ouadfeul et al., 2012).In mathematics, the CWT of a continuous, square-integrable function x(t) at a scale a > 0 and translational value b ∈ R is expressed by the following integral: where ψ(t) is a continuous function in both the time domain and the frequency domain called the mother wavelet and * represents the operation of complex conjugate.The main purpose of the mother wavelet is to provide a source function to generate the daughter wavelets which are simply the translated and scaled versions of the mother wavelet.To recover the original signal x(t), inverse continuous wavelet transform can be exploited: is the dual function of ψ(t), and the dual function should satisfy: where ψ(t) = C −1 ψ ψ(t), in which: is called the admissibility constant and ψ is the Fourier transform of ψ.For a successful inverse transform, the admissibility constant form has to satisfy the admissibility condition: 0 < C ψ < +∞.It is possible to show that the admissibility condition implies that ψ(0) = 0, so that a wavelet must integrate to zero (Keller, 2004).Based on CWT, the wavelet  power of a time series x(t) at the time scale space is called the scalogram and is simply defined as the squared modulus of W x (a, b).
The choice of a particular wavelet mother function can influence the time, the scale, and the frequency resolution of the time series decomposition (Mallat, 1998;Torrence and Compo, 1998).The Morlet wavelet provides a good balance between scale (frequencies) and time localizations (Grinsted et al., 2004;Mi et al., 2005).It is one of the best mother functions in terms of reproducing the frequency decomposition of the signal (Kirby and Swain, 2004).It has often been successfully used in the study of environmental variables (Torrence and Compo, 1998;Torrence and Webster, 1999;Maraun et al., 2007) and ecological variables (Ménard et al., 2007;Rouyer et al., 2008).Accordingly, Morlet wavelet was selected to be applied in the present study.Namely, the Morlet wavelet is defined as This wavelet is the product of a complex sinusoidal e −i2ω 0 t by a Gaussian envelope e −t 2 /2 where ω 0 is the central angular frequency of the wavelet.The term π −π/4 is a normalization factor to ensure unit variance.Moreover.the relation between frequencies and wavelet scales is given by 1/f = (4π/a)/(ω 0 + 2 + ω 2 0 ) (Cazelles et al., 2007).When ω 0 ≈ 2π, the wavelet scale a is inversely related to the frequency, i.e. a = 1/f .This greatly simplifies the interpretation of the wavelet analysis and one can replace, on all equations, the scale a by the period 1/f (Cazelles et al., 2007).The vertical components of GPS time series are very noisy.The values of the wavelet transform are generally corrupted as the wavelet approaches the edges of the time series, creating a boundary effect.Further, the affected region increases in extent as the scale (period) parameter increases.This region is known as the cone of influence (COI) and the spectral information within this cone is likely to be less accurate.Therefore, the introduction of a COI is suggested in which the transform suffers from these edge effects.The COI is defined so that the wavelet power for a discontinuity at the edges decreases by a factor e −2 and ensures that the edge effects are negligible beyond this point (Torrence and Compo, 1998).We consider only the information above the COI in order to avoid spurious features caused by the wavelet method.

Cross wavelet transform (XWT) and the wavelet coherence (WTC)
While CWT is a general tool for analyzing localized intermittent oscillations in a time series, in many applications, it is desirable to quantify statistical relationships between two non-stationary signals.In Fourier analysis, the coherence is used to determine the association between two signals, x(t) and y(t).The coherence function is a direct measurement of the correlation between the spectra of two time series (Chatfield, 1989).To quantify the relationships between two nonstationary signals, the following quantities can be computed: the XWT and the WTC.
The CWT of two time series x(t) and y(t) at scale a and position b are W x (a, b) and W y (a, b), respectively.Their product, namely XWT, is defined as For real-valued time series, the XWT provides a mean to indicate the coincident events over frequency, for each time in the signals x(t) and y(t).The WTC is defined as the crossspectrum normalized by (Grinsted et al., 2004) where S denotes a smoothing operator in time-scale space.
The smoothing could be obtained by a convolution with a constant-length window function in time-scale axis: denotes a smoothing operation with the weight function U δ, (τ, s) that satisfies U δ, (τ, s)dτ ds = 1 (Chatfield, 1989).The values of R x,y (a, b) are thus bounded by 0 ≤ R x,y (a, b) ≤ 1.The WTC is equal to 1 when there is a perfect A unimodal distribution of the phase difference (for the chosen range of scales or periods) indicates that there is a preferred value of ψ x,y (a, b) and thus a statistical tendency for the two time series to be phase locked.Conversely, the lack of association between the phase of x(t) and y(t) is characterized by a broad and uniform distribution (Cazelles et al., 2007).

Numerical analysis
Each vertical position time series was cleaned by excluding data based on the robust interquartile range (IQR) statistic, defined as the difference between its 75th and 25th percentiles (Nikolaidis et al., 2001).All known offsets from the GPS time series (according to the headers of time series files given by NCC) are removed by the Heaviside step function (see the top sides of Figs.2-5).
The most significant subsidence rate was observed up to −120.06 ± 1.97 mm yr −1 in the north-west of the lake (Gharagheshlagh).Moreover, the subsidence value in the western part of the region (Urmia) was −70.73 ± 0.77 mm yr −1 , in the south-east of the region (Bonab) was −11.17 ± 1.34 mm yr −1 and in the north-east of the region (Tasuj) was −4.65 ± 0.72 mm yr −1 .
Frequently, the continuous GPS measurements contain gaps and irregular sampling intervals originating from failures in the measuring equipment or the upgrade of a receiver.However, since analyzing the non-stationary data using time-frequency analysis requires equally spaced values, the gaps in the used data sets are filled using an appreciate interpolation algorithm such as multi-layer feed-forward backpropagation neural network method before the data analysis.
This method can adequately represent any arbitrary nonlinear function when a properly trained neural network is used (Moghtased-Azar and Zaletnyik, 2008).In addition, useful relationships among different inputs and outputs can be clarified by using this method which is commonly used for training the neural networks in many applications.The performance of this algorithm is reported to be satisfactory in  the interpolation and prediction of the values in time series (Jwo and Chin, 2002;Schuh et al., 2002;Erol, 2011).
So, the prepared time series of vertical components of four GPS stations in the study area were analyzed using CWT techniques.The CWT power spectra of the time series of the vertical components (associated with four GPS stations in the study area) are displayed at the bottom of Figs.2-5.Color scales indicate whether the periodicities have strong (red) or weak (blue) power and the black line contours show 95 % level of confidences for the existence of periodicities for certain scale-periods and time.This illustrations show the high positive values of the wavelet power spectrums at the time positions of A comparison among Figs.2, 3, 4, 5 and 6 reveals similar patterns in the high period bands.However, in order to search for statistically relevant correlation between the periodicities found in GPS data and daily temperature data, XWT analysis was carried out of pairs of GPS data and daily temperature time series.
Typical result of this analysis is shown in Fig. 7  Additionally, the WTC of the GPS time series and temperature time series is shown in Fig. 8. Compared with the XWT, large sections stand out as significant areas.The WTC finds regions in time frequency space where the two time series co-vary (but does not necessarily have high power).Consequently, water level variations of Lake Urmia (influenced by drought conditions in the region) in terms of simultaneous time intervals of GPS and daily temperature measurements was illustrated by the satellite altimetry data.Namely, the preferred time interval for the altimetry data was considered from October 2005 to September 2010.Figure 9 shows relative Lake Urmia height variations computed by Jason-1 and Jason-2/OSTM altimetry while Jason-1 data was from February 2002 to January 2009 and Jason-2/OSTM data was from July 2008 to August 2011 (Crétaux et al., 2011).All classical corrections (orbit, ionospheric and tropospheric corrections, polar and solid Earth tides) were applied.
Further, a specified time interval, namely from October 2005 to September 2010, was separated into three subintervals and afterwards we have calculated the linear trends of water level variations over the individual intervals.Namely,

Conclusions
The objective of this study was to detect and specify the nonlinear correlation of land subsidence and temperature activities in the region from 2005 to 2009 around Lake Urmia.The time-frequency maps of GPS and temperature time series demonstrated similarity between the portrayed patterns in the high period bands.In addition, in order to search for statistically relevant correlation between the periodicities found in GPS data and daily temperature data, XWT analysis was ).This event was detected by XWT as a critical interval to be holding the strong correlation between the land subsidence phenomena and surface temperature.

Fig. 1 .
Fig. 1.Map showing location of study area and of four continuous GPS stations around it.

Fig. 1 .
Fig. 1.Map showing location of study area and of four continuous GPS stations around it.

Fig. 2 .Fig. 3 .
Fig. 2. The scalograms of the GPS time series of Bonab.This illustration shows the high positive values of the wavelet power spe the time positions of: 72-490 days band in the period from 7 Apr 2006 to 2 Aug 2010, 77-210 days band in the period from 26 M to 23 Feb 2009 and 560-1024 days band in the period from 11 Jul 2007 to 12 Dec 2009.The black contour line shows the 95% co interval.

Fig. 2 .
Fig. 2. The scalograms of the GPS time series of Bonab.This illustration shows the high positive values of the wavelet power spectrum at the time positions of 72-490 days band in the period from 7 April 2006 to 2 August 2010, 77-210 days band in the period from 26 May 2007 to 23 February 2009 and 560-1024 days band in the period from 11 July 2007 to 12 December 2009.The black contour line shows the 95 % confidence interval.

Fig. 2 .Fig. 3 .
Fig. 2. The scalograms of the GPS time series of Bonab.This illustration shows the high positive values of the wavelet power spec the time positions of: 72-490 days band in the period from 7 Apr 2006 to 2 Aug 2010, 77-210 days band in the period from 26 Ma to 23 Feb 2009 and 560-1024 days band in the period from 11 Jul 2007 to 12 Dec 2009.The black contour line shows the 95% con interval.

Fig. 3 .
Fig. 3.The scalograms of the GPS time series of Gharagheshlagh.This illustration shows the high positive values of the wavelet power spectrum at the time positions of 144-235 days band in the period from 20 February 2008 to 13 September 2008, 148-1024 days band in the period from 17 January 2008 to 9 April 2010, 35-88 days band in the period from 1 May 2009 to 20 January 2009 and 98-134 days band in the period from 25 April 2009 to 15 July 2009.The black contour line shows the 95 % confidence interval.
the GPS time series of Urmia.This illustration shows the high positive values of the wavelet power spectrum at ositions of: 183-1024 days band in the period from 15 Oct 2007 to 30 Apr 2010.The black contour line shows the 95% confidence 29 Dec 2005 29 Oct 2006 29 Aug 2007 28 Jun 2008 28 Apr 2009 26 Feb 2010 the GPS time series of Tasuj.This illustration shows the high positive values of the wavelet power spectrum at the tions of: 85-645 days band in the period from 23 Oct 2006 to 18 Dec 2007, 134-610 days band in the period from 14 Jan 2008 to 09, 228-445 days band in the period from 30 Oct 2009 to 18 Aug 2010, 167-260 days band from 22 Apr 2010 to 23 Aug 2010 and 8 days band in the period from 15 Sep 2007 to 1 Jan 2010.The black contour line shows the 95% confidence interval.

Fig. 4 .
Fig. 4. The scalograms of the GPS time series of Urmia.This illustration shows the high positive values of the wavelet power spectrum at the time positions of 183-1024 days band in the period from 15 October 2007 to 30 April 2010.The black contour line shows the 95 % confidence interval.

Fig. 4 .Fig. 5 .
Fig. 4. The scalograms of the GPS time series of Urmia.This illustration shows the high positive values of the wavelet power spe the time positions of: 183-1024 days band in the period from 15 Oct 2007 to 30 Apr 2010.The black contour line shows the 95% co interval.

Fig. 5 .
Fig. 5.The scalograms of the GPS time series of Tasuj.This illustration shows the high positive values of the wavelet power spectrum at the time positions of 85-645 days band in the period from 23 October 2006 to 18 December 2007, 134-610 days band in the period from 14 January 2008 to 19 February 2009, 228-445 days band in the period from 30 October 2009 to 18 August 2010, 167-260 days band from 22 April 2010 to 23 August 2010 and 1040-2048 days band in the period from 15 September 2007 to 1 January 2010.The black contour line shows the 95 % confidence interval.
i. Bonab: 72-490 days band in the period from 7 April 2006 to 2 August 2010, 77-210 days band in the period from 26 May 2007 to 23 February 2009 and 560-1024 days band in the period from 11 July 2007 to 12 December 2009.ii.Gharagheshlagh: 144-235 days band in the period from 20 February 2008 to 13 September 2008, 148-1024 days band in the period from 17 January 2008 to 9 April 2010, 35-88 days band in the period from 1 May 2009 to 20 January 2009 and 98-134 days band in the period from 25 April 2009 to 15 July 2009.iii.Urmia: 183-1024 days band in the period from 15 October 2007 to 30 April 2010.iv.Tasuj: 85-645 days band in the period from 23 October 2006 to 18 December 2007, 134-610 days band in the period from 14 January 2008 to 19 February 2009, 228-445 days band in the period from 30 October 2009 to 18 August 2010, 167-260 days band from 22 April 2010 to 23 August 2010 and 1040-2048 days band in the period from 15 September 2007 to 1 January 2010.Moreover, Fig. 6 illustrates time series of daily temperature values of four GPS stations in the study area (in which thermometers are located near the GPS receivers).This illustration shows that all stations have common high positive values of the wavelet power spectrums at the 285-437 days band in the period from 26 January 2007 to 22 November 2008.
for the GPS time series and temperature time series.In all panels, the curved line in the bottom indicates the cone of influence under which the signals must be discarded.The color level bar indicates the intensity of each signal.The contour black line determines a 95 % level of confidence for the existence of a periodicity for certain scale/period and time.The arrows in each plot indicate the phase difference between the components of the two time series.The right arrow indicates an in-phase relation while the left arrow means anti-phase.The results illustrate significant common powers in: i. Bonab: 230-256 days band from 8 February 2007 to 9 March 2008 and 157-218 days band from 18 August 2007 to 1 March 2009.ii.Gharagheshlagh: 90-235 days band from 15 October 2007 to 19 June 2009.iii.Urmia: 180-256 days band from 16 September 2007 to 4 February 2009.iv.Tasuj: 62-67 days band from 7 January 2007 to 10 April 2007, 109-134 days band from 27 October 2006 to 17 March 2007 and 118-248 days band from 14 May 2007 to 19 March 2009.Clearly, significant common patterns in the XWT of four pairs of time series are illustrated in the high period bands from 180-218 days band from September 2007 to February 2009.

Fig. 6 .
Fig. 6.The scalograms of the temperature time series.This illustration shows that all stations have common high positive values of the wavelet power spectrums at the 285-437 days band in the period from 26 January 2007 to 22 November 2008.The curve delimits the cone of influence (COI), which indicates that the region was not influenced by edge effects.The black contour line shows the 95 % confidence interval.

Fig. 8 .
Fig. 8.The WTC between the GPS time series (associated with vertical components) and temperature time series.The WTC finds regio in time frequency space where the two time series co-vary (but do not necessarily have high power).The black contour line shows the 95 confidence interval.

Fig. 8 .
Fig. 8.The WTC between the GPS time series (associated with vertical components) and temperature time series.The WTC finds regions in time frequency space where the two time series co-vary (but do not necessarily have high power).The black contour line shows the 95 % confidence interval.

Fig. 9 .
Fig. 9. Relative Lake Urmia height variations computed by Jason-1 and Jason-2/OSTM altimetry in which Jason-1 data is from February 2002 to January 2009 and Jason-2/OSTM data is from July 2008 to August 2011 (available on the web site HYDROWEB: http: //www.LEGOS.obs-mip.fr/soa/hydrologie/HYDROWEB).The linear trends of water level variations over the individual intervals are calculated: before Sep 2007 with a rate −23.08 cm yr −1 , from September 2007 to February 2009 with a rate −70.5 cm yr −1 and after February 2009 with a rate 21.38 cm yr −1 .