Long Term Subsidence Analysis and Soil Fracturing Zonation Based on InSAR Time Series Modelling in Northern Zona Metropolitana del Valle de Mexico

Abstract

In this study deformation processes in northern Zona Metropolitana del Valle de Mexico (ZMVM) are evaluated by means of advanced multi-temporal interferometry. ERS and ENVISAT time series, covering approximately an 11-year period (between 1999 and 2010), were produced showing mainly linear subsidence behaviour for almost the entire area under study, but increasing rates that reach up to 285 mm/yr. Important non-linear deformation was identified in certain areas, presumably suggesting interaction between subsidence and other processes. Thus, a methodology for identification of probable fracturing zones based on discrimination and modelling of the non-linear (quadratic function) component is presented. This component was mapped and temporal subsidence evolution profiles were constructed across areas where notable acceleration (maximum of 8 mm/yr²) or deceleration (maximum of −9 mm/yr²) is found. This methodology enables location of potential soil fractures that could impact relevant infrastructure such as the Tunel Emisor Oriente (TEO) (along the structure rates exceed 200 mm/yr). Additionally, subsidence behaviour during wet and dry seasons is tackled in partially urbanized areas. This paper provides useful information for geological risk assessment in the area.

1. Introduction

The Valley of Mexico is a highlands plateau in Central Mexico which is surrounded by mountains of volcanic origin and characterized by the common presence of faults and high seismic activity [1]. Subsidence in the region occurs as a result of intensive groundwater extraction and constitutes one of the most harmful and costly hazards impacting the Valley of Mexico. Furthermore, in addition to the regional systems of faults identified by various scientists (e.g., [2,3,4]), fracturing risk zones might also be induced by differential sinking rates. Subsidence and its associated phenomena consequently cause damage to infrastructure, worsen the water quality and increase flooding risks (e.g., [1,5,6]).

The Valley of Mexico subsidence is a complex process that varies spatially and temporally. Thus, techniques that enable immediate data generation for supporting decision making at low cost are required. Several multidisciplinary studies have been focused on the south of the Zona Metropolitana del Valle de Mexico (ZMVM) to measure, analyse and understand the subsidence problem (e.g., [7,8,9,10,11,12,13]). However, a few studies considered the entire northern ZMVM (e.g., [14,15,16]), characterized by significant urban sprawl over areas of geological risk (either natural and/or human-induced) and potential natural recharge zones for the aquifer.

Radar remote sensing offers excellent opportunities (large spatial coverage at relatively low cost) to rapidly derive information useful to assess risks associated with several hazards (e.g., volcanos, landslides, earthquake, subsidence). The standard Synthetic Aperture Radar Interferometry (InSAR), for example, has been proved to be a valuable tool for deformation studies (e.g., [17,18,19]); however, limitations due to geometric and temporal decorrelation and degradation due to effects from atmospheric artefacts exist. This technique particularly fails for motion detection in vegetated and cropping land areas, as found in the northern part of the ZMVM where ~30% of the territory corresponds to this type of soil coverage. Moreover, the sparse temporal coverage from conventional InSAR does not provide enough information about the temporal evolution of the deformation, i.e., the number of highly coherent interferograms spanning a short period of time is limited and combination of images that produces optimal results for deformation studies is not always possible [20]. Thus, the application of multi-temporal InSAR (MTI) methods such as Permanent Scatterers (PS) (e.g., [21,22]), Small Baseline Subset (SBAS) (e.g., [23,24]) and their combination (e.g., [25]) provides the possibility to partially overcome the restrictions of the conventional InSAR approach, enabling the detection and measurement of deformation with sub-centimetre precision.

A previous work in the northern ZMVM [26] presented a subsidence analysis limited to the 2002–2010 period using ENVISAT data only. Here, the subsidence and associated fracturing investigation in northern ZMVM by means of multi-pass InSAR is extended to approximately cover the 1999–2010 period, using both ERS and ENVISAT acquisitions. In addition, an approach is herein proposed for the automatic detection of areas prone to fracturing in the north part of the ZMVM based on the modelling of the non-linear component. This methodology includes the construction of subsidence evolution profiles where notable accelerated or decelerated deformation is observed. The analysis focuses on two profiles along and approximately perpendicular to the hydraulic drainage infrastructure of the Tunel Emisor Oriente (TEO), with a 7.5 m diameter, 62 km long and with 24 vertical shafts distributed along its length (see Figure 1). Correlation between the identified fracturing areas, groundwater extraction and location of different geotechnical units is also addressed. Influence of the dry and wet seasons in the subsidence behaviour is likewise reviewed.

