Application of Insar and Gravimetry for Land Subsidence Hazard Zoning in Aguascalientes, Mexico

In this work we present an application of InSAR and gravimetric surveys for risk management related to land subsidence and surface ground faulting generation. A subsidence velocity map derived from the 2007–2011 ALOS SAR imagery and a sediment thicknesses map obtained from the inversion of gravimetric data were integrated with a surface fault map to produce a subsidence hazard zoning in the city of Aguascalientes, Mexico. The resulting zoning is presented together with specific recommendations about geotechnical studies needed for further evaluation of surface faulting in these hazard zones. The derived zoning map consists in four zones including null hazard (stable terrain without subsidence), low hazard (areas prone to subsidence), medium hazard (zones with subsidence) and high hazard (zones with surface faulting). InSAR results displayed subsidence LOS velocities up to 10 cm/year and two subsidence areas unknown before this study. Gravimetric results revealed that the thicker sediment sequence is located toward north of Aguascalientes City reaching up to 600 m in thickness, which correspond to a high subsidence LOS velocity zone (up to 6 cm/year).


Introduction
Land subsidence induced by groundwater extraction is a man-induced geological hazard affecting many cities in the word.One of the main hazards on ground subsiding areas is the development of subsidence-related surface faults and earth fissures, because they damage housing and other infrastructure, decreasing their real estate value.In Aguascalientes (see Section 3.1), surface faults can develop displacement across their escarpment (Figure 1), with a variable width of an active zone up to ten meters in which housing structures are easily damaged.Surface faults usually develop over the bounds of the subsidence zones, but they are not rare toward the central part of subsiding areas.Another subsidence-related problem is the increased flood likelihood due to the disruption of sewage utilities and changes in the surface drainage.In this work we use the terms "ground fault", "surface fault", "surface crack", or "fissure" to refer to subsidence-related terrain discontinuities, and the terms "quaternary fault", "pre-existing fault" or simply "fault" refer to pre-subsidence tectonic faults.
Land subsidence due to groundwater extraction is a slow and gradual process whose effects are usually observed long after subsidence has begun.Because there is no loss of life usually associated with its occurrence, subsidence may not be considered a major disaster.Hence, land subsidence is usually not considered in risk mitigation public policies in Mexico as opposed to other hazardous events such as earthquakes, volcanic eruptions, hurricanes and floods.Unlike subsidence, these events are characterized by short duration and high intensity with consequences and catastrophic effects that are immediately observable.
The Disaster Risk Index (DRI) defined by the United Nations Development Programme [1] provides a methodological tool to assess the impact of catastrophic events [2].Currently, the DRI is solved for three natural hazards: earthquake, tropical cyclones and flooding, but risk assessment of other hazards including land subsidence due to groundwater extraction remains as scientific challenges and as works in progress.
In general, risk assessment involves three basic elements: (a) characterization of hazardous event, which is the natural or anthropogenic event able to cause life-loss or property damage (b) quantification of physical exposure, which refers the number of lives or assets exposed to the hazardous event; (c) determination of vulnerability, which is the features that make the physical exposure able to absorb the impact of the hazardous event [1,2,30,31].Once the elements of risk have been determined, risk is calculated through a model which relates these elements of risk.
Risk subsidence assessment has been addressed in different ways.For example some works propose methodologies for the determination of areas prone to develop subsidence-related faulting [23,26,27,[32][33][34].Other works proposed methodologies to determine the areas that are likely to develop subsidence [29,35].In all these approaches, the physical exposure has been the terrain but not constructions.
In this work we propose an approach for producing subsidence hazard maps, where man-made structures are the assets exposed to the hazardous event, using the city of Aguascalientes, Mexico as a study case.We processed InSAR data to derive subsidence velocity maps which were combined with geophysical and geological information to produce a zoning map of subsidence hazard.
The presented work covers only the first part of the subsidence risk assessment process: the characterization of the hazardous event.Nevertheless, the resulting maps from this approach provide valuable information for damage prevention in cities subject to subsidence.

Identification of Factors Contributing to Subsidence Hazard
Hazard is usually characterized in terms of the probability that a hazardous event occurs with certain intensity during a given time period [30].This approach works for phenomena such as earthquakes or hurricanes, which cannot be avoided and whose magnitude cannot be reduced.In those cases, historical data of the events and their magnitudes are processed with probabilistic techniques in order to estimate the magnitude of an event in a particular time period.
However, land subsidence due to groundwater extraction is a predictable phenomenon, whose occurrence can be avoided if groundwater is not extracted, and whose intensity can be reduced by changing the ground water extraction policies.The mechanism leading to the generation of subsidence due to groundwater extraction and associated faulting is well-known [7,[36][37][38][39][40].Nevertheless, the ability to calculate time, location, speed and magnitude of subsidence is limited due to the difficulty to determine the parameters that control the process [41,42].
Subsidence occurrence is a combination of two factors: (a) the existence of unconsolidated or poorly consolidated sediments deposits that comprise the aquifer system and (b) lowering of the groundwater level [37,38,43].The existence of large thickness of sediments prone to consolidation in the subsoil, which contain water susceptible to be pumped, is by itself the geological environment potentially prone to subside, whilst water table lowering is the triggering factor of subsidence.However, if the aquifer system is not comprised of unconsolidated sediments, then subsidence will not develop, even if a groundwater level reduction takes place.Thus, the existence and spatial distribution of consolidation-prone basin-fill sediments is the first independent factor to be considered for analyzing land subsidence hazard zoning, which is the mapping of the areas prone to subsidence and the lowest level of hazardous for constructions.
In areas where subsidence began long ago, subsidence can evolve in two manners: (a) subsidence develops in a uniform way generating a vertical displacement field with low horizontal gradients; (b) subsidence develops differentially, then the horizontal gradients of the vertical displacement field are large enough to develop horizontal stresses leading to ground failure development including surface faults and cracks [7,23,44,45].
Uniform subsidence changes the surface slope and modifies the slope of the sewer utilities, developing new flooding-prone areas.Hence, zones where uniform subsidence occurs are the second hazard level for civil structures and urban infrastructure.
Surface faulting is the main concern for property owners and local governments in subsiding urban areas not only for their damage potential but also because ground failures can severely reduce real estate value.As a result, areas where surface faults have developed pose the highest hazard level to civil structures and other urban infrastructure.
In summary, the key elements to assess the hazard level for man-made structures due to land subsidence by groundwater extraction are: (1) determination of existence of deformable sediments containing groundwater; (2) detection of sediment zones with subsidence in progress; and (3) location of affected areas by surface faulting.

InSAR and Gravimetry for Determining Areas Having Subsidence Hazard Factors
Both gravimetry and InSAR are useful tools for determining two of the three subsidence hazard factors previously discussed.Gravimetric analysis can be used to determine the distribution of sediments prone to be consolidated, while InSAR detects the affected areas by subsidence.The presence of surface faults, which is the other subsidence hazard contributing factor, can be obtained through field based cartography.
Gravimetric measurements have been used successfully to determine the distribution and thickness of granular filling in sedimentary basins [23,34,46].In central Mexico, areas undergoing land subsidence due to groundwater extraction are usually comprised of sequences of unconsolidated sediment overlying a more dense volcanic or sedimentary rock formation.Consequently, the density contrast between the underlying consolidated rock and the overlying sediments is significant, which is very favorable for a gravimetric study.
The result of a gravimetric survey is a gravimetric anomaly map, which is the difference between the measured and theoretical gravity field.This difference is attributed to the density heterogeneities in subsurface: low gravimetric anomalies suggest presence of low density material close to the surface, and high values of the anomaly indicate denser material strata close to the surface.Hence, a gravimetric analysis provides valuable information about the spatial distribution and thickness of sediments.
Furthermore, some studies have shown that in geological settings susceptible to developing subsidence by groundwater extraction, there is an inverse correlation between sediment thickness and the gravimetric anomaly [23,34,47].This relationship can be used as a qualitative way to characterize areas with large sediment thickness.
Additionally, gravimetric anomaly data can be inverted or directly modeled in order to elaborate models of the sediment thickness distribution.The modeling may be enhanced by constraining the model with lithological and other direct observations, allowing the generation of detailed maps of sediment thickness distribution Because land subsidence due to groundwater extraction may affect large areas, its detection and quantification are quite suitable for satellite remote sensing techniques.interferometric synthetic aperture radar (InSAR) techniques have been successfully used to characterize subsiding areas [11,48,49].

Study Area
The city of Aguascalientes is located within the Aguascalientes graben in central Mexico, 430 km NW of Mexico City.Close to one million inhabitants live in the city and suburban municipalities, 725,000 in the city of Aguascalientes, and the other 225,000 in the 7 surrounding municipalities [50]: Cosio, Jesús María, Rincón de Romos, Pabellón de Aretaga, San Francisco de los Romo, San Pedro Piedra Gorda and Luis Moya (Figure 2a).Intense groundwater extraction initiated in the early 1970's due to an increase in agricultural and industrial activities, triggering land subsidence and development of surface faults [51][52][53] and even the reactivation of tectonic faults [7].
According to SIFAGG (Sistema de Información de Fallas y Grietas) [54], currently 208 surface faults and fractures have been mapped throughout the entire Aguascalientes valley with an accumulated length of 290 km, affecting 1865 buildings mainly housings, from which 1438 of those are located within the city of Aguascalientes.Figure 2b shows only those surface faults within Aguascalientes City.

Determination of Deformable Sediment Distribution
The total thickness of the unconsolidated sediments in the study area (Figure 2b) was determined through a gravimetric study, which consisted in the surveying and processing of 339 ground-based gravimetric measurements, according to Telford et al. [55] and using a Scintrex CG5 gravimeter.
Figure 3a shows the Bouguer's anomaly of the study area.Gravity data were modeled according to Singh and Guptasarma [56] using PyGMI software developed by Cole [57] in order to determine the shape of the sedimentary package and underlying rocks.The software works out the gravimetric anomaly produced by an initial model of the rocky stratum and the sediments distribution.That is the calculated anomaly.The initial model is modified until the calculated anomaly matches the anomaly determined by processing field data (measured anomaly).The model is more realistic if there are data for constraining some points of depth to rocky stratum, or density information of the subsoil materials.We used 1500 ˆ1500 m and 500 m high blocks for the gravimetric modeling.Figure 3b shows the calculated Bouguer's anomaly from the obtained model of bedrock and sediment thicknesses distribution.
The geology of the Aguascalientes Valley was broadly described by Aranda-Gómez [51] and Loza-Aguirre et al. [58].For the purpose of this work, we used a simplified stratigraphic column of the subsiding area elaborated with lithological logs wells information (Figure 4).The densities associated to each stratigraphic unit used for the gravimetric modeling were those reported by Pacheco-Martínez et al. [59].The non consolidation-prone units are composed of a polycmitic conglomerated, rhyolite and ignimbrite, and the consolidation-prone unit is a quaternary alluvial sequence with a diverse content of silt, sand and gravel.The top surface of the rocky stratum (Figure 5a) was constrained in several locations through information of logs wells lithology and rock outcrops on both sides of the graben.Finally, an isopach map (Figure 5b), which is a sediment thickness variations map, was calculated from the resulting elevation difference between the top bedrock surface of the conglomerate and volcanic rocks and the ground surface.In order to characterize the subsidence field within the Aguascalientes graben, a subsidence LOS (line of sight) velocity map was obtained from 34 ALOS scenes (ascending track 191, frames 420, 430) acquired between August 2007 and March 2011.We used the ROI_PAC software developed by NASA's Jet Propulsion Laboratory [60].Imagery was processed with the differential InSAR technique (D-InSAR) to obtain 28 interferograms to derive an InSAR subsidence LOS velocity map.The InSAR time series analysis was obtained using the Short Baseline Subset Analysis (SBAS) technique [61][62][63].Topographic correction was applied according to Fattahi and Amelung [64].A threshold of 0.7 temporal coherence was used for the resulting velocity map (Figure 6a).The high calculated temporal coherence of the study area indicates a minimal effect from unwrapping errors.

Determination of the Subsidence-Affected Area
Figure 6a shows three subsiding areas.The largest one corresponds to Aguascalientes Valley, where the rate of subsidence has reached 10 cm/year in two locations (Rincón de Romos and José Gómez Portugal).Both locations are agricultural areas with intense groundwater pumping.The other two subsidence areas located at north Loreto and east NISSAN II (Figure 6a), also correspond to agricultural zones.These two subsiding areas located outside of Aguascalientes Valley had not been documented before this study.Subsidence velocity in both areas is up 6 cm/year.
Figures 2a and 6a show that subsidence is developing in flat areas surrounded by topographic elevations of rocky outcrops.The velocity map of Figure 6a shows that during the period from 2007 to 2012 the area of rocky outcrops encircling at the Valley did not develop subsidence.This is consistent with other works which indicated that subsidence occurs only within the Valley where a significant thickness of sediments exists.Hence, the rocky outcrops located outside the valley (dark blue areas in Figure 6a and surrounding mountains in Figure 2a) can be considered stable zones and reference points for monitoring areas with subsidence.Subsidence in Aguascalientes has a nonlinear trend as shown by a subsidence time series from the INEGI GPS station (Figure 6c).This graph shows that subsidence velocity is decreasing more and more since the beginning of the records, but it has been at a steady pace since 2005.

Surface Faults Mapping
An updated subsidence-related surface fault cartography was generated using both information available in previous surface fault maps and field verification of recently generated surface faults not included in the previous maps.The resulting surface fault location was integrated into Figures 2b and 6b.
Figure 6b shows that subsidence-related surface faults preferentially develop within the subsiding area and along its limit where differential subsidence favors their formation.Nevertheless, some segments of them (1, 2 and 3 in Figure 6b) are located away from the limit of subsidence area, in locations where the subsidence velocity map shows null subsidence.The location of these segments of surface faults coincides with pre-subsidence tectonic faults [7] which are cutting the sediment unit.These segments of surface faults located outside of the subsiding area suggest that preexisting tectonic faults are reactivated by the influence of the subsidence stress field.

Subsidence Hazard Zoning
The subsidence hazard zoning map shown in Figure 7c was calculated from the hazard matrix described in Table 1.The hazard level is directly dependant on the existence of those factors discussed in Section 2.1.For example, the null hazard level (last column in Table 1 and green tones in Figure 7c) corresponds to the area where there is not a single hazard factor: consolidation-prone sediments do not exist (Figure 7a), and subsidence and faults have not been observed (Figure 7b).Conversely, the high hazard zone (second column in Table 1 and red lines in the hazard map of Figure 7c) is the area affected by surface faults, and in which the three hazard factors are present, except for the surface faults segments 1, 2 and 3 discussed in Section 3.3.In each of the resulting hazard zones, the effect of ground subsidence on the existing civil structures and urban infrastructure is different, and the prevention and mitigation actions must be in accordance with the type of potential threat.We provide recommendations for geotechnical studies that are needed to assess subsoil conditions according to the hazard level, both for any future real estate development and for those existing structures.We list the recommendations from lowest to highest hazard level.
Null hazard zone: This corresponds to the area in which land subsidence was not detected during the period of InSAR data acquisitions (subsidence LOS velocity = 0 in Figure 6), and surface faults have not been observed yet.Furthermore, in this zone the total thickness of unconsolidated sediments is negligible.Null hazard zone corresponds to the rocky outcrops area and those areas in which the bedrock is covered by a thin stratum of sediments.Standard geotechnical studies are sufficient for this zone.
Low hazard zone: This is the zone of unconsolidated sediments forming the overexploited granular aquifer system (sediment thickness > zero in Figure 5), which are bounded by rocky outcrops, and in which subsidence LOS velocity was zero for the observation period (dark blue tones in Figure 6) and surface faults have not been observed.
In this zone, standard geotechnical studies are needed, plus a superficial geological study searching for evidence of tectonic faults and cracks in the sediment unit, which could be reactivated if a decline in groundwater table take place in this zone.
Medium hazard zone: This corresponds to the zone of unconsolidated sediments forming the overexploited granular aquifer system (sediment thickness > zero in Figure 5), and in which subsidence was detected through the LOS velocity map (subsidence LOS velocity > 0 in Figure 6).Although this zone is undergoing subsidence, it is presented uniformly, such that terrain ruptures have not developed.
Geotechnical studies similar to those suggested for low hazard zone are needed, plus a geophysical survey for detecting any incipient or blind surface faulting.Resistivity profiles have been used with relative success for detecting surface faults in their initial stage when they may not be visible on the surface yet [7].In case that resistivity profiles show an anomaly that could be related to an incipient surface fault, direct excavation of a trench may be needed in order to confirm its existence.Another threat to buildings in this area is the modification of surface drainage as a result of subsidence, which could develop flooding areas and also renders sewage systems ineffective.As this is a regional effect of subsidence, the municipal agency responsible for operating the sewage system along with the agency responsible for urban planning should take this issue into consideration for developing timely mitigation measures.
High hazard zone: This zone corresponds to areas where subsidence-related surface faults and fissures have been detected due to the damages they cause to constructions and terrain surface along their trace.The zone includes a band of terrain in each side of the trace of the surface fault in which subsidence is developing in a differential manner.
Geotechnical studies similar to those suggested for the medium hazard zone are needed, plus a detailed analysis to determine the influence width of subsidence-related faults in which civil structures may be damaged due to uneven subsidence.

Discussion
The approach of this work was to elaborate a zoning map of the hazard to which constructions are exposed in subsidence areas.The resulting hazard zoning map includes areas prone to subsidence, those currently undergoing ground subsidence and existent surface fault traces.The hazard map does not include the zones prone to surface faults generation or future zones of surface faults growth.
Although the subsidence-induced faulting mechanism due to groundwater extraction is very well understood, examples of determination of parameters involved in surface fault generation are scarce, and their field and laboratory measurements are still not a common practice for many real estate developers.As a consequence, subsidence-related fault modeling customized for specific civil structures that may be derived in individualized construction design is still an unattainable goal.Some authors have proposed solutions for determining the zones prone to develop subsidence-related ground faulting [15,23,27,32,33].These methodologies are based on the determination of the horizontal gradients of a measured parameter.The proposed solutions correlate ground faulting areas to high values of horizontal gradients of gravimetric anomalies [23], soil vibration frequency [32,33] and high values of horizontal gradients of subsidence magnitude [15,27].Their results show that these methodologies can forecast shallow ground fault generation, but further research is still needed to determine critical gradient threshold values related to ground faulting developing in specific cases.Also, further work is still needed to improve the temporal and spatial accuracy for surface fault generation predictions.
The most hazardous zone for constructions is on ground failures.The zoning includes only the linear feature representing the trace of the terrain rupture.However, ground failures have an active zone defined by the width of the trace in which differential subsidence is significant to induce damage to the constructions.Studies to determine the width of the active zone or the influence zone of a ground failure are practically nonexistent.More research is needed to obtain a reliable and practical methodology to determine accurately enough this parameter.
The recommendations in this work for geotechnical studies in subsidence areas are to explore subsoil in order to prevent effects of subsidence, mainly those related to surface faults.The current practice in Aguascalientes City and other Mexican cities where surface faults have damaged constructions is the continuous repair of buildings as long as ground failure does not affect their structural stability.Otherwise, buildings are demolished and the terrain is used for parking, green area or other uses except for building construction.Performance investigations of constructions built over ground failures are incipient [65,66].Preliminary results show that a combination of flexible materials and special structures, including a device for restoring the level of construction at certain intervals, could be a solution to prevent damage to constructions.
As subsidence is a progressing deformation and rupture subsoil process, hazard levels could change with time.If subsidence continues, then new surface faults and fissures will be generated, and the existing ones will enlarge, producing new zones of high hazard.Conversely, if subsidence stops, then terrain surface will not experience differential subsidence.Surface faults and fissures will become inactive, and the terrain surface will be stable for constructions.In any case, the subsidence hazard zoning map should be periodically updated to prevent the implementation of wrong measures in risk management.
Parameters to determine hazard factors were obtained directly by on-site measurements (gravimetry, mapping of subsidence-related surface faults and rocky outcrops) and by remote sensing (InSAR).Hence, accuracy of the zoning depends on the resolution of the methods used for processing of measured field data.
In the case of the determination of rocky stratum and sediment distribution (low hazard zone), we used 1500 ˆ1500 m and 50 m high blocks for the gravimetric modeling.Hence, the resolution achieved for the map production was 1500 m.Therefore, features with sizes lesser than 1500 m are not represented in the model of rocky stratum (Figure 5a).Hence, the model of rocky stratum might not include detail sufficient to identify topographic features that are triggering differential subsidence and causing increasing subsidence in specific locations.
However, the resolution of sediment thickness map distribution does not affect the accuracy of hazard zoning, because the limit of zones of low and null hazard was defined from a GPS surveying on the bounds of rocky outcrops with an accuracy of ˘4 m.Likewise, the limit of zones of low and medium hazard is not influenced by the resolution of sediment thickness map because this limit was determined by subsidence measured by InSAR.
The ALOS-InSAR images used to produce the subsidence velocity map have a resolution of 30 m.Hence, as we used a GPS with sub-meter resolution for surveying of rocky outcrops and surface faults, the greatest error in zoning map limits could come from the subsidence velocity map.However, this error is negligible for scale maps 1:10,000 which is the scale more frequently used for urban planning management.
The zoning map does not include either the effect of lateral variations in the sediment thickness or the influence of subsidence magnitude developed in specific locations.The insertion of these factors would improve the zoning process, which would result in a map with more zones than the map derived from the factors analyzed in this work.
Each of such zones would indicate the hazard related to subsidence more accurately.Nevertheless, the simplification of the hazard zoning in four zones as was proposed in this work allows specific and practical recommendations for geotechnical exploration.
A map with many zones or a continuum zoning map could be an interesting scientific contribution.However, the resulting zoning map could lose its practicality and usefulness for urban planning and development.Additionally, zones with higher subsidence can be addressed as special and specific cases of medium hazard level (as determined in this work).
Risk assessment involves three stages: (1) hazard characterization; (2) quantification of exposed elements at the hazard and determination of their vulnerability; and (3) calculation of risk.Although the characterization of hazard factors is really relevant only in the context of risk assessment, our work represent a significant contribution to the development of a further complete methodology, mainly because stages 2 and 3 are still unresolved issues and lines of research in progress.

Conclusions
A hazard subsidence zoning map is a necessary element in risk management of subsidence effects on man-made structures.Gravimetric measurements and InSAR-derived subsidence velocity maps provide valuable information for this hazard map.Gravimetric surveys are quite suitable for determining sediment thickness which may be prone to consolidation and consequently develop ground subsidence, while InSAR techniques are most suitable for a precise characterization of the subsidence field.
A combination of both techniques, along with reliable subsurface surface faulting information, was used to derive a subsidence hazard map for the city of Aguascalientes.This map will allow state and municipal government agencies to clearly specify specific geotechnical and geophysical studies according to the hazard level of the zone in which new constructions are planned.
Due to the dynamic nature of the subsidence process, hazard zoning maps, such as the one described in this work, need to be continuously updated.Future InSAR analysis using newly acquired data or enhanced processing techniques along with other geodetic infrastructure development such as continuously operating GNSS stations will provide updated ground subsidence velocity maps which, in turn, will allow periodic updates of the subsidence hazard maps.

Figure 1 .
Figure 1.Surface faults due to differential subsidence in Aguascalientes.(a) Surface fault affecting man-made structures; (b) Surface faulting has developed escarpments up to 1.80 m high.The location of the subsidence-related surface faults is shown in Figure 2b.

Figure 2 .
Figure 2. Location map of the study area: (a) Aguascalientes Valley and urban areas within; (b) Study area of Aguascalientes City.

Figure 3 .
Figure 3. Gravimetric anomaly of the study area.(a) Measured anomaly; (b) calculated anomaly from model.

Figure 4 .
Figure 4. Simplified stratigraphic column of the study area.

Figure 5 .
Figure 5. (a) 3D model of the surface topography and bedrock topography; and (b) isopach map of sediment thickness in the city of Aguascalientes.

Figure 6 .
Figure 6.Subsidence velocity maps for (a) Aguascalientes Valley; (b) Aguascalientes City including surface faults and (c) Subsidence graph from INEGI GPS station.White circled numbers indicate surface faults outside of subsidence area.

Figure 7 .
Figure 7. (a) Total sediment thickness map; (b) subsidence LOS velocity map; (c) subsidence hazard zoning map of Aguascalientes City.Hazard zones are shown over a shaded relief for reference.

Table 1 .
Hazard factor matrix of subsidence effects on constructions.