Subsidence activity maps derived from DInSAR data: Orihuela case study

Introduction Conclusions References


Introduction
Subsidence caused by water withdrawal is a well-known phenomenon which affects areas worldwide .Structures built in these areas must withstand the vertical (and sometimes also horizontal) ground displacements, resulting in widespread damages when they cannot support differential settlements under their foundations (Skempton and Mac-Donald, 1956;Bjerrum, 1963;Burland and Wroth, 1974;Boscardin and Cording, 1989;Namazi and Mohamad, 2013).Thus, the monitoring of the displacements affecting the structures plays an essential role on their serviceability and, exceptionally, on their safety.Classical surveying techniques (e.g.topographical techniques and extensometers) are commonly used for monitoring structures placed on subsiding areas (Poland et al., 1984;Galloway et al., 1999;Galloway and Burvey, 2011;Tomás et al., 2014).However, in recent decades, differential SAR (synthetic aperture radar) interferometry (DInSAR) has become an alternative method for measuring infrastructure displacements, showing a good capability for measuring ground surface displacements over wide areas (Galloway and Burvey, 2011;Tomás et al., 2014).
This work focuses on the utility of A-DInSAR techniques to monitor and characterize the phenomenon of subsidence both regionally and locally, adopting a geotechnical approach for the local scale.A preliminary analysis was shown in Sanabria et al. (2014).In this paper, a much deeper analysis and discussion is presented.The ground surface settlements, caused by groundwater withdrawals, in the city of Orihuela are evaluated by ERS (European Remote-Sensing Satellite)-1/2 and Envisat ASAR (Advanced Synthetic Aperture Radar) sensors covering two different periods July 1995-December 2005 and January 2004-December 2008.The A-DInSAR displacement data obtained are interpolated based on conditional Sequential Gaussian Simulation (SGS) to generate subsidence activity maps both at regional and local scales.The conditional SGS allows quantifying the spatial variability of the interpolation and provides a confidence level on the interpolation.From the interpolated maps, the serviceability limit states of 27 heritage buildings (16th-19th century) of the city of Orihuela (Fig. 1) are studied by means of geometrical-geotechnical criteria (i.e.differential settlements and angular distortions).A qualitative evaluation of the results is made with 10 buildings where damages have been reported.Finally, a comparison between the serviceability parameters obtained for the Santas Justa and Rufina Church (Iglesia de las Santas Justa y Rufina) and a detailed damage study (Tomás et al., 2012) is performed.

Description of the study area
The city of Orihuela is located in the Vega Baja of the Segura River (VBSR) (province of Alicante, SE Spain).The basin is filled by Neogene-Quaternary sediments deposited by the Segura River.The substratum of the VBSR basin (Fig. 1a) consists on Permo-Triassic rocks and Tertiary sediments that outcrop in the north and south of the basin (de Boer et al., 1982).These materials constitute the geotechnical substratum of the city of Orihuela, being the Pleistocene to Holocene sediments the most compressible ones.The spatial distribution of soft soils in the VBSR increases towards to the centre of the valley reaching a maximum thickness of up to 50 m (Delgado et al., 2000).These data have been digitized and completed with the geotechnical information derived from new boreholes located west of the city of Orihuela and interpolated with the kriging method (Matheron, 1963) to generate the soft soil thickness map (Fig. 4).Tomás et al. (2010) characterized these soft sediments as sediments with moderate to high compressibility, exhibiting compression indexes (Cc) varying from 0.07 to 0.29 and with an average value of 0.18.
From a hydrogeological point of view the study area belongs to the Guadalentín-Segura aquifer system (IGME, 1986), which is divided in two units.The first one is an unconfined shallow aquifer with low conductivity (sand, silts and clays) with the water table a few metres below the ground surface.The deep aquifer is formed by gravels, usually interbedded with marls, showing a greater hydraulic conductivity than the upper aquifer (IGME, 1986).
In the past decades, the recent subsidence phenomenon in the city of Orihuela is related to the excessive water pumping from the shallow and deeper aquifers during drought periods : 1993: -1996: and 2006: -2008: (Tomás et al., 2010)).During these periods groundwater extraction increased , in authorized and illegal wells, while groundwater recharges were reduced due to low precipitations.This groundwater withdrawal caused an important piezometric level drop (Fig. 1b).During the years 1993-1996 the piezometric level dropped about 5 m as is shown in piezometer no. 4. For years 1998-2002 the decrease was about 5-8 m.Finally, the more pronounced decline occurred during the years 2004-2008, where the piezometers 1 and 3 show a drop of 18 and 10 m, while piezometers 2 and 4 show a decline of about 8 m.The drop of piezometric levels entailed an important soil consolidation due to pore water pressure decrease (Mulas et al., 2003;Tomás et al., 2007aTomás et al., , b, 2010)), producing moderate damages to structures and infrastructures (Martínez et al., 2004) in the whole Segura River basin.Subsidence-related damages were reported in the western part of the city of Orihuela in 1995 by the local press.Furthermore, several heritage buildings affected by settlements have been repaired since the 1990 s (Louis, 2005;Maciá, 2005;Louis et al., 2012).Recently, a subsidence damage assessment of a Gothic church in the city of Orihuela, using advanced differential interferometry for the 1995-2008 period and field data, was performed by Tomás et al. (2012).