Figure 1. (a) Location of the Mexico Basin; (b) SBAS area of analysis (red square) within the Mexico Basin (white area). The ZMVM is represented in grey and the capital Mexico City indicated by the red star; (c) Main features of the study area are represented: Mountainous range, hills and ancient lakes (Texcoco, Xaltocan and Zumpango). Railways and the TEO vertical shafts are depicted, as well as, the Registro Publico de Derechos de Agua (REDPA) and the Pozos de Accion Inmediata (PAI) well branches.

Figure 1. (a) Location of the Mexico Basin; (b) SBAS area of analysis (red square) within the Mexico Basin (white area). The ZMVM is represented in grey and the capital Mexico City indicated by the red star; (c) Main features of the study area are represented: Mountainous range, hills and ancient lakes (Texcoco, Xaltocan and Zumpango). Railways and the TEO vertical shafts are depicted, as well as, the Registro Publico de Derechos de Agua (REDPA) and the Pozos de Accion Inmediata (PAI) well branches.

2. Study Area

The Valley of Mexico belongs to the Trans Mexican Volcanic Belt and it encompasses the ZMVM, which is the most important cultural, economic and industrial centre in Mexico. The ZMVM, with a surface of approximately 4715 km² and which is located in the western part of the Mexico Basin at an average altitude of 2240 m above sea level (m.a.s.l.), is surrounded by several mountainous ranges and volcanoes that reach up to 5500 m.a.s.l. The area of study in the present paper is limited to the northern ZMVM, between the meridians 98°57′W, 99°18′W and the parallels 19°49′N, 19°29′N (Figure 1). This area, similar to other parts of the ZMVM, has been experiencing an accelerated demographic growth (~2.13% for the period 2000–2010 with a population of ~5,742,000 inhabitants for the area of study by 2010; [27]) and sizable volumes of water from the Mexico Basin aquifer system (~1830 hm³/year; [28]) have been extracted to satisfy the local water needs, leading to its overexploitation. Consequently, a loss of the pore pressure is induced producing the compaction of the clay rich deposits in the lacustrine area where the city and the study area are found. Particularly in the north of the Valley of Mexico, the water balance of the aquifer (difference between extraction and recharges) results in a deficit estimated by 2011 in −194.97 hm³/year and it is expected to increase for 2021 to −236.29 hm³/year [15], thus possibly increasing subsidence in the region.

Furthermore, multidisciplinary studies in the area are essential since ground sinking and associated fracturing can negatively impact important hydraulic projects such as the TEO. This tunnel is aimed to conduct the excess of water out of the Mexico Basin, thus reducing the risk of flooding and helping maintenance of existing drainage structures.

3. Advanced InSAR Analysis

3.1. SAR Dataset and MTI Processing

Fifty one ENVISAT images and 11 ERS 1/2 acquisitions were used in this study, both acquired along descending orbits. The ERS images span the time interval between February 1999 and December 2000 and the ENVISAT data set covers from 2002 to 2010. An area of about 30 × 30 km was cropped from all images, corresponding to northern ZMVM.

The set of available raw SAR images were focused to get the Single Look Complex (SLC) product by using the Repeat Orbit Interferometry PACkage (ROIPAC) processing software [29], and the interferogram formation was carried out by the Delft object-oriented radar interferometric software (DORIS) [30]. Orbital effects in the ERS and ENVISAT interferometric pairs were diminished by using precise orbits [31] provided by the Delft Institute for Earth-Oriented Space Research (DEOS) and the European Space Agency (ESA), respectively. Subtraction of topographic fringes from each interferogram was performed using the 3-arc second Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM), version 4 [32].

The MTI analysis was carried out using the SBAS module of StaMPS [22,25], which was adapted for selection of an appropriated spatial reference selected on the basis of geological information. Most methods propose an initial linear model of the deformation [21,33,34,35], limiting their capacity to identify complex earth processes that do not follow the presupposed behaviour. On the other hand, the StaMPS algorithm makes no assumptions about the temporal evolution of the deformation, characteristic that makes it capable of detecting pixels in natural terrains or undergoing non-steady deformation. This property is critical in the present analysis since subsidence and associated soil fracturing in the study area impacts both urban structures and agriculture crops.

The SBAS analysis was applied to the set of available ERS and ENVISAT data. All SLC ERS and ENVISAT radar scenes were first coregistered to the 3 December 1999 and 12 December 2005 acquisitions, respectively. Note that in both cases the “master” image belongs to the dry season (November–May) to reduce possible atmospheric contributions. In order to reduce the geometrical and temporal decorrelation, interferograms with perpendicular baselines lower than 600 m and temporal baselines lower than 350 days were constructed (Figure 2). A total number of 18 ERS and 98 ENVISAT interferograms that provide links between all the corresponding acquisitions were built. The flat-earth and topographic effects were computed and subtracted from the interferometric phase, and spectral filtering was applied to the remaining signal to further reduce the geometric decorrelation.

