Natural Hazards and Earth System Sciences Advanced interpretation of subsidence in Murcia ( SE Spain ) using A-DInSAR data – modelling and validation

Subsidence is a natural hazard that affects wide areas in the world causing important economic costs annually. This phenomenon has occurred in the metropolitan area of Murcia City (SE Spain) as a result of groundwater overexploitation. In this work aquifer system subsidence is investigated using an advanced differential SAR interferometry remote sensing technique (A-DInSAR) called Stable Point Network (SPN). The SPN derived displacement results, mainly the velocity displacement maps and the time series of the displacement, reveal that in the period 2004–2008 the rate of subsidence in Murcia metropolitan area doubled with respect to the previous period from 1995 to 2005. The acceleration of the deformation phenomenon is explained by the drought period started in 2006. The comparison of the temporal evolution of the displacements measured with the extensometers and the SPN technique shows an average absolute error of 3.9 ±3.8 mm. Finally, results from a finite element model developed to simulate the recorded time history subsidence from known water table height changes compares well with the SPN displacement time series estimations. This result demonstrates the potential of A-DInSAR techniques to validate subsidence prediction models as an alternative to using instrumental ground based techniques for validation.


Introduction
Ground subsidence is a natural hazard that affects wide areas causing important economic damages and a high social Correspondence to: G. Herrera (g.herrera@igme.es)alarm.It may be due to several causes as deep material dissolution, excavation of ground tunnels or mining galleries, fluid withdrawal (petroleum, water, gas), deep erosion (piping), lateral soil creep, compaction of soil materials or tectonic activity.All these causes are manifested at terrain surface as vertical deformations that can vary from a few millimetres to several meters during periods that vary from minutes to years.In this case subsidence due to aquifer system consolidation induced by ground water overexploitation is investigated.This phenomenon involves the settlement of the ground surface and affects wide areas.Human structures built on these areas must withstand this subsidence, and damage occurs when differential settlements cannot be accommodated by their foundations.It is estimated that there are over 150 cities in the world with serious problems of subsidence due to excessive groundwater withdrawal (Hu et al., 2004).For instance, well-known examples of subsidence include the Po Valley (Italy), Mexico DC, Antelope, Santa Clara and San Joaquin Valleys (USA), Bangkok (Thailand) and many other areas in the world.Subsidence has occurred in the metropolitan area of Murcia City (SE Spain) as a result of excessive pumping of groundwater.Martínez et al. (2004) documented damages in buildings and other structures with an estimated cost of more than 50 million euros that generated a significant social impact.
Differential Interferometry with Synthetic Aperture Radar (DInSAR) has become over the last decade an important remote sensing tool for the estimation of the temporal evolution of ground surface displacements.This technique is able to analyze wide areas periodically and at a low cost, monitoring deformation of the ground surface at millimetre level.The standard DInSAR technique, which is based on a single interferogram generated from a pair of SAR images (Rosen et Published by Copernicus Publications on behalf of the European Geosciences Union. al., 2000), has been improved in the last years by several Advanced DInSAR techniques (A-DInSAR), which are based on the processing of multiple interferograms derived from a large set of SAR images (Ferretti et al., 2000;Berardino et al., 2002;Mora et al., 2003;Arnaud et al., 2003;Werner et al., 2003;Hooper et al., 2004).An interesting standard DIn-SAR technique application to subsidence analysis was described by Galloway et al. (1998), Hoffman (2003) and Galloway and Hoffmann (2007).Advanced A-DInSAR applications are found in Ferretti et al. (2004), Zerbini et al. (2007) and Bell et al. (2008).
The temporal evolution of the subsidence in Murcia city has been monitored by an extensometer network since 2001 (Mulas et al., 2001(Mulas et al., , 2005;;Peral et al., 2002), and since 1975groundwater table variations (Aragón et al., 2004) have been recorded through a spatially dense piezometric network.This phenomenon has been also analysed through the application of different Advanced DInSAR techniques in previous works (Mulas et al., 2004;Tomás et al., 2005;Herrera et al., 2008).In this case 129 SAR images of the study area acquired from descending orbits from July 1995 to December 2008 are analysed for two different time periods; 1995-2005 and 2004-2008 with the Stable Point Network (SPN) technique (Arnaud et al., 2003).
The SPN derived displacement results, mainly velocity displacement maps and time series of displacement, are analysed in detail and compared with previous available information on the subsidence phenomenon in Murcia.The relationship between the temporal evolution of the displacement and the water table in representative sites of the metropolitan area of Murcia is analysed.Then the displacement velocity maps of the whole basin are compared with the thickness map of compressible deposits.Subsequently, the temporal evolution of the displacement is compared with the extensometer measurements of the deformation to estimate the error of the SPN technique.Finally, a finite element model is used to simulate the recorded subsidence using known water table height changes.Results are compared and validated with the extensometer and SPN displacement time series estimations.

Description of the study area
This work is focused on a section of the Segura River Valley (SE Spain) that corresponds to the metropolitan area of the city of Murcia, which is the most populated city of the basin with more than 400 000 inhabitants.This area has developed by occupying the flood plain of the Segura River, also known as the Vega Media, covering a surface of 206 km 2 , with approximately 50% being cultivated land (Gumiel et al., 2001).The Segura River runs from SW to NE along the valley with height changes that vary from 60 m o.s.l in Alcantarilla to 30 m o.s.l in Beniel (Fig. 1b) From a geological point of view the study area is located in the oriental sector of the Betic Cordillera.A compressive stress field has acted since the Upper Miocene in this sector and has led to the development of a basin bounded by active faults, the Lorca-Alhama, to the north, and the Carrascoy-Bajo Segura, to the south (Montenat et al., 1990).The sedimentary record of the basin has been deformed simultaneously to its deposition, creating a broad syncline in which progressively younger sediments have been deposited.
The basement is made up by old (Permian and Triassic), deformed materials corresponding to the Internal Zones of the Betic Cordillera.These materials also crop out along the border of the Segura River Valley.The basin fill comprises Upper Miocene to Quaternary sediment deposits that can be divided into three units (Fig. 1b).Older materials (Upper Miocene) composed mainly by a thick sequence (more than 600 m) of marls (Cerón and Pulido-Bosch, 1996;Mulas et al., 2003).Above them, Pliocene-Quaternary rocks comprise marls and clays interbedded with several levels of gravels and sands, deposited in a continental environment that can reach 250 m at some places (Aragón et al., 2004).At the ground surface, recent continental (meander, channel, oxbow lakes, flood plain, alluvial fans, etc.) sediments are found.Silts and clays are abundant in flood plain and oxbow deposits, while sand is common in channel areas and in the alluvial fans formed in the reliefs that surround the valley.Thickness of recent sediments varies between 3 and 30 m (Rodríguez-Jurado et al., 2000).Additionally, at some places anthropic deposits can be found on top.
The aquifer system is constituted by two units, a superficial aquifer and a deep aquifer (Cerón and Pulido-Bosch, 1996;Aragón et al., 2004).The superficial aquifer comprises the first 5 to 30 m of recent clay, silt and sands facies (Fig. 1c).The deep aquifer system, located under the previously described aquifer, is composed of a sequence of gravels and sands alternating with marls, clay and silts (Pleistocene to Pliocene materials).In the upper part of this aquifer, just beneath the superficial aquifer, there is a gravel layer, 10 to 30 m thick that has been intensively exploited since 90 decade as an important water resource.The materials that constitute the superficial aquifer are the most compressible ones.They constitute a layer of medium to soft sediments, susceptible to suffer consolidation due to variations (increase) of effective stresses acting on them.In the following sections they will be referred as the compressible deposits.On the contrary, the underlying gravels and sands are more rigid and represent the geotechnical substratum of the zone, used as the support level for deep foundations.
Our subsidence model suggests that when water is pumped from the upper gravels of the deep aquifer, a gradient is created resulting in the lowering of the water table.Consequently, pore pressure in the superficial aquifer declines and the aquifer undergoes a consolidation process.In the past 20 years two drought periods have occurred: 1992-1995 and 2006-2009.In the first period a lowering of the piezometric level of 8 m approximately caused a maximum estimated ground settlement of 8 cm (Martínez et al., 2004).).The direct cost was estimated in 50 million euros (Martínez et al., 2004).The local authorities requested a detailed study of the aquifer system including a considerable number of boreholes and laboratory and "in situ" tests (Mulas et al., 2000).