Methodology
In this section a new proposed methodology to obtain subsidence activity maps from persistent scatterer interferometry (PSI) displacement estimates is described.A more detailed description of the PSI techniques can be found in Arnaud et al. (2003) and Duro et al. (2005).Sansosti et al. (2010) and Tomás et al. (2014) also present a review of PSI applications for subsidence research.In this work, ground subsidence measurements were obtained using a PSI technique called the stable point network (SPN).The SPN software (Arnaud et al., 2003;Duro et al., 2005) uses the DIAPA-SON interferometric algorithm for all SAR data handling, e.g.co-registration task and interferogram generation.The SPN procedure generates three main products starting from a set of single-look complex (SLC) SAR images: (a) the average deformation velocity along the LOS of a single persistent scatterer (PS), (b) a map of height error, and (c) the LOS displacement time series of individual PS.The proposed methodology is based on the third product, i.e. the total cumulated displacement over a period covered by a set of SAR images.
In the first stage the relation between ground displacement and the available geothematic layers is studied (Fig. 2, first step).First of all, a spatial analysis is performed between the PSs and geothematic map.This analysis is performed in order to determine if different populations of PSs can be distinguished attending to geothematic layers.However, the Pearson correlation coefficient (i.e. a measure of the strength of a linear association between two variables; McCullagh and  Nelder, 1989) is calculated between the PSs and the geothematic layers.This coefficient allows discriminating if the geothematic layers should be considered as a second input when performing the interpolation.
The second step consists on the interpolation of the cumulative displacement along the satellite LOS; taking into account previously identified PS populations (Fig. 2, second step).For this task, geostatistical tools, coupled with a geographic information system (GIS) (Burrough and McDonnell, 1998), allow us to perform spatial interpolation of scattered measurements and to obtain an assessment of the corresponding accuracy and precision.Kriging (Matheron, 1963) is a geostatistical interpolation technique that considers both the distance and the relation between sampled data points when inferring values at unsampled locations (Journel and Huijberegts, 1978;Isaaks and Srivastava, 1989;Goovaerts, 1997).The kriging smoothing effect on the interpolated maps may be a disadvantage because the reality is expected to be more variable.As a consequence of this effect the variance of the estimated values is lower than the variance of the real val-ues.Geostatistical simulation (Goovaerts, 1997) allows generating multiple equiprobable realizations of the attribute under study, rather than simply estimating the mean.This is a key property of this approach, since a series of realizations representing a plausible range is generated, not just one best estimate.Hence, accuracy can be estimated through distributions of inferred values at unsampled locations using the series of simulated realizations.A large number of geostatitiscally based algorithms exist for the simulation of realizations (i.e, spectral simulation, sequential Gaussian simulation method, Boolean simulation, turning bands, etc.).The conditional SGS method (Gómez-Hernández et al., 1993) has been used in this work because of its long history and wide acceptance for environmental modelling applications.
First of all, each PS population is analysed and, if necessary, transformed into a Gaussian or normalized distribution (Fig. 2, first step).Then, for each normalized PS population, the experimental variograms (Goovaerts, 1997) are established.Variograms are widely used to quantify the spatial variability of spatial phenomena (e.g.Journel and Huijbregts, 1978;Armstrong, 1984;Olea, 1994;Goovaerts, 1997).In this work, the plot of the semivariances as a function of distance from a PS is referred to as variogram: where N is the number of PS pairs separated at a distance h, δ(x i ) are PS cumulated displacement values and δ(x i + h) are all the PS cumulated displacement values at a distance h away from the PS x i .The analysis of this function for different separation distances (h values) allows the characterization of the spatial variability of the PSs.After the calculation of the experimental variogram, a variogram model must be inferred.The parameters of the variogram model are used as input for the SGS.This interpolation method generates multiple equiprobable surfaces of the displacement reproducing the observed data at their locations.To preserve resolution of the satellite images, the chosen pixel size of the interpolated surfaces should be the same as the ground resolution of SAR images.
In the third step from the multiple equiprobable surfaces (Fig. 2, third step), variance and percentile surfaces are retrieved (mean, 68th percentile and 95th percentile).The percentile surfaces are analysed in order to generate the subsidence activity map.The subsidence activity map consists of two maps: the subsidence map and the variance map.The first one, the chosen percentile surface, is a subsidenceinterpolated surface with a confidence level on the interpolation.The second one, the variance, is an estimate of the spatial variability of the subsidence-interpolated surface.Thus, for every pixel, even for those areas in which there were no PS data, the method provides a subsidence value, an associated range of variation and a level of confidence for it.In the fourth step (Fig. 2, fourth step) subsidence activity maps are used to identify buildings that can be damaged by ground subsidence, according to the serviceability limit state (SLS) criterion.The presence of damages mainly depends on the structure's typology and the settlement's magnitude and distribution.The SLS are those conditions that make the structure unsuitable for its projected use.In foundation's design, the most common serviceability limit states are differential settlements (δ s ; i.e. unequal settling of a building's foundation that is computed as the maximum difference between two points from the foundation) and angular distortions (β max ; i.e. the ratio of the differential settlement between two points and the horizontal distance between them), which must be less or equal than the corresponding limiting value stated for them.A vast number of limiting criteria for settlements and angular distortion values are available in the geotechnical literature and standards (e.g.Terzaghi and Peck, 1948;Skempton and McDonalds, 1956;Burland et al., 1977;EN, 1990;CTE, 2006).In this work, considering the cohesive character of the available soils and the high rigidity of the heritage masonry building studied in the city of Orihuela, we have adopted the values shown in Table 1 for angular distortion and differential settlements following the aforementioned authors' recommendations.Note that the limiting values established in the literature correspond to vertical displacements (δ vs ).Therefore, these vertical displacement values were multiplied by the cosine of the satellite look angle (θ = 23 • for ERS and Envisat) (Fig. 3), in order to compute the settlements projected along the LOS (δ vs-LOS ) for a direct comparison between A-DInSAR and settlement-limiting values: (2) Then, according to Eq. ( 2) the vertical differential settlement (δ s ) between two pixels (i and j ) can be expressed as Consequently, considering Eq. ( 3) and an allowable vertical differential settlement (δ s max ), the maximum allowable differential settlement between two pixels (i and j ) along the LOS (δ s-LOS ) is Additionally, the vertical angular distortion (β) between two pixels (i and j ) can be also computed by means of the expression Consequently, considering Eq. ( 5) and adopting a maximum allowable angular distortion caused by vertical settlements (β max ), the maximum angular distortion between two pixels (i and j ) along the LOS (β max-LOS ) can be expressed as The values provided by Eqs. ( 4) and ( 6), considering the limiting vertical values shown in Table 1, are used as limiting values throughout the analysis.Note that the main advantage of projecting the maximum allowable settlement and angular distortion along the LOS is that these values can be directly compared with the displacement values provided by DInSAR.
To calculate the angular distortion and the differential settlement along the LOS for a building, a buffer is performed around each building.The size of the buffer can be obtained by calculating the root square of the diagonal of a pixel.The buffer area is used to extract the cumulative displacement values from the subsidence activity raster maps, and represents the subsidence influence area of the building where damages can be induced.From the selected cumulative displacement values, both the maximum differential settlement and the maximum angular distortion along the LOS are computed for each building (Fig. 3) and compared with the maximum allowable values shown in Table 1.