Figure 2. Temporal (X-axis) and perpendicular (Y-Axis) baselines. Black dots represent the SAR acquisitions and the grey lines the interferometric links. (a) ERS data; (b) ENVISAT data.

Figure 2. Temporal (X-axis) and perpendicular (Y-Axis) baselines. Black dots represent the SAR acquisitions and the grey lines the interferometric links. (a) ERS data; (b) ENVISAT data.

The study area was divided into several patches and slowly-decorrelating filtered phase (SDFP; [25]) pixels with amplitude difference dispersion D_ΔA [25] larger than 0.62 were initially weeded. Afterwards, a phase stability analysis was performed and SDFP pixels with a standard noise deviation lower than 1.2 and 1.7 for the ENVISAT and ERS data sets, respectively, were selected for further processing. Larger decorrelation was observed in the ERS pairs, thus a larger standard noise deviation was chosen to increase the number of selected pixels from this data set, particularly over some areas of interest (e.g., cropping land areas).

The unambiguous phase was solved by a 3D unwrapping algorithm [36] which takes advantage of the temporal information. The method applies the SNAPHU algorithm [37] for unwrapping the phase in space, being first the irregular network of the SDFP pixels reduced to a regular grid [38]. Finally, the cumulated phase delays are retrieved by using Weighted Least Squares (WLS) with weights estimated from the phase information of each small baseline (SB) interferogram. The SBAS analysis is complemented by maps of system misclosure built from the computation, for each pixel, of the Root Mean Square (RMS) between the observed interferometric phase and the one reconstructed from the estimated phase delays (as in Cavalié et al. [39]). These maps allow identification and location of unwrapping errors, and evaluation of the quality of the ERS and ENVISAT time series. The latter is relevant since correct estimation of the time series is necessary for the methodology described in Section 4. We highlight that these maps are used as a global measure of quality of the unwrapping results.

Existing studies in the Mexico Basin (e.g., [10,11,12,13,40]) suggest that the main component of the deformation is vertical. Thus, the line of sight (LOS) velocities (v_LOS) were converted into vertical motion v_s using the simple geometrical expression v_s = v_LOS /cosθ, where θ represents the satellite look angle, and used in our post-processing methodology.

3.2. SBAS Results

3.2.1. Global ERS and ENVISAT Results

The ERS and ENVISAT SBAS rates in the LOS direction, overlaid in Google Imagery, are shown in Figure 3. More than 55,000 ENVISAT and 40,000 ERS SDFP pixels were identified in the study area with a mean point density of 52 ENVISAT SDFP pixels/km² and 44 ERS SDFP pixels/km². Overall, the SDFP pixel density in both cases tends to be higher in urban areas and decreases in vegetated or in cropping land areas. Note that the slightly larger ENVISAT frequency may allow a better identification of certain phase-stable pixels (e.g., corresponding to manmade structures), while ERS might be more sensible to identify those in vegetated areas.

The times series of the areas impacted by the highest subsidence rates are shown in Figure 4a,b, being rescaled to the first acquisition date in each case, although the master database was the image from December 1999 and December 2005 for the ERS and the ENVISAT, respectively. They indicate nearly constant ground subsidence rates in almost all the locations over the periods of study.

The localities of Ecatepec, Tecamac, Coacalco and Jaltenco present the largest subsidence rates. Areas in the municipality of Tepotzotlan were identified to be seriously affected by subsidence (~70 mm/yr in LOS) for the period 1999–2000 (see Figure 4a); however, for the period 2002–2010, the deformation rates decreased (~53 mm/yr in LOS; see Figure 4b); probably, a slight reduction in the water extraction rates may have contributed to the recovering of the groundwater levels. The rapid subsidence rates observed in these towns generally correlate with the population growth (Figure 4c). Nevertheless, the intensive water extraction intended for agriculture and cattle, such as in Jaltenco [41] and in Tepotzotlan [42] may produce the clay dehydration and subsequent subsidence in these sites.