Data processing and results
The SPN algorithm has been applied to images acquired by the European Space Agency (ESA) ERS-1/2 and Envisat ASAR sensors covering two periods July 1995-December 2005 and January 2004-December 2008.A similar crop of about 20 km×8 km was selected from the 100 km×100 km acquired SAR images, corresponding to the Vega Media of the Segura River (Fig. 2).ach interferometric pair has been selected with a perpendicular spatial baseline smaller than 800 m, a temporal baseline shorter than 6 years in the case of 1995-2005 and 3 years in the case of 2004-2008, and a relative Doppler centroid difference below 400 Hz.The DEM of the Shuttle Radar Topography Mission (SRTM) has been used.The pixel selection for the estimation of displacements was based on a combination of several quality parameters including low amplitude standard deviation and high model coherence.
The maps of the cumulative displacement in August 2005 and December 2008, estimated along the line of sight (LOS), obtained by the SPN technique for both periods are presented in Fig. 2. The deformation maps have been geocoded and superimposed on an aerial image, showing 9 colours that represent intervals of LOS displacement.Note that reliable information is available only in urban areas and rock outcrops.In this sense, displacements are observed mainly on the urban and sub-urban areas of the central part of the valley, corresponding with the ground surface beyond unconsolidated sediments.In the period 1995-2005 two main areas of deformation have been detected.The area covered by the extensometers around the urban area of Murcia is affected by total displacements that vary between −1 and −4 cm (Fig. 2).To the northeast a subsidence area (H5 in Fig. 2) shows total displacements between −3 cm and −6 cm.In the period 2004-2008 a generalised subsidence in the valley is observed of greater magnitude than in the former period.In fact, the areas indicated above also underwent the greatest deformation of the study area in the second period.In the area covered by the extensometer network subsidence varies between −3 and −6 cm, and in the northeast cumulative subsidence varies between −4 and −12 cm.