Data analysis
In this section the above mentioned methodology is applied to the city of Orihuela, where SAR images covering two different periods from July 1995 to December 2005 and from January 2004 to December 2008 have been processed using the SPN algorithm.First of all, the resulting PSI data for both processed periods are described.Then, following the first two stages of the methodology, the results obtained in the city of Orihuela for both periods are presented.

PSI data
PSI displacement maps were retrieved from the SPN processing of 110 C-band SLC SAR images from descending orbits acquired by the European Space Agency (ESA) ERS-1/2 and Envisat ASAR satellites covering two different periods: July 1995-December 2005 and January 2004-December 2008 (Table 2).A similar crop of about 20 km × 8 km was selected from the 100 km × 100 km SAR images, corresponding to the westernmost sector of the Vega Baja of the Segura River.Interferograms were generated from pairs of SAR images with a perpendicular spatial baseline smaller than 800 m, a temporal baseline shorter than 6 and 3 years for the 1995-2005 and 2004-2008  Even though the reference points used for both periods are not strictly the same, they were both located in nearby stable areas of the city of Murcia, 20 km west of the city of Orihuela.Considering the whole data set for both periods, the average displacement rate of the PSs included within the stable lithologies (i.e. the mountain ranges) is below 2 mm year −1 , which is the common stability threshold adopted in the scientific literature for C-band satellite sensors.Moreover, the validation experiments, performed for these SPN data sets with the extensometric network in the city of Murcia (Herrera et al., 2009a, b;Tomás et al., 2011), provided a similar cumulative error (±5 mm) for both periods analysed.Thus, even though the reference point is different for each period, both data sets define the same unstable areas with the same error.
In both periods, 5730 and 4922 PS (Table 3) were detected within the city of Orihuela (280 km 2 area), measuring cumulated displacements between −119 and +67 mm (Fig. 4).Note that negative values represent subsidence and positive values surface uplift.Taking into account the validation experiment performed by Herrera et al. (2009a, b) and Tomás et al. (2011), we define ±5 mm as the stability range for cumulated displacements.Even if the cumulated displacements means are similar for both periods (Table 3), the displacement rate for the period 2004-2008 is almost double that for the period 1995-2005 (Table 3).Additionally, the statistical skewness (i.e. a measure of the degree of asymmetry in

Spatial analysis
Once the preliminary statistical analyses were completed, the relationship between the geothematic layers and the cumulative displacement was studied.Tomás et al. (2010) pointed out that the geology and the thickness of alluvial soft sediments were two important subsidence conditioning factors.
In order to discriminate if the geology is a factor to be considered for the A-DInSAR data interpolation, the PS data were classified in two different ways and then compared one to each other.The first one separates each period of PS in two normal populations (first and second components in Fig. 5).The probability of every PS to belong to the first or second component was computed through the expectationmaximization algorithm for mixtures of univariate normal function from the mixtools R package (Benaglia et al., 2009).This is an iterative process implemented in R software (R Development Core Team, 2010) which allows distinguishing the presence of subpopulations that follow a Gaussian distribution (components) within a global population.The second separates the PS located within the alluvial sediments (Holocene) and those located on the rock outcrops (Tertiary and Permo-Triassic) forming the edge of the basin (Fig. 5).In the period 1995-2005 the PS cannot be separated into two different populations attending to the component analysis (Fig. 5a), as the stable PSs are heterogeneously distributed over the study area.On the contrary, in the period 2004-2008 two different normal populations can be separated (Fig. 5b)  from the component analysis.The first component encloses stable PS on the stable rock outcrops, while the second component includes active PS on the younger and more compressible alluvial sediments.Consequently, through this analysis it is demonstrated that the geology is a subsidence conditioning factor.Therefore, to perform PS data interpolations two different PS populations should be distinguished in each period: one in the stable sediments and another in the more deformable soft sediments.
The second conditioning factor pointed out by Tomás et al. (2010) was the thickness of alluvial soft sediments.Fig. 4 qualitatively shows a spatial correlation between soft soil thickness and cumulative ground displacements, in agreement with Tomás et al. (2010) statements.Consequently, if a high linear correlation exists between the soft soil thickness and the ground displacements, the thickness could be used as a secondary variable to obtain a more accurate ground displacement interpolation.To evaluate if the soft soil thickness could be included as a secondary variable the Pearson correlation coefficient (McCullagh and Nelder, 1989) was calculated between the displacements and the thickness of soft soils.For this purpose the elaborated soft soil thickness map described in Sect. 2 was used.Note that the PSs not included within this map were not considered.The Pearson correlation coefficient is −0.22 and −0.52 for periods 1995-2005 and 2004-2008, respectively (Appendix A).Using Kaiser's (1974) scale, the values obtained indicate that our variables have a poor linear correlation and/or other variables might also be involved.Therefore, the thickness of the soft soils cannot be considered as an additional source of information for improving the displacement data interpolation.The observed lack of high linear correlation can be explained considering that the ground displacements are the result of the combined and superposed effect of piezometric level changes and the soil thickness deformability that do not follow a linear relationship.

Variogram analysis
Before obtaining the interpolation surfaces for each period, the experimental variograms and fitted models should be calculated.The libraries of R software (R Development Core Team, 2010), Gstat (Pebesma, 2004) and geoR (Ribeiro et al., 2001), combined with SGeMS software (Remy et al., 2009) were used for the displacement variogram visualization and analysis.The analysis of the experimental directional variograms, lead to the detection of anisotropy (N60E), which is assumed as the maximum spatial continuity direction.This direction coincides with the basin's axis direction, this fact reinforces the idea that geology is a conditioning factor of subsidence.
The experimental and fitted variogram models obtained for the four PS populations previously defined are shown in Fig. 6 and described in Table 4. Experimental variograms for the 1995-2005 period present higher fluctuations around the sill than the ones computed for the 2004-2008 period.This behaviour can be attributed to the fact that during the period 1995-2005 the groundwater level was under a recuperation phase.Hence, measured subsidence rates are low and sparsely distributed and correspond to the residual consolidation of the 1992-1995 drought period.During the 2004-2008 drought period there was an intensive exploitation of the aquifer in the vicinity of the city of Orihuela, which induced a regional subsidence over the area.This process explains the lower spatial variability of PS displacement values, which corresponds to a lower slope of the variogram models near the origin during this period.Consequently, the analyses of the variogram models confirm a different subsidence spatial variability for each period, which corresponds to different phases (piezometric level recuperation and fall) of the subsidence phenomena.

Results and discussion
In this section, following the last two stages of the methodology presented in Sect.3, the regional and local subsidence activity maps are generated and discussed for both analysed periods in the city of Orihuela.First all, the different percentile maps are analysed.Then, the SLS for 27 historical buildings are calculated and evaluated with the reported damages inventory.Finally a comparison is made between a damage study and the results obtained for the Santas Justa and Rufina Church.

Regional subsidence activity maps
The parameters of the four fitted models (Table 4) are introduced in the SGeMS software to generate 100 equally likely realizations of subsidence by the SGS method.As the ground resolution of ERS1/2 and Envisat images is 20 m, the pixel size of the resulting surfaces has been set to 20 m.Subsidence activity maps are obtained from the SGS results.The mean of these simulations represents the expected subsidence value, which should be very similar to the result of the kriging interpolation (Figs.7a, 8a).The spatial variability of the subsidence is evaluated calculating the variance for the 100 simulations performed (Figs.7b, 8b).The confidence level for subsidence-interpolated values is established calculating the percentiles from the 100 simulations.In our test site, 68th and 95th percentile maps are analysed for both periods (Figs. 7c,d,8c,d).These maps display the threshold subsidence value for which there is a probability of 68 % and 95 %, respectively, that the true subsidence is greater or equal to that threshold.Note that greater subsidence means a more negative value.Visually, during the period 1995-2005, the 68th and 95th percentile maps are the same, being their average cumulated subsidence −17 mm.Note that the mean, which represents a 50 % confidence level, exhibits an average cumulated displacement of −36 mm that is closer to the average of the PS data set (−47 mm).This fact suggests that the mean is the best estimator to elaborate subsidence activity maps in this period.
Concerning the period 2004-2008, 68th and 95th percentiles differ substantially from each other, visually and numerically, with an average cumulated subsidence of −47 and −31 mm, respectively.The latter represents very low subsidence values compared to the −47 mm average cumulated subsidence of the PSs (sampled data), which is closer to the value of the 68th percentile map.Consequently, the 68th percentile is selected to elaborate the subsidence activity map during the period 2004-2008, providing a 68 % confidence level, and the mean for the period 1995-2005 that provides a 50 % confidence level on the interpolation.
Analysing the subsidence activity maps generated with the highest confidence level (Figs.7a, 8c), a lower and a more localized subsidence is shown in the period 1995-2005 than in the period 2004-2008, when a greater and more extended subsidence occurred.This observation is explained by the higher spatial variability (Sect.4.3) estimated during the first period (variance between 0 and 473.07 mm 2 ) than during the second one (variance between 0 and 387.83 mm 2 ).Overall, the geostatistical analysis agrees with the fact that during the period 2004-2008 there was a widespread subsidence due to water withdrawal, whereas for period 1995-2005 a more heterogeneous subsidence was due to the groundwater level recuperation.

Local subsidence activity maps
The purpose at this stage is to identify those buildings susceptible to suffer damage induced by subsidence based on subsidence activity maps.Following the SLS assessment methodology proposed in Sect.3, maximum differential settlements and maximum angular distortions along the LOS are calculated for the 27 historical buildings of the city of Orihuela and compared with allowable values.The buffer area used in this case study is 14 m as the pixel size of the subsidence activity maps is 20 m.Note that for this resolution only five pixels are considered for the SLS calculation of two buildings (22 and 23 in Table 5) due to their small area.In both cases, predictions might not be accurate enough and a higher resolution would be required.The results calculated for the remaining buildings indicate if the adopted SLS criteria exceed or do not exceed the values shown in Table 1; they represent the local subsidence activity map (Fig. 9).In these maps it is observed that the direction of both the maximum angular distortion (Fig. 9c, d) and the maximum differential settlement (Fig. 9a, b) are mainly radial to the Mesozoic reliefs dipping towards the thickest compressible sediments.
According to SLS results (Table 5), the expected damage in the period 1995-2005 affects almost twice the number of buildings than in the period 2004-2008.This is due to the more heterogeneous and variable spatial distribution of the 1995-2005 period displacements, which favours the existence of greater differential settlements and angular distortion.Unfortunately, the lower confidence degree (50 %) obtained in the 1995-2005 period's subsidence activity map could bias the calculation of the maximum differential settlement and the maximum angular distortion along the LOS (δ s max-LOS , β max-LOS ).However, the greater confidence degree (68 %) and the lower variance obtained for the 2004-2008 period's subsidence activity map provides a more reliable calculation of the service-limite state variables.
In the past two decades, damages have been reported at least in 10 historical buildings.For the rest of the buildings there is no information available about suffered damages .A description of reported damages can be found in López (1992), Louis (2005), Louis et al. (2012), Lara (2003), Maciá (2005) and Tomás et al. (2012).This damage description has been complemented with field visits to damaged buildings during the period 2011-2013.The year of the repair and/or stabilization works (Table 5) is a key parameter because after a repair action previous damage cannot be identified in field surveys.Most of the repair actions consist on the replacement of the pavement, the sealing and filling of cracks and the elimination of capillary ascent moisture and leaks (Louis et al., 2012;Lara 2003).However, these actions do not strengthen the foundation of the buildings and further settlements could reopen previous existing cracks or generate new ones.Taking into account that some buildings (e.g.buildings 1 and 6 in Table 5) have been completely or partially underpinned, low displacements and damages are  Considering the heterogeneity of the available information and the previously mentioned constrains that can distort the interpretation of the results, reported damages (Table 5) have only been used for a qualitative evaluation of local subsidence activity maps.Nevertheless, 100 % of the buildings (over 10) where damages have been reported exhibit "medium" to "high" levels of expected damages according to SLS calculations.Note that during the 1995-2005 period 90 % of the buildings (over 10) exceeded the allowable SLS while for the 2004-2008 period this rate dropped to 50 % of the buildings (over 10).In this context, damages observed on several buildings (2, 3, 4, 6, 9, 10, 13 and 20) can be attributed to cumulated subsidence during both periods, since no restoration activities were performed.For these buildings subsidence activity maps show a good agreement with the observed damage since the SLS threshold values were exceeded in at least one of the two periods.
SLS values calculated for building 1 were only exceeded during the period 1995-2005 and coincide with damages reported in 2003 (Table 5).Recent field work permitted to attest that damages in this building did not persist in time since gypsum plaster markers located in the cracks have been intact since the building was underpinned.In the case of building 17, allowable SLS values were exceeded during both periods, even if a good overall condition was recognized during recent field work.This mismatch is explained by the building's restoration in 2009.Additionally López (1992) recognized numerous damages that could be related to the expected "medium" and "high" levels of damage calculated for the 1995-2005 and 2004-2008 periods, respectively.Analyzing those buildings where no information about their SLS state could be gathered, buildings 12 and 15 (Table 5) show in both periods, an expected "high" or "medium" level of damage and should be considered as priority targets in future field works.A second target would include nine buildings (7, 8, 11, 15, 19, 21, 22, 25, and 26 in Table 5) with an expected "high" or "medium" level of damage.The remaining buildings (5,16,18,23,24 and 27 in Table 5) do not show a potential of being damaged and should be inspected in the last stage.Thus with this methodology, two groups of problematic buildings could be identified, establishing priorities for further inspections.

Analysis of the Santas Justa and Rufina Church
A more detailed analysis of the local subsidence activity map is carried out in the Santas Justa and Rufina Church, where a forensic analysis is available.The Santas Justa and Rufina Church is located on the north bank of the Segura River in the centre of the city of Orihuela (Fig. 7 and 8).It was declared a Spanish National Monument in 1971.During the last 20 years damages induced by regional subsidence have affected the structural elements of this church (Fig. 10).Tomás et al. (2012) performed a detailed analysis of the damage suffered by the building, summarized here.In the 1970 s the principal chapel suffered an important tilt (Fig. 10).In the 1990 s, important settlements affected the San José Chapel and the foundation was reinforced in 2002 with micropiles.In the first decade of the present century, displacements affected the whole north zone of the church and the principal chapel area.Multiple cracks where identified in the north, west, and east walls, affecting as well La Comunión Cupola (Fig. 10a), the sacristy and the antesacristy (Fig. 10c).Several plaster markers were placed in the cracks in 2006 and controlled in 2008 (Fig. 10b).Figure 10 shows that multiple cracks grew during this period as a consequence of the sinking of the foundation walls caused by ground subsidence.
Subsidence activity maps, maximum differential settlement and maximum angular distortion along the LOS, of the Santas Justa and Rufina Church are shown in Fig. 11.The average cumulated displacements are −39.2 and −31.9 mm for the 1995-2005 and 2004-2008 periods, respectively, but a greater variance is appreciated for the first period.The maximum angular distortion shows medium (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) to high (2004)(2005)(2006)(2007)(2008) values.Note that during the first period, the maximum angular distortion is located NE of the church, mainly affecting the main chapel.However, it is observed that during the second period the maximum angular distortion is located SW of La Comunión Chapel.The orientation of the two vectors representing the maximum angular distortion is E-W, which is in agreement with the location of reported damages (Fig. 10).Concerning the maximum differential settlement, the vector direction (NE-SW) is the same for both periods.This would indicate that the church would have tilted towards the SW and as a consequence the directions of cracks (normal to tension stresses) are expected to be oriented from NW to SE, coinciding with the reported damages (Fig. 10).Even if the 23.01 mm threshold has not been exceeded in each period, and taking into account that the direction of both vectors is similar, the combination of them yields a 34 mm maximum differential settlement.One can say that the reported damages match, spatially and temporally, with the calculated SLS.

Conclusions
This work proposes a novel method to produce subsidence activity maps based on the geostatistical analysis of persistent scatterer interferometry (PSI) displacement data.These maps permit us to identify widespread subsiding areas and buildings where damages could be produced, providing a useful tool for the planning and management activities of the local authorities.However, the derived results cannot be exploited in an isolated way.They need to be combined with in situ field data to confirm the potential levels of damage.The methodology has been tested in 27 historical buildings of the city of Orihuela (SE Spain), which have been damaged in past decades due to subsidence triggered by groundwater overexploitation.
The spatial analysis and directional variograms revealed that geology is a subsidence conditioning factor that should be taken into account when performing a geostatistical analysis.The variogram-fitted models differ for each period, confirming a different spatio-temporal subsidence behaviour, which is explained by groundwater level changes during droughts.The conditional sequential Gaussian simulation (SGS) provided realistic subsidence estimations including the evaluation of its confidence degree level and its spatial variability.
In this case study, regional subsidence activity maps revealed a lower and heterogeneous subsidence for the period 1995-2005, being more intense and homogeneous for the period 2004-2008.However, local subsidence activity maps permitted identifying buildings susceptible to suffer damages through the calculation of the serviceability limit state (SLS) parameters: the differential settlement and the angular distortion.In 27 historical buildings of the city of Orihuela, SLS results have been compared with reported damages and field checks available for 10 buildings, showing a 100 % success rate.Additionally, SLS allowed distinguishing two groups of potentially damaged buildings, enabling us to establish an inspection priority in agreement with the expected level of damage.
The accuracy of the SLS prediction is limited by the satellite resolution and the area of the targeted buildings.Therefore, even though average-resolution satellite PSI displacement data was used in this work (ERS and Envisat), better results are foreseen using high-resolution satellite sensors (TerraSAR-X, Cosmo-SkyMed, Sentinel).

Figure 1 .
Figure 1.(a) Location and geology of the study area.(b) Piezometric level evolution for the study period.

Figure 2 .
Figure 2. Proposed methodology for the elaboration of subsidence activity maps.

Figure 3 .
Figure 3. Schematic explanation of the serviceability limit state parameters adopted in the analysis.
periods, re-spectively, and a relative Doppler centroid difference below 400 Hz.The digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM) has been used.The PS selection for the estimation of displacements was based on a combination of several quality parameters including lowamplitude standard deviation and high model coherence.