Ecatepec, where the largest deformation occurs, is the highest populated city within the study area (see Figure 4c). The increment of subsidence rates in Tecamac and Coacalco might also have been potentiated by the demographic growth (i.e., increment in the water demand). Until 2000 the population in Coacalco was much larger than in Tecamac, corresponding to the observed mean deformation velocities which were higher in the former city for the 1999–2000 period. However, according to the plot in Figure 4c, the difference in population between them decreased, and in 2010, the population in Tecamac exceeded that of Coacalco by approximately 87,000 people, possibly manifesting as more rapid subsidence due to larger groundwater pumping rates to satisfy the local needs (see Figure 4a,b). Note that the location of the maximum subsiding points in these cities differs for the two different periods of study (see Figure 3), indicating that land subsidence evolves temporally and spatially, affecting new areas or becoming more critical. Considerable settling velocities (up to 40 mm/yr and 60 mm/yr, for the periods 1999–2000 and 2002–2010, respectively) were also detected in agricultural zones.

Figure 3. LOS deformation maps. (a) ERS; (b) ENVISAT. Time series for selected SDFP pixels (red stars) in each map are shown in Figure 4. Examples of time series presenting a non-linear component and corresponding to pixels Aq, Bq, Nq and Dq are presented in Section 4.1.

Figure 3. LOS deformation maps. (a) ERS; (b) ENVISAT. Time series for selected SDFP pixels (red stars) in each map are shown in Figure 4. Examples of time series presenting a non-linear component and corresponding to pixels Aq, Bq, Nq and Dq are presented in Section 4.1.

Figure 4. Displacement time series in LOS of SDFP pixels E, T, J, C and TY (Figure 3) corresponding to: (a) ERS and (b) ENVISAT data sets; (c) Population evolution for the period 1990–2010 [27]. Dashed lines in panels (a–c) indicate the corresponding tendency line.

Figure 4. Displacement time series in LOS of SDFP pixels E, T, J, C and TY (Figure 3) corresponding to: (a) ERS and (b) ENVISAT data sets; (c) Population evolution for the period 1990–2010 [27]. Dashed lines in panels (a–c) indicate the corresponding tendency line.

The histogram in Figure 5 presents the LOS deformation rates for the ENVISAT and ERS SDFP pixels. Most phase-stable pixels correspond to negative values while a small number to positive or almost zero values, indicating a prevailing and continuous subsidence in the area for the period 1999–2010. The smaller percentage of the SDFP pixels (~30% for both ERS and ENVISAT) apparently experiencing a relative uplift, mainly correspond to the adjacent highlands surrounding the lacustrine plain; a rebound effect, suggesting as well elastic behaviour of the aquifer, could explain these values. On the other hand, they can also be interpreted as a possible remaining effect from (stratified and turbulent) atmospheric artefacts, particularly affecting the signal coming from mountainous areas. These observations may also indicate that these far field velocities need to be better estimated.

ERS and ENVISAT RMS system misclosure maps, used to evaluate the unwrapped phase, are compared in Figure 6. In both cases larger unwrapping errors occur in steep mountainous areas prone to being affected by geometrical distortion or where subsidence gradients are large. However, from the RMS maps, unwrapping errors in the ENVISAT results are larger than in those obtained from the ERS data. A larger number of noisy pixels or higher deformation gradients in the former case (ENVISAT) may explain these results. Inspection of deviation maps of each ERS and ENVISAT interferogram indicates presence of significant discrepancies in some of them (particularly those with larger temporal and/or perpendicular baselines). Notwithstanding, both RMS maps suggest overall good performance of the unwrapping algorithm and that reliable deformation time series can be used in the methodology for identification of fracture-prone areas described in Section 4.

Figure 5. Histogram of the ERS and ENVISAT LOS velocities in northern ZMVM.

Figure 5. Histogram of the ERS and ENVISAT LOS velocities in northern ZMVM.

Figure 6. RMS system misclosure maps overlaying Landsat Imagery computed from the difference between the observed interferometric phase and the one reconstructed from the estimated phase delay by using: (left) ERS images, and (right) ENVISAT acquisitions.

Figure 6. RMS system misclosure maps overlaying Landsat Imagery computed from the difference between the observed interferometric phase and the one reconstructed from the estimated phase delay by using: (left) ERS images, and (right) ENVISAT acquisitions.

3.2.2. Seasonal Analysis

An analysis of the subsidence was conducted considering the rainy (June–October) and the dry (November–May) seasons in the ZMVM [43] and their possible influence on the subsidence behaviour. For this purpose, Tultepec and Tultitlán in the northern part of the ZMVM were investigated which present a linear and non-linear subsidence behaviour, respectively. Both places are characterized by expansive soils [44,45] and a high percentage of crop areas (i.e., Tultepec: 57%, Tultitlan: 30%) which allow an important recharge of the aquifers at these specific locations (conditions that may differ compared to other sites in the ZMVM). These sites are about 6 km separated so that similar subsidence is expected; however, closeness to Cerro Tultepec (natural recharge area) might affect the subsidence behaviour.

The initial and final subsidence values at each rainy/dry season were considered in order to account for a straightforward comparison between the “seasonal” behaviour and the global tendency (Figure 7). It is noticed that, in general, the subsidence at these locations continuously increases without a significant seasonal effect. However, the overall results at the end of each season exhibit slightly lower subsidence during the rainy season. When comparing the line of seasonal subsidence and the global tendency in the period of analysis, it is observed that, in general, at the end of the rainy season the subsidence is above the tendency, while in the dry season it is below (see Figure 7). This effect is smaller in Tultepec than in Tultitlán, and it is particularly observed for the period 2002–2006.

Figure 7. Seasonal time series analysis. The wet and rainy seasons are distinguished in white and blue, respectively. The coloured red area correspond to the area defined from the “seasonal tendency” (light brown line) and the global subsidence trend (dashed black line) of the data (violet line).

Figure 7. Seasonal time series analysis. The wet and rainy seasons are distinguished in white and blue, respectively. The coloured red area correspond to the area defined from the “seasonal tendency” (light brown line) and the global subsidence trend (dashed black line) of the data (violet line).

In addition to the possible aquifer recharge associated with the rainy season, water infiltration into existing discontinuities (especially during flooding events) may also explain the subsidence variations in the area of study. Further research is needed to estimate in more detail the possible seasonal subsidence effect in the entire area of study and to find out if a possible aquifer recharge or reduction of water exploitation from the aquifers might reduce the subsidence.

4. Post -Processing Tool for Soil Fracturing Zonation

4.1. Methodology for Fracturing Zonation Based on InSAR Time Series Modelling

One of the main hazards associated with subsidence is soil fracturing, which influences and even destroys important infrastructure. The accelerated demographic growth of the ZMVM involves an increasing water demand and improvement of the waste and potable water services. In northern ZMVM, important projects to overcome the limitations of the actual drainage system are being carried out (e.g., TEO). Moreover, the uncontrolled urban sprawl over natural recharge areas or those prone to geological risks, suggest the need for a better urban planning. Thus, detailed subsidence and soil fracturing zonation is relevant in the area to mitigate their effects.

Conventional in-situ methods such as Ground Penetrating Radar (GPR), Global Positioning System (GPS) or geodetic levelling would be very effective for detection and monitoring of areas prone to ground rupture in the ZMVM, particularly in its northern part where to our knowledge no previous detailed related studies have been performed. In fact, campaigns were executed in several regions in Mexico and in the south part of the ZMVM, showing their potential for soil fracturing identification (e.g., [9,12,46,47,48,49,50]). Nevertheless, their usage would be costly both in terms of time and economics for complete fracturing investigation over the entire region of interest. Moreover, the depth of penetration of the GPR is constrained by the presence of conductive clays or high conductivity pore fluid, and the data interpretation requires a skilled operator or investment in training.

On the other hand, InSAR time series can provide useful information about anomalies probably associated with processes such as fracturing activity or related to an elastic aquifer’s response (e.g., [51,52,53]). Nonetheless, a manual evaluation of the time series in the entire area result impractical considering the density of stable-phase pixels (more than 55,000 for an area of approximately 900 km²). Thus, a methodology for automatic detection of areas prone to fracturing in the north of the ZMVM based on the identification and modelling of the non-linear deformation is presented here (Figure 8).

As a first step to the proposed methodology, a database with the accumulated displacement for the period of study is built (see Figure 8). The time series of each SDFP pixel, which quality is previously assessed by the RMS maps, are analysed to detect whether or not they adjust to a linear deformation model. For that purpose, residuals between the linear model and the actual deformation behaviour inferred from the InSAR time series are computed as [54]:

$$\text{r}=\sqrt{(1/\text{M}){\displaystyle \sum _{\text{t}=1}^{\text{M}}{\text{(}{\varphi }_{\text{t}}-{\varphi }_{\text{t}}^{\text{mod}}\text{)}}^2}}$$

(1)

where M is the number of images, ϕ_t the estimated phase delays and ϕ_t^m°^d the deformation modelled as at + b.

Figure 8. Simplified flowchart of the methodology for soil fracturing zonation.

Figure 8. Simplified flowchart of the methodology for soil fracturing zonation.

The map of residuals is used to identify areas subjected to non-linear deformation. From this map pixels with high residuals were selected and their time series evaluated to identify the pattern of the anomalies observed (see Figure 9). In principle, pixels with r > 9 mm were observed to adjust a quadratic polynomial function; however, the non-steady component may also fit other functions (see also Section 5). Based on this and a preceding analysis [11,26], the non-linear deformation is modelled using a quadratic model and mapped by extracting for each SDFP pixel the term b from the following equation

$$at+{bt}^2+c$$

(2)

This so-called quadratic coefficient map is later used as an indicator of areas prone to potential fracturing risk. Subsidence evolution profiles are then constructed where notable acceleration (or deceleration) is found and probable fracture triggering mechanisms are evaluated. For the latter, a database containing various information (e.g., seismic activity, extraction rates, flooding) is organized.

Figure 9. Time series of points showing accelerated (a) and decelerated (b) motion patterns. Top: deformation adjusted with a linear model (solid coloured line) Bottom: residuals (coloured dots) between the linear fit and the observations, adjusted with a quadratic polynomial function (solid coloured line). See location of points Aq, Bq, Dq and Nq in Figure 3.

Figure 9. Time series of points showing accelerated (a) and decelerated (b) motion patterns. Top: deformation adjusted with a linear model (solid coloured line) Bottom: residuals (coloured dots) between the linear fit and the observations, adjusted with a quadratic polynomial function (solid coloured line). See location of points Aq, Bq, Dq and Nq in Figure 3.

4.2. Soil Fracture-Prone Zones and Impact on Infrastructure

Due to its more effective temporal sampling, the methodology for soil fracturing zonation was applied to the available ENVISAT data set. The quadratic coefficient map indicates the presence of notable acceleration and deceleration in several areas (Figure 10a). Overall observation of the surface geology (Figure 10b) and of Figure 10a suggests that areas prone to ground fractures may occur where contrasting materials are found and/or there is a high density of wells.

The largest acceleration zone is located in Ecatepec the Morelos, near the Gran Canal del Desagüe (GCD) and the TEO. Notable decelerated motion has been found in other places as well. The fracture-prone zone identified in Tecamac, for example, coincides with records of prominent subsidence and fracturing [55,56]. These phenomena impacted the fresh drinking water services and the drainage system [55]. A cavity-liked shape deformation identified in this area could have been produced by the presence of an underground cavern (volcanic cones have been mapped in this sector; [55]) or induced by water leaks in the hydraulic system beneath [57,58].

Larger acceleration is also observed near the Teoloyucan Pozos de Accion Inmediata (PAI) well branch. Spatially non-homogeneous distribution of the pumping rates might lead to complex and different stress-strain behaviours within short distances (see also [59]), thus making the area susceptible to fracturing. On the other hand, the presence of collapsible soil (i.e., vertisol) [60,61] and frequent flooding may interact with the differential settlement mechanism intensifying the propagation of fractures.

Figure 10. (a) Map of quadratic coefficient overlaid on SRTM V4 shaded relief. Three different zones, defined in terms of the geotechnical properties of the soil [76], are approximately represented in the study area: (i) Zone 1- Lacustrine, formed by clay deposits of high compressibility which are found within the area of the now-extinct lakes of the Valley of Mexico; (ii) Zone 2- Transition, formed mainly by alluvial deposits and (iii) Zone 3- Foothills, constituted by very compact but heterogeneous volcanic soils (mainly incompressible deposits). Names of main communities are depicted as well as the TEO vertical shafts and the profile E02-E02ʹ shown in Figure 11. Location of water treatment plants and waste disposal sites can also be observed; some of them close to areas prone to soil fracturing. (b) Geological formations (after [59]) related to zones 1, 2, and 3 in (a) and location of groundwater extraction wells.

Figure 10. (a) Map of quadratic coefficient overlaid on SRTM V4 shaded relief. Three different zones, defined in terms of the geotechnical properties of the soil [76], are approximately represented in the study area: (i) Zone 1- Lacustrine, formed by clay deposits of high compressibility which are found within the area of the now-extinct lakes of the Valley of Mexico; (ii) Zone 2- Transition, formed mainly by alluvial deposits and (iii) Zone 3- Foothills, constituted by very compact but heterogeneous volcanic soils (mainly incompressible deposits). Names of main communities are depicted as well as the TEO vertical shafts and the profile E02-E02ʹ shown in Figure 11. Location of water treatment plants and waste disposal sites can also be observed; some of them close to areas prone to soil fracturing. (b) Geological formations (after [59]) related to zones 1, 2, and 3 in (a) and location of groundwater extraction wells.

In Jaltenco, an area located in agriculture crops and where deformation appears to be accelerating has also been distinguished. Damage in buildings, railroads and hydraulic infrastructure due to fracturing has been reported to affect most areas where notable acceleration or deceleration occurs (e.g., [55,56,62,63,64,65,66]).

As stated by several researchers (e.g., [1,57,67,68]), ground ruptures associated with subsidence are the most disastrous; nevertheless, several mechanisms may interact producing their apparition and propagation. Particularly, the fracturing zone in Ecatepec de Morelos suggested by our methodology occurs near the abrupt transition zone described by Ovando Shelley et al. [69], where a sharp change in the clay thickness may trigger the occurrence of fractures. This area is also subjected to seismic activity [70] and influenced by frequent flooding. Temporal subsidence evolution profiles for the 2002–2010 period were constructed to further analyse the Ecatepec area.

The profile of the TEO in Figure 11, between the vertical shafts L0 and L5 is coincident with the location of the GCD. Ground subsidence has affected hydraulic infrastructure such as the GCD, completely reducing its capacity for wastewater discharge outside the city. Thus, the TEO has been constructed to overcome the limitations of the GCD and to increase the existing capacity of the drainage hydraulic system in the ZMVM. Nevertheless, continuous subsidence also impacts this infrastructure.

Figure 11. 3D temporal evolution deformation profiles across the Ecatepec de Morelos area. Profile E02-E02ʹ crosses several sites (dotted lines) reported to be affected by fractures and intersects Section L0–L5 along the TEO structure. Probable fracturing activity (PFA) may develop at the location of displacement jumps. Different materials are denoted with colours (yellow: mainly alluvial deposits; orange: primarily compressible deposits).

Figure 11. 3D temporal evolution deformation profiles across the Ecatepec de Morelos area. Profile E02-E02ʹ crosses several sites (dotted lines) reported to be affected by fractures and intersects Section L0–L5 along the TEO structure. Probable fracturing activity (PFA) may develop at the location of displacement jumps. Different materials are denoted with colours (yellow: mainly alluvial deposits; orange: primarily compressible deposits).

The maximum subsidence for the period of study (~1490 mm) along the TEO profile analysed is located close to the vertical shaft L0. From this zone, the subsidence decreases continuously up to the location close to the shaft L3. Major subsidence is particularly located between shafts L0 and L2, if compared to that between L2 and L5. The correlation between the accumulated subsidence in this profile and the thickness of the clay deposits [71] might explain the different subsidence variations between sections L0–L2 and L3–L5, thus becoming more important than the large pumping rates in the area. Several displacement jumps are found in this subsidence evolution profile that might be related to the location of possible fracturing areas, probably induced by changes in the clay deposits depth (smaller in section L2–L5 than in section L0–L2; see also Juárez Camarena et al. [72]). Furthermore, this profile suggests that probable fracturing activity (PFA) might have developed at approximately 2800 m from the vertical shaft L0 (between shafts L2 and L3), where contrasting displacement behaviour is identified. A concentration of pumping wells, tentatively compounding ground failures, is observed at this location.

The profile E02-E02ʹ across Ciudad Azteca (in Ecatepec de Morelos) and other areas, where the presence of soil fractures and differential subsidence was registered [73,74] is also analysed and described (see Figure 11). The subsidence in the profile section located in the Transition Zone (mainly alluvial deposits; see also Figure 10b) is not as significant as that in the lacustrine so that, at this boundary fracturing activity may generate. Some displacement jumps observed along profile E02-E02ʹ in the lacustrine zone (in orange) coincide with the location of reported fractures [74,75].

Note also that most of the inferred fracturing areas are located near waste disposal sites and/or water treatment plants (see Figure 10a). Thus, this map can also be used to improve waste and water management by proving critical spatial information for location of garbage dumps and treatment plants.

5. Discussion

Our ERS and ENVISAT results indicate the presence of extending subsidence in northern ZMVM for the period of study, whereby complex spatial and temporal distribution can be controlled by several factors acting or interacting. Most of the study area is located within the Cuautitlan-Pachuca aquifer zone, which is the main water supply source for the northern part of the Mexico Basin [15]. The sector located to the south of the Sierra de Guadalupe mainly corresponds to the Texcoco system. The maximum subsidence rates found in the latter area are in agreement with the approximate 290 mm/yr reported by Cabral-Cano et al. [77], despite the use of a different InSAR technique and period of analysis. The overexploitation of these aquifer units has produced dehydration of the compressible deposits (mainly clay and silt), reduction of the pore pressure and their compaction. The subsequent effect of this process is the surface subsidence. More than 70% of the water wells within the study area are located to the north of Sierra de Guadalupe (see Figure 1), where also the largest decline of the groundwater level has been recorded (~2.6 m/yr; [59]). From these observations, larger subsidence rates could have been expected to the north; nevertheless, opposed results have been found indicating that the magnitude of the land sinking in the entire area is not strictly controlled by the pumping rates and/or the drawdown of the piezometric levels (e.g., [78]), but by other parameters such as thickness of the compressible materials. Juárez Camarena et al. [72] characterized the materials of the northern ZMVM based on in-situ geotechnical measurements. They suggested that the presence of thicker and high water content soft soils decreases from the centre of plain (i.e., lacustrine environment) towards the surrounding rigid soils. Furthermore, they indicated that the depth of the lacustrine deposits varies decreasing from the Texcoco to the Xaltocan lakebed area, correlating with our findings about the spatial pattern of the subsidence. Larger pumpage could explain the observed rates in certain areas; however, the volume of water imported from outer Basins needs also to be considered to precisely evaluate if the ground sinking observed can be directly attributed to the local extraction practices. Besides, the presence of underneath paleo-channels or caverns in the area [14] can also produce sudden collapse.

As previously stated, the quadratic coefficient map indicates the presence of several areas potentially affected by fracturing. Nonetheless, precise identification of the triggering fracturing mechanisms is quite complex. In certain sectors, the presence of contrasting geological materials (e.g., near volcanic structures such as Sierra Guadalupe, Tultepec Hill; see Figure 10b) and the observed deformation rates suggest that differential deformation may be the main factor acting there. Continuous compaction of the lacustrine deposits surrounded by rocky soils reduces their thickness [69], thus compounding soil fracturing. Other agents might also add or interact such as hydraulic fracturing particularly after heavy rain events. A seismic component should not be precluded, since several epicentres have been identified in the area [79,80], and several houses have been reported to be severely impacted by the 2012 earthquake with an epicentre in Guerrero state [81]. Moreover, the reactivation of buried structures such as grabens [4,14] could also represent a potential source of fracturing and subsidence in the area.

It is clear that assuming a quadratic function for areas presumably affected by soil fracturing might not always be correct. Thus, we emphasize that more detailed geotechnical, geo-hydrological and geophysical information as well as field surveys (i.e., GPS and GPR) and reconnaissance of damages are still necessary to better define the threshold of the residuals r and to prove the assumption of quadratic deformation behaviour in all areas subjected to the apparition or propagation of soil fractures. On the other hand, non-linear deformation identified in some sites was observed to better adjust a polynomial function that describes a decoupled primary and secondary consolidation [82]. Nevertheless, the proposed quadratic approximation of the non-linear component is useful for a preliminary zonation and forecasting of possible hazardous fracturing areas. Accordingly, our methodology contributes to hazard management by ameliorating the vulnerability maps of the area.

The seasonal subsidence analysis needs to be extended to other areas to achieve more precise conclusions about the behaviour observed. In addition, because the deformation time series depends on the filtering settings, they need to be carefully evaluated as the selection of other parameters can lead to different interpretation of the results. As a remark, the effect of the rainy season is expected to be more significant for the southern area of the ZMVM than for the present study area due to the precipitation rates: <700 mm/year in the northern and 700–1000 mm/year in the southern and west area of the ZMVM [43]. However, comprehensive studies are required to verify this behaviour and, in general, to better understand the complex dynamics of the subsidence and associated fracturing in the area.

6. Conclusions

A detailed deformation analysis in the north of the ZMVM for the 1999–2010 period by means of multi-pass InSAR has been presented. Our results show increasing subsidence rates during the approximate 11-year-analysis, reaching up to 285 mm/yr. In general, highly populated areas are subsiding more rapidly than the rural zones, although considerable deformation has also been identified in the latter areas.

The exceptional soil characteristics of the ZMVM manifest in the complex subsidence and associated processes. Most land sinking in the study area (~90%) adjusts to a linear model; however, non-steady deformation was also identified. The proposed quadratic coefficient map suggests that areas where deformation appears to be accelerating or decelerating are affected by ground failures; these outcomes are in agreement with existing records. Moreover, subsidence affecting certain areas in expansive soils presents a particular behaviour: relatively larger ground sinking rates during the dry seasons in comparison to those corresponding to the rainy ones. Detailed evaluation and characterization of the subsidence in expansive soils are useful, particularly for the artificial recharge projects planned for the area and because they cause significant structural damage, affecting crop production as well [83].

Our interpretation and appraisal of ground rupture based on the InSAR technique can essentially contribute to geological risk management, providing valuable information for urban planning and strategic waste management, particularly over the northern ZMVM, where a few related studies have been undertaken. Moreover, the fracturing methodology herein proposed could be adapted to other areas and updated hazard mapping can be achieved by extending the deformation times series using data from on-going satellites such as TerraSAR-X, COSMO-Sky, RadarSAT -2 , Sentinel-1.