Results interpretation
The displacements measured with the SPN technique has been geocoded and compared with the temporal evolution of piezometric levels and the thickness of the compressible layer, which are the key parameters that control subsidence due to aquitard consolidation of the metropolitan area of Murcia.

Comparison with piezometric level
In general, piezometric level is closely related with annual precipitation because rain infiltration and irrigation are the most important sources of recharge of the aquifer.The historical evolution of the piezometric level shows a stable situation, with only small annual oscillations lower than 2-4 m (H1 in Fig. 3).During the drought periods water extractions in authorized and illegal wells from deep aquifer increase and infiltration decrease.The consequent piezometric level decline whose effects have been observed during the historical droughts of 1980-1983, 1993-1995and 2006-2008 (H1 in Fig. 3 (H1 in Fig. 3).The drought in the 80's derived in piezometric falls higher than 5 m during 1983 in the Segura Valley.However, there is no information about ground subsidence during this period.The 1993-1995 drought caused a generalized lowering of 8 m (Mulas et al., 2000), with maximum values higher than 15 m.The 2005-2008 drought period caused a water table lowering in the range between 6 and 13 m.
The temporal evolution of A-DInSAR measured ground deformations shows a good correlation with piezometric level changes in the deep aquifer.The analysis of these results has identified two main patterns of ground subsidence during this period for the study zone.Figure 3 shows the relationship between A-DInSAR measured surface deformation and piezometric level for representative sites of the metropolitan area of Murcia.Piezometers H2, H3 and H4 are located in the urban area of Murcia, on recent deformable sediments 12 to 22 m thick.Piezometer H1 is in an G. Herrera   area southwest of the city where wells exploit deep aquifers during drought periods.In this case clay layers alternate with sands and gravel beds.Finally piezometer H5 and H6 are located to the northeast following the Segura River on deformable silts and clays that reach up to 28 m thick.Ground deformation occurred between 1993 and 1995 could not be detected with the SPN processing because the first reliable SAR image in this area dates from July 1995.The piezometric level recuperated in 1997 and it has remained almost constant, showing small seasonal fluctuations of a few meters and a general descending trend until 2005.The A-DInSAR displacements measured in this period show slow velocities in agreement with the measured small piezometric decrease.Between January 2005 and December 2008 the generalised lowering of the water table coincides with an increase of the A-DInSAR estimated deformation velocities.The piezometers located in the urban area of Murcia (H2, H3 and H4) show a piezometric decrease from 1.5 to 7 m and a subsidence rate from 2 to 5 mm/year for each of the above mentioned periods.On the other hand, the piezometers H5 and H6, located towards the northeast from the urban area, measured a water table decline up to 8.2 m and the corresponding displacement velocity increase was as much as 1 cm/year during this period.In this sector the velocity of the deformation almost doubled with respect to the urban area, probably due to the thicker compressible layer observed in this area that varies from 15 to 28 m.A different behaviour is not observed in piezometer H1 where the estimated displacement velocity was 4.7 mm/year, despite a 13 m water table decline (double than in the rest of piezometers) and the presence of a greater than 20 m thick compressible layer (Fig. 4a).Nevertheless the new available geotechnical boreholes in this area shows that the thickness of the compressible layer around piezometer H1 is similar to that observed in the urban area of Murcia (between 10 and 20 m deep) instead of greater than 20 m as observed in Fig. 4a, being in this case clay materials more consolidated.

Comparison with the thickness map of compressible deposits
A thickness map of compressible strata was elaborated by Mulas et al. (2000) interpolating the thickness data available from geotechnical boreholes (Fig. 4a).This map allows analyzing the relationship between the SPN velocity maps and the cumulative thickness of sediments susceptible to suffer consolidation.In Fig. 4b and c the displacement velocity maps are superposed to the thickness map of compressible deposits.By comparing these maps it is observed that the greatest velocity is found at the northeast where the thickest compressible deposits are found.
In order to perform a quantitative comparison for both periods, the mean and standard deviation of the velocity of all the pixels (SPs) located within the different compressible thickness areas have been calculated.The bars in Fig. 5 represent the number of SPs found in each of these compressible areas and the interpolated points represent the mean and standard deviation values of the velocity of all the pixels found within each of the compressible thickness areas.For both periods it is observed that the thicker the compressible layer the greater the estimated displacement velocity.

Comparison with the extensometer time series
The SPN accuracy has been assessed by comparing retrieved deformations against measurements from extensometer boreholes.The basic idea is to compare the displacement time series achieved from the SAR data analysis for both periods with those available from the extensometer network measurements projected along the Line of Sight (LOS).In order to enable the comparison LOS-projected extensometer time series have been interpolated via a linear regression within the interval common to the SAR data, i.e. from 2001 to 2007.
The comparison with the extensometers has been done calculating two errors: (1) the mean and standard deviation of the absolute difference between SPs and LOS-projected extensometer deformation time series, (2) the mean and standard deviation of the difference between SPs and LOSprojected extensometer deformation time series.Figure 6 shows the location of the selected SPs around each extensometer.The temporal evolution of the displacements measured with eight representative extensometers and estimated with the SPN technique is plotted in Fig. 7.A good agreement can be observed between the ground based measurements and the SAR based estimations of deformation.
Comparing the complete extensometer network with the SPN displacement time series, the mean (µ and standard deviation (σ values of the differences are µ=t1.4mm and σ =4.9 mm, and the mean and standard deviation values of the absolute differences are µ=3.9mm and σ =3.8 mm.The most recent validation experiments available in the scientific literature (Strozzi et al., 2001;Hanssen, 2003;Colesanti et al., 2003;Casu et al,, 2006;Crosetto et al., 2008;Herrera, 2008;Tomás 2009) show that the error of the deformation time series, estimated with different A-DInSAR techniques, compared with ground truth measurements is within the ±6.9 mm interval.These values show the realistic numbers achieved in this analysis (Table 1).
Nevertheless, some extensometers (ei1, ei4 and v3) have recorded smaller displacement velocities than those estimated with the SPN technique since October 2004, coinciding with the start of the most recent drought period.A possible explanation of this difference can be that the SPN technique is measuring ground deformations whereas the extensometer network is measuring deformations only in the first 20 m below the surface.Deformation of compressible layers deeper than 20 m would not be detected by the extensometer network (Tomás, 2009).

Subsidence forecasting model using SPN data
A finite element model is developed to simulate the recorded time history subsidence considering known water  (Universidad Politécnica de Madrid) has been used.This code simulates the time dependent coupling behavior between the soil and the interstitial fluids solving the Biot (1941) equations.

Model description
The metropolitan area of Murcia covers several square kilometers, being the considered depth of the study around 30 m, with too few data to simulate a full three dimensional  analysis.For this reason several one-dimensional analyses have been preferred.The considered columns of soil coincide with geotechnical boreholes drilled near the extensometers where materials are well defined (Fig. 8).The water table height given by piezometer H-1 begins in 1975, whereas measurements in piezometers H-2, H-3 and H-4 start in 1994.
In order to take into account water level changes of these piezometers since 1975, data of piezometer H-1 from 1975 to 1994 has been used for the rest of piezometers.
The different soil columns have been discretized using a mixed eight node quadrilateral element with reduced integration (8 nodes for displacements, 4 nodes for pore pressure) that permits to simulate the coupling between the solid skeleton and the interstitial water.We assume that the col-umn is fixed at the base and that the horizontal displacement is zero (one dimensional analysis).Concerning the hydraulic boundary conditions, fluctuations of water table height measured by the piezometers have been imposed prescribing the pore pressure in the base of the column.The bottom of the column is always placed in the gravel materials where changes in pore pressure are almost immediate.As no additional information concerning suction on the surface is available other boundaries are supposed to be impermeable.
The constitutive relation is a key point because the model will not be able to reproduce features that are not contemplated in this relation.In this case an elastic constitutive relation has been chosen for three main reasons: (1) Mohr-Coulomb criteria does not permit to compute plastic Concerning the hydraulic constitutive relations, permeability is assumed to be constant and the influence of the degree of saturation (Sr) has been studied assuming three different retention curves using the simple relation proposed by (Lloret and Alonso, 1985).The first one corresponds to the saturated case with Sr1=100%, the second one corresponds to S r =1−m tanh(n s) and the third one corresponds to Sr3=1−0.8 tanh(20 s), being Sr and s the degree of saturation and the suction respectively.Therefore the density of the material above the water table height will vary depending on the selected retention curve.
As it is shown in the columns of Fig. 8 three main materials have been considered for the aquifer system simulation.The In some cases sand layers alternate with clay units.The material parameters adopted in the calculations (Table 2) were taken from Mulas et al. (2000).It can be seen that vertical permeability coefficients are quite small but in agreement with those given by Neuzil (1994).Therefore we expect a very slow consolidation process.
A first static analysis was computed in order to obtain an initial effective stress state and an initial pore pressure G. Herrera et   distribution.This calculation corresponds to the hydrostatic equilibrium for the water table height measured the 12 February 1975.Before calculating the temporal evolution of the displacement we have analyzed the influence of the retention curve defined for the partially saturated zone on the subsidence results.In this case a time dependent consolidation analysis showed that the retention curve has a minor influence on the computed vertical settlement.Consequently the saturated case has been adopted for all the computations.The model simulates the slow dissipation of the excess pore pressure in the clay material due to the very small value of vertical permeability.This process was modeled using a time dependent analysis with a time step equal to one day from the 12 February 1975 until the last available record of water table height.

Model results
The temporal evolution of the displacement computed with the model and projected in the LOS is compared with the displacement time series measured with the SPN technique for the SPs near extensometer locations (Fig. 9). Figure 8 shows the distribution of the different facies in the soil columns chosen for each of the extensometers, as well as the name of the geotechnical and piezometric borehole used for each of them.The planimetric location of the extensometers with respect to the above mentioned boreholes as well as the location of the SPs used for the comparison is shown in Fig. 6.These plots demonstrate the strong similarity between the computed deformations and those measured from the SAR images dataset.Moreover two errors have been calculated: This comparison between model and SPN displacement time series shows that the mean values of the difference of the time series for all the extensometers location are µ=−0.5 mm (average) and σ =6.0 mm (standard deviation), and the average values of the absolute difference are µ=5.5 mm and σ =4.7 mm (Table 3).The fourth and fifth columns from the left in Table 3 show that the average maximum and minimum values of the difference of the SP and the model time series for the extensometer locations are µ=16.4mm and µ=−14.4mm.Values for individual extensometric locations indicate that there are significant differences between the compared time series.These differences are concentrated in the 1995-1997 period (Fig. 9) where the SPN technique was unable to monitor the uplift associated with the recuperation of the water table computed by the model until 1997.The SPN could not estimate the subsidence experienced in the drought period (1992)(1993)(1994)(1995) due to the lack of reliable SAR images before 1995.Nevertheless the deformations ocurring in the most recent drought period (2006)(2007)(2008)(2009) were accurately estimated with the SPN technique.

Conclusions
The SPN derived displacement results, mainly the velocity displacement maps and the time series of the displacement, reveal that in the period 2004 -2008 the rate of subsidence in Murcia metropolitan area doubled with respect to the previous period from 1995 to 2005.The acceleration of the deformation phenomenon is explained by the drought period started in 2006.An average piezometric lowering of 7 m has been correlated with a displacement velocity increase from 2 to 5 mm/year between the first and second monitoring periods.A velocity increase up to 1 cm/year has been observed in the north-western part of the metropolitan area where the thickest compressible deposits are found.
The comparison of the SPN velocity maps for both analysed periods with the thickness of the compressible deposits confirms a direct relationship between both variables.The comparison of the temporal evolution of the displacements measured with the extensometers and the SPN deformation time series validates the DInSAR results and provides and assessment of accuracy of the remote sensing technique.The average (µ and standard deviation (σ values computed for the absolute difference between SPN deformations and extensometers are µ=3.9mm and σ =3.8 mm.
Finally, results from a finite element model developed to simulate the recorded time history subsidence from known water table height changes compares well with the SPN displacement time series estimations, showing that the average values of the absolute difference between the numerical model and the SPN time series of displacement for all the extensometers locations are µ=5.5 mm and σ =4.7 mm.This result demonstrates the potential of A-DInSAR techniques to validate subsidence prediction models as an alternative to using instrumental ground based techniques for validation.

Fig. 1 .
Fig. 1.Situation and geology of the Vega Media of the Segura River.

Fig. 2 .
Fig. 2. Top: Total deformation estimated with SPN technique between July 1995 and December 2005.Bottom: Total deformation estimated with SPN technique between January 2004 and December 2008.Labels have been placed in those piezometers that will be analysed in the next sections.

Fig. 3 .
Fig.3.Temporal evolution of the water table compared with the temporal evolution of the deformation estimated for the periods1995-2005  and 2004-2008.

Fig. 4 .Fig. 5 .
Fig. 4. (a) Thickness map of compressible deposits (modified from Mulas et al., 2000).Measured LOS displacement rate between 1995 and 2005 (b), and between 2004 and 2008 (c) superposed over the thickness map of compressible deposits.Labels have been placed in those extensometers (a) that are analysed in the next sections.

Fig. 6 .
Fig. 6.Location of the selected SPs around each extensometer.The white arrow indicates the distance to the closest piezometer.

Fig. 7 .
Fig. 7. Time series of the LOS displacements estimated with the SPN technique and with the LOS-projected measurements from the extensometers.

Fig. 9 .
Fig. 9. LOS projected time series of the deformation estimated with the model and with the SPN technique.

Table 1 .
Comparison between SP's and LOS projected extensometer time series.
al.: Advanced interpretation of subsidence in Murcia (SE Spain) using A-DInSAR data

Table 3 .
Comparison between SP's and LOS projected model time series.