Figure 5 .
Figure 5.PS classification according to normal distributions (first and second components) and geological criteria for both studied periods: (a) 1995-2005 and (a) 2004-2008.

Figure 6 .
Figure 6.Experimental variograms with the fitted theoretical model in 60 and 150 • directions (angles are counted from the north, increasing clockwise).

Figure 9 .
Figure 9. Local subsidence activity maps for both periods studied.Notice that the interpolation surfaces used are the mean and the 68th percentile for 1995-2005 and the 2004-2008 periods, respectively.

Figure 10 .
Figure 10.Location of the structural damage and restoration and reinforcement actions performed on the Santas Justa and Rufina Church (modified from Tomás et al., 2012).

Table 1 .
Adopted SLS criterion for the performed analysis.

Table 2 .
Main characteristics of the processed data stacks of the westernmost sector of Vega Baja of the Segura River: λ, wavelength; θ , look angle; Desc, descending.

Table 3 .
Descriptive statistical measurements for both periods.

Table 4 .
Parameters of the variogram models for each normalized PS data group.

Table 5 .
Note that angular distortions (β max-LOS ) and (δ s max-LOS ) values provided by InSAR for both periods have been classified using the font format set in Table1(bold for negligible, underlined for medium and italic for high).NA: not available; R: Renaissance; G: Gothic; B: Baroque; N: neoclassical; U: unknown.( * ) Total or partial underpinned foundation.