Wide-Area InSAR Survey of Surface Deformation in Urban Areas and Geothermal Fields in the Eastern Trans-Mexican Volcanic Belt, Mexico

Abstract

This paper provides the first wide-area Interferometric Synthetic Aperture Radar (InSAR) survey of the whole eastern Trans-Mexican Volcanic Belt (42,200 km²). The aims are to identify ground deformation hotspots within major urbanized areas and rural valleys, establish baselines in geothermal exploration sites, and analyze deformation at geothermal exploitation sites and its relationship with energy production. The whole 2003–2010 ENVISAT C-band SAR archive available over the region was processed with the Small BAseline Subset (SBAS) InSAR method to retrieve over 840,000 coherent targets and estimate their ground displacement rates and time series. Land subsidence hotspots due to aquifer drawdown are found within the city of Puebla (up to −53 mm/year vertical rates, groundwater pumping for industrial use), Tlaxcala and Apizaco (−17 mm/year, industrial and public), the valley of Tecamachalco (−22 mm/year, agricultural), Tulancingo (−55 mm/year, public, industrial and agricultural), and in the eastern Mexico City metropolitan area (−44 mm/year, agricultural). The baseline for the Acoculco caldera complex shows widespread ground stability. Conversely, localized subsidence patterns of −5 to −10 mm/year exist around Las Derrumbadas and Cerro Pinto in the Serdán-Oriental basin, due to intense groundwater pumping for agriculture. A well-defined land subsidence area with −11 mm/year maximum rates is found at Los Humeros volcanic complex within Los Potreros collapse, correlating well with energy production infrastructure location and historical steam production rates. Field surveys carried out in Acoculco and Los Humeros in 2018 provide supporting evidence for the identification of hydrothermal manifestations, and understanding of the landscape and surface deformation patterns within the geothermal fields.

1. Introduction

There is an increasingly growing literature on the use of remote sensing for geothermal exploration and research to gather geological, geophysical and geochemical knowledge on the geothermal resources to be potentially exploited (e.g., [1]). Remote sensing methods in such context exploit optical, thermal, hyperspectral and microwave data, the latter mainly based on satellite Synthetic Aperture Radar (SAR) imagery and advanced processing using differential Interferometric SAR (InSAR) methods. These can be based on either single-pairs of SAR scenes (e.g., [2,3]) or multi-temporal approaches to identify persistent scatterers or coherent targets and derive their ground deformation history with up to millimeter precision (e.g., [4,5]).

Mexico is among the top producers of geothermal energy globally [6], with 26 geothermal exploration permits in force as of the end of 2018, and six utilization concessions granted for the exploitation of geothermal resources and production of electricity, as reported by the Ministry of Energy [7]. Reflecting the geology and tectonic setting of the country, four of the concessions are located within the Trans-Mexican Volcanic Belt (TMVB) in central Mexico, i.e., Los Azufres and Los Humeros, and more recently Domo de San Pedro and Cerritos Colorados. The other two are in Baja California and Baja California Sur, i.e., Cerro Prieto and Las Tres Virgenes, respectively.

Many research studies based on InSAR targeted several sites and regions of interest in Mexico, to characterize land motion and instability, mainly due to exploitation of groundwater resources, volcanic activity and tectonics (e.g., [8,9,10,11,12,13,14,15,16]). Somehow rarer are, on the other hand, InSAR studies on the analysis and monitoring of surface deformation in Mexican geothermal fields. ERS-1/2 SAR imagery acquired in the 1990s over Cerro Prieto was used to model subsidence due to geothermal fluid withdrawal [17]. The same geothermal field was also the geographic focus of an InSAR study based on ENVISAT scenes acquired in 2003–2006 [18], another based on RADARSAT-2 data in 2011–2016 [19] and a more recent investigation exploiting Sentinel-1 imagery [20].

Los Humeros volcanic complex was surveyed with ENVISAT differential interferograms for the period 2003–2007 [21,22]. The same data were further processed using a single-master persistent scatterers InSAR method to model volume changes in the reservoir [23]. Other tests using Sentinel-1 imagery for the period 2016–2017 were also initiated [24]. These studies confirmed the presence of non-negligible topography-related tropospheric artifacts in the differential interferograms, which caused difficulties in the retrieval of deformation signals and needed to be compensated. More recently, Sentinel-1 pairs were used to model co-seismic deformation that occurred during the Mw 4.2 earthquake recorded within the caldera at 1.5 km depth on 8th February 2016 [25,26].

These initial InSAR observations at Los Humeros demonstrate that this geothermal field offers a scientifically interesting test site to investigate the potential intermingle between ground deformation patterns and local and regional tectonics, seismicity, and energy production activities. An investigation scenario as articulated as that of Los Humeros, also suggests that InSAR observations should be considered even prior to the exploitation of geothermal resources, in order to derive a surface deformation baseline that can provide information on the status of the site before any activity initiates. This is the case of the Acoculco, Cerro Pinto and Las Derrumbadas geothermal fields, within which exploration was approved to investigate their geothermal potential [27,28,29].

The appraisal of the environmental scenario of ground stability in areas of authorized geothermal exploration contributes, indeed, to collecting state-of-the-art knowledge about the landscape and the geohazard status of the sites of potential exploitation. Exploration typically consists of a series of geological, geophysical and geochemical investigations aimed at a better characterization of the geothermal field and corroboration of the existence and potential of the geothermal resource. In such context, ground displacement data from InSAR can feed into the exploration process by providing information on wide-area ground stability (or instability), and presence of any apparent ground deformation phenomena occurring at the local scale, either connected or not with the presence of the geothermal resource.

In the specific case of the above mentioned Mexican geothermal exploration sites, the environmental scenario of ground stability cannot be retrieved reliably without accounting for the wider geological context of the TMVB. The InSAR analysis needs to start from a regional scale approach, aiming to identify first the presence of any deformation patterns due to wide-scale processes that may represent a background deformation field. Secondly, such regional analysis allows an inventory of local scale deformation patterns due to other anthropogenic activities (e.g., groundwater pumping), not to be confused with those due to geothermal exploitation. This approach suits well the case of the eastern TMVB, a region where different geological processes and human activity coexist, given the presence of extended agricultural areas, one of the largest urban areas of Mexico (i.e., the city of Puebla), and numerous minor villages and towns scattered across the valleys. Among all the processes, subsidence due to groundwater pumping is renowned to be widespread across central Mexico (e.g., [12,15]).

Therefore, to achieve an improved understanding of land deformation processes occurring at Los Humeros in response to historical steam production rates in the 2000s, and at the geothermal fields of Acoculco, Cerro Pinto and Las Derrumbadas still under exploration, in this work a multi-temporal InSAR processing approach is used to process the whole archive of ENVISAT C-band Advanced SAR (ASAR) imagery. The study region encompasses the whole eastern sector of TMVB, for a total of 42,200 km² (Figure 1), extending from the state of Hidalgo, through the eastern sector of México and Morelos, including the state of Tlaxcala, a large portion of the state of Puebla and the western periphery of the state of Veracruz (Section 2).

By exploiting a number of interferometric pairs and their temporal combination, a robust estimation of ground deformation occurring in the region is achieved (Section 3). Significant deformation rates are found within the urban footprints of Puebla, Tlaxcala and Apizaco, and the valleys of Tecamachalco and Tulancingo, where groundwater abstraction for public, agricultural and industrial use induces ground compaction at significant annual rates (Section 4.1). Ground stability is then investigated at Acoculco, Cerro Pinto and Las Derrumbadas, with the goal to characterize for the first time the baseline surface deformation scenario of the Tulancingo-Acoculco caldera complex and the Serdán-Oriental basin, respectively (Section 4.2). The InSAR results are finally discussed in relation to the energy production history and infrastructure distribution at Los Humeros, considering also natural factors such as topography, morphology, geology, and seismicity (Section 4.3). The conclusions outline future research directions on the use of InSAR to achieve a thorough geological interpretation of subsidence mechanisms observed in the eastern TMVB, and particularly at the exploitation site of Los Humeros (Section 5). Among future activities, Copernicus Sentinel-1 SAR data will be used to extend the baselines and the assessment to the 2010s and up to present.

2. Study Area

Spanning east-west across central Mexico, from the Gulf of Mexico to the Pacific Ocean, the 1000 km-long, Plio-Quaternary continental volcanic arc of the TMVB rests on the southern edge of the North American Plate, beneath which subduct the Rivera and Cocos oceanic Plates along the northern end of the Middle American Trench (inset in Figure 1, and Figure 2).

The area of investigation lies within the eastern sector of the TMVB and is bounded by Sierra Nevada to the west, Sierra Madre Oriental to the north, and Sierra Madre del Sur to the south. The area extends from the town of Pachuca in the north, towards the city of Puebla, the natural reserve of Sierra de Tentzo, and up to the town of Acatlán de Osorio to the south. Terrain elevation generally ranges between 1000 and 2500 m a.s.l. (Figure 1), with peaks at volcanoes such as Popocatépetl (5426 m a.s.l., active), La Malinche (4461 m a.s.l., inactive) and Cofre de Perote (4282 m a.s.l., inactive).

Land cover is dominated by agricultural fields, although significant is the proportion of desert scrub and pasturelands, especially in the eastern sector of the investigated region. Forested areas are mainly found approximately between 2600 and 4000 m a.s.l. over the Sierra Nevada, La Malinche, Sierra Madre Oriental and Sierra Madre del Sur (Figure 1). Excluding the small portions of the densely populated metropolitan area of Mexico City, the largest urban area within the surveyed region is the city of Puebla (550 km², 1.5 million inhabitants), capital of the homonymous state, followed by the city of Pachuca (150 km², 0.27 million inhabitants) in the state of Hidalgo. Many minor villages and towns also concentrate in the valleys.

Figure 2 provides a 1:250,000 scale geology of the region, while the local geological setting of each site is described in Section 2.1 for Los Humeros caldera complex, Section 2.2 for Las Derrumbadas complex, and Section 2.3 for Tulancingo-Acoculco caldera complex. For the sake of brevity, and to avoid duplication with published literature, the geological settings in this section are intentionally not accompanied with detailed maps (the reader can refer to the cited references). Section 2.4 also describes the valleys of Puebla and Tecamachalco, to pave the way for the discussion of the InSAR results in the context of groundwater exploitation and aquifer drawdown. Further information about aquifer composition and monitoring is also provided in the context of the data discussion across Section 4.

2.1. Los Humeros Caldera Complex

Los Humeros volcanic complex, located in the Serdán-Oriental basin, has been the focus of a number of geological, geochemical and geophysical studies since the 1960s, due to its geothermal potential (e.g., [34,35]).

In a sector of the TMVB with stronger alkaline influence, at the northern border of the graben of Libres-Oriental, Los Humeros caldera is a calc-alkaline volcanic complex, with an origin associated with the subduction of the Cocos Plate [36]. The regional basement of the caldera emerges in the Sierra Madre Oriental (Figure 1), formed by Paleozoic metamorphosed granodioritic rocks, covered in discordance by a Mesozoic sequence of the Tethysian type, where the continental opening facies from the Triassic-Jurassic are clear, followed by a sequence of sediments and precipitates of mixed origin, composed of gypsum and terrigenous sequences, sometimes carbonated (e.g., [37]). The Tethysian Sea with marine facies was established (Jurassic to Upper Cretaceous), interrupted by the Laramide deformation and later by granodioritic Miocene intrusions. The sequence is covered by Pleistocene andesites and basalts, and voluminous pyroclastic deposits produced by Los Humeros and covering most of the Serdán-Oriental basin (Figure 2).

The current morphology of the caldera complex is the consequence of two main collapses, the first caused by an explosive rhyolitic eruption (~0.46 Ma) that formed the quasi-circular 15–20 km-wide Los Humeros caldera, which was followed by the eruption of the Zaragoza ignimbrite (~0.1–0.14 Ma) that allowed the development of the 7–10 km-wide Los Potreros caldera, nested within Los Humeros caldera (e.g., [37,38,39]). The much smaller Xalapazco crater also formed (~0.02–0.04 Ma) within Los Potreros caldera, followed by more recent minor eruptions and production of extensive, morphologically youthful, basaltic lava flows (see Section 4.3 for a detailed geological map of the geothermal field).

The structures of Los Humeros caldera can be grouped into three large sequences: the regional structures that manifest in the caldera, the caldera rim structures, and the structures that developed inside the caldera. The regional features show larger annular structures, especially in the north-eastern, eastern and western sectors of the caldera. Within Los Potreros caldera, where the geothermal field activity takes place, the intra-caldera structures concentrate asymmetrically, mainly in the southwestern sector. The main reservoir of the geothermal field is hosted by the old volcanic succession of pre-caldera Cuyuaco, Alseseca and Teziutlan andesites (10.5–1.55 Ma) of low to medium permeability, which are in the form of several blocks mostly separated by normal faults, underlain by the low-permeability basement of Cretaceous-Jurassic limestone, laterally confined by the calderas ring faults, and sealed upwards by low-permeability Quaternary ignimbrites [38,40].

Two main systems of faults can be distinguished, i.e., a NE-SW to E-W striking and north-dipping faults (e.g., Las Papas and Las Viboras), and another with NW-SE to N-S orientation and west-dipping including Los Humeros and Maxtaloya faults. In the southern sector of the caldera, there are two NE-SW and NW-SE oriented eruptive fractures that can be considered as regional structures due to the alignment of monogenetic centers or effusive activity, which concentrate the most recent and basic (sometimes alkaline) volcanism of Los Humeros caldera (e.g., [41]). Recent field observations provided evidence of neotectonic deformation of the caldera floor and movements along the volcano-tectonic faults inside Los Potreros caldera, likely originated by an active resurgence of the caldera floor [40].

Los Humeros is currently the third-largest geothermal field in Mexico, after Cerro Prieto and Los Azufres [42]. Studies suggest that the magma chamber is now partially solidified, with present-day temperatures of 600–650 °C at 5 to 7–8 km depth (e.g., [35]).

The Mexican Federal Commission for Electricity (CFE) started exploration in 1978, drilled the first well in 1981, and began operations in 1991 with the first geothermoelectric generation unit operating at 5 MW [33]. Since then, the infrastructure significantly developed and now the geothermal field comprises over 50 wells for production, injection and exploration, and 11 generation units with a total installed capacity of 93.6 MW and operating capacity of 68.6 MW [43]. Production wells are mainly located along the main active faults NNW-SSE (e.g., Maxtaloya-Los Humeros faults) or near fault fracturing zones in the northern sector of Los Potreros caldera [40]. Works are ongoing for further development of the energy production infrastructure, with perforation of some new wells and construction of new geothermoelectric units, as part of the Los Humeros III project [44].

2.2. Las Derrumbadas Complex

At ~30 km to the south of Los Humeros and 10 km to the east of Totolcingo endorheic lagoon, the rhyolitic twin-dome volcano complex of Las Derrumbadas overlooks the Serdán-Oriental basin (Figure 2). The complex consists of two domes, i.e., the Derrumbada Roja or Derrumbada de Fuego (southern dome, 3480 m a.s.l.) and Derrumbada Azul or Derrumbada de Agua (northern dome, 3420 m a.s.l.), surrounded by debris-avalanche deposits and expanses of lahar and pyroclastic deposits [45]. The Pleistocene rhyolite tuff ring-dome complex of Cerro Pinto or Derrumbada Blanca (2992 m a.s.l.) [46], located only a few kilometers north, is often associated with and included in Las Derrumbadas complex (see also Figure 1 and Section 4.2.2).

Las Derrumbadas domes are emplaced in a regional NW-SE lineament, interpreted as a front of overthrust affecting the marine sedimentary sequence [35]. Together with Cerro Pizarro, these domes are the most prominent features across the Serdán-Oriental basin (e.g., [47,48]).

The permission for geothermal exploration at Las Derrumbadas was granted to the company Geotérmica Derrumbabas S.A. de C.V. [28]. The area involved is 150 km² wide, and extends from the Alchichica lagoon to the north, through the villages of Techachalco, San Luis Atexcac and up to Los Álamos Tepetitlán to the south. Instead, the exploration permission area granted to the company Geotérmica Cerro Pinto S.A. de C.V. for Cerro Pinto extends 149 km² [29], from the town of Tepeyahualco to the north, towards the Derrumbada Roja to the south, and includes the village of Jalapasco (see Section 4.2.2 for detailed delineation of the exploration sites). High enthalpy, deep geothermal systems that controlled the location of eruptions and the distribution of heat sources exist [28,29,35]. Thermal water springs and other superficial geothermal manifestations (e.g., fumaroles and hydrothermally altered rock) have been identified in the area, confirming the geothermal potential of this field.

Las Derrumbadas complex is located within the Libres-Oriental aquifer, extending 3500 km² from Los Humeros caldera to the north, towards the south up to La Malinche volcano. The basin is endorheic, and surface water flow is limited. The main input comes from the west from Tlaxcala through the Barranca La Malinche river that flows into the Totolcingo lake (or El Carmen), and there is no outflow to external water bodies [49]. The aquifer mainly consists of extrusive igneous rocks, such as pyroclasts and lavas, except in the lower sector, where its upper section includes unconsolidated alluvial deposits, and in the mountain massifs, where it is constituted by fractured lava effusions. In some sectors, such as at the Totolcingo and Tepeyahualco (or El Salado) lakes, the aquifer is semiconfined by clayey deposits of low permeability and thickness of various tens of meters. In other zones, the aquifer can be locally semiconfined by fine-grained alluvial deposits.

2.3. Tulancingo-Acoculco Caldera Complex

In the northern sector of the state of Puebla, the Tulancingo-Acoculco calc-alkaline caldera complex shares the same high-level regional tectonic regime of many other volcanic complexes of the TMVB, with origin associated with the Cocos Plate subduction beneath the North American Plate.

Superficial lithology is dominated by a volcanic sequence of rhyolites, dacites, andesites, basalts, cineritic domes, and cones, onto faulted calcareous sediments composing the Cretaceous basement of the Sierra Madre Oriental (Figure 2), which is exposed to the southeast of the caldera complex close to the towns of Chignahuapan and Zacatlán, and to the north of the municipality of Tulancingo de Bravo (see also Section 4.2.1). As the Tulancingo–Acoculco caldera complex lies at the intersection of two regional fault systems, both NW–SE and NE–SW-trending regional systems can be observed, i.e., the Tulancingo–Tlaxco and the Apan–Piedras Encimadas fault systems.

The volcanic complex started to develop during late Pliocene with formation of the 32 km-wide Tulancingo caldera structure (3.0–2.7 Ma), within which formed the 18 km-wide Acoculco caldera during Pleistocene (1.7–0.24 Ma) [50]. Volcanic activity persisted until 0.06 Ma, with several domes, cinder cones, fissure lava flows and two ignimbrite eruptions (1.2–0.65 Ma) [51].

The Acoculco volcanic complex is considered to be a super field of high enthalpy. Pervasive hydrothermal alteration, meter-scale craters formed by hydrothermal explosions, and mounds of hydrothermal debris and breccias, as well as gas discharges and acid-sulfate springs of low temperature, can be observed [52].

The permit granted to CFE for the site of Acoculco involves a 150 km² area [27], east of the homonymous town. Exploration started in the 1980s and is still ongoing. Two deep exploratory wells were drilled by CFE in 1994 and 2008, revealing temperatures above 300 °C at 2 km depth—which suggest the site as a potentially interesting location for exploitation—but absence of fluids and low permeability [44], reason for which the complex has been classified as a hot-dry rock reservoir. The magmatic heat supply, magma reservoir, and processes at Acoculco were recently investigated in [53], to better characterize the potential of the geothermal field.

2.4. Puebla and the Tecamachalco Valleys

The fifth most populated city of the whole Mexico, Puebla (or also, Heroica Puebla de Zaragoza, or Puebla de los Ángeles) is the largest urban area of the valley of Puebla, with a population of 1.5 million inhabitants in 2015 [54], or 2.9 million if considering the larger agglomeration of the metropolitan area of Greater Puebla. Its historic center was inscribed on the UNESCO World Heritage List in 1987 for its major religious buildings such as the 16th–17th century cathedral, fine buildings like the old archbishop’s palace, and typical houses with walls covered in colorful tiles (azulejos). With ~10 million of Mexican and foreign visitors recorded in 2018 only [55], Puebla is one of the most touristic cities of Mexico, and this constant presence sums up with permanent residents.

The Cuetlaxcoapan valley, where the city was founded in 1531, is surrounded by the peaks of Sierra Nevada (i.e., Popocatépetl, Iztaccíhuatl and Mount Tlaloc) to the west, and La Malinche to the east (Figure 1), which act as the main recharge areas for the aquifer system, with annual rainfall of ~1000 mm.

A highly permeable, unconfined, shallow aquifer of Plio-Quaternary volcanic sands and gravels (130 m thick), overlies a low permeability semi-confined aquifer made up of Pliocene lacustrine deposits onto fractured Neogene andesite, pyroclastic deposits and conglomerates (up to 200 m thick), and Mesozoic carbonate rocks [56]. A deeper aquifer system is also present below the latter, and consists of Cretaceous limestone, overlying gypsum and anhydrite. The thickness of this sequence is controlled by the regional structures, such as the E-W Tlaxcala and Tetlatlahuca faults delimiting the Iztaccíhuatl-Malinche graben, and the NE–SW Malinche fault crossing the city of Puebla (Figure 2).

The whole aquifer of the valley of Puebla extends 2025 km² (Figure 3) and, as of 2010, included 1200 inventoried groundwater extraction sites, 735 of which wells, with primary use for public-urban needs and the total extracted volumes of 327 million m³/year recorded in 2014 [57].

A water deficit of about 700 l/s (i.e., 22 million m³/year) was estimated for the valley in 2006 [58], followed by an inversion of the trend and average annual availability of 18 million m³/year later in 2008 [59], and 24 million m³/year in 2015 [60]. One among the aquifers with denser network of groundwater wells of central Mexico (Figure 3), the unconfined Puebla valley aquifer system acts, indeed, as the primary source of potable water for the whole basin. Agricultural activities, population growth, and touristic presence contribute to increasing demands for hydric resources. Over the last 20 years industrial activities also developed significantly due to decentralization from Mexico City, and the overexploitation of aquifer resources—especially within the city of Puebla (where a number of new wells were drilled, and water was pumped at intensive extraction rates)—changed the regional hydrology and historical patterns of groundwater flows [61]. As a consequence, the need for groundwater quality monitoring and hydro-geochemical studies became stronger and stronger [56].

To the east of the aquifer of Puebla, the National Water Commission (CONAGUA) identifies another aquifer system with intensive water extraction activities, i.e., the valley of Tecamachalco (Figure 3) and the Esperanza-Oriental hydrological system (located between Pico de Orizaba and La Malinche volcanoes), which recently experienced a drop in groundwater levels [64]. Local stratigraphy consists of a Plio-Quaternary shallow aquifer of unconsolidated carbonate-rich volcano-sedimentary units overlying a 200 m thick unit of consolidated conglomerates that confines a deep aquifer of Cretaceous limestone, shale, and sandstone (Figure 2).

In addition to the supply from the Manuel Ávila Camacho reservoir (also known as Valsequillo or Balcón del Diablo reservoir) south of the city of Puebla, since 1975 water demands for irrigation in the valley of Tecamachalco started to be addressed by abstracting groundwater from the aquifer. With the declining water levels, extraction depths progressively increased [65,66]. Aquifer recharge takes place primarily via rainfall infiltration and horizontal underground flow from the highlands at the north and west of the basin, part from irrigation itself, and part from leaks in the distribution pipes, but there is no direct input from major rivers [66].

An increase in the extracted volumes of groundwater was observed from 1988 (228 million m³/year), to 1999 (283 million m³/year) and 2011 (343 million m³/year), with the majority of pumping occurring in the Distrito de Riego (i.e., west of Tecamachalco town, see sector with denser network of wells within the Tecamachalco aquifer in Figure 3), and minor in Esperanza, Tepeaca and Palmar de Bravo [66]. Over 80% of the aquifer exploitation is for agricultural purposes, 17% for domestic use, and less than 1% for industrial use [66]. A water deficit of 32 million m³/year was estimated in 2002 [67], rising up to over 68 million m³/year in 2008 [59]. In the framework of the yearly nationwide assessment carried out by CONAGUA, in 2018 the aquifer was still in deficit and also classified as overexploited [68].

3. Materials and Methods

3.1. Satellite SAR Data

The whole available archive of Advanced SAR (ASAR) Image Swath IS2 raw scenes acquired in C-band (i.e., 5.6 cm wavelength, 5.3 GHz frequency) by the ENVISAT satellite over the region of interest was used. These are characterized by a Line-Of-Sight (LOS) with ~23° look angle at scene center, VV co-polarization, ~20 m ground resolution and nominal site revisit of 35 days. A summary of the six data stacks that were processed is provided in Table 1.

In particular, twenty-two acquisitions along descending track 212 covering a period of four years between 8th April 2003 and 13th March 2007 were used to study the eastern sector of the region of interest (i.e., ID #1 and #2 in Table 1). The two processed stacks of standard frames cover overall an area of approximately 20,000 km², spreading from the city of Teziutlán in the northern sector of the state of Puebla, towards the south up to Tehuacán (see also Figure 1). A 7 year-long data stack acquired along descending track 484 between 16th February 2003 and 21st March 2010 (i.e., ID #3 and #4 in Table 1) was then used to cover the western sector of the region of interest, including the Acoculco complex, for a total extent of the processed area of around 21,100 km², from the village of Tulancingo in the north, towards the city of Puebla, and up to the town of Acatlán to the south.

Additionally, another seven year-long data stack acquired along ascending track 148 between 24th January 2003 and 2nd April 2010 (i.e., ID #5 and #6 in Table 1) was used to cover the western sector of the investigated area, with an extent of approximately 21,200 km², from the city of Pachuca in the state of Hidalgo, towards the south up to the town of San Juan Ixcaquixtla in the state of Puebla (see also Figure 1). It is to be noted that no data stacks populated with sufficient numbers of scenes acquired along ascending orbits were available for the eastern part of the area of interest, which was covered by only six images along track 105 during the whole ENVISAT mission lifetime.

3.2. Multi-temporal InSAR Analysis

After focusing the ENVISAT raw (level 0) scenes to single look complex data and precise co-registration of each stack to a common geometry, differential interferograms were formed for each of the six stacks, based on multi-master SAR pairs characterized by small baselines (i.e., a maximum of 400 m for the normal baseline, and 4 years for the temporal baseline), in order to minimize the effects of phase decorrelation (e.g., [69]). Phase noise was also reduced by employing 20 by 4 multi-looking in the azimuth and range directions, resulting in square pixels of around 90 m on the ground.

The multi-temporal InSAR image processing workflow that was used to identify coherent targets and estimate their annual motion velocities and temporal histories of deformation along the satellite LOS, followed the Small BAseline Subset (SBAS) approach as developed by [5], and recently parallelized by [70,71,72]. In particular, the Parallel SBAS (P-SBAS) algorithm integrated on ESA’s Grid Processing On Demand (G-POD) platform and its high-performance computing resources were exploited [73]. This implementation was successfully adopted in InSAR applications aimed to map and monitor natural and anthropogenic processes and hazards (e.g., [14,74,75,76,77]), and is particularly suited to survey wide regions with multiple data stacks [71], that can be combined together to achieve a comprehensive overview of surface processes affecting the study regions.

Using a temporal coherence threshold of 0.7, more than 840,000 coherent targets were identified overall from the analysis of the six stacks of images (Figure 4 and Table 1), corresponding to average densities of up to 21 targets/km² across each image frame. As expected, the distribution of coherent targets throughout the region varies with land cover. Greater densities of targets are found over urbanized areas (e.g., Puebla, Tlaxcala, Tulancingo, see Section 4.1) as well as onto rocky terrain and areas covered by desert scrub (e.g., onto basaltic andesite south of Los Humeros caldera, see Section 4.3), where the interferometric coherence as observed in the differential interferograms was generally high. On the other side, due to temporal decorrelation much sparser distributions are found in rural zones, such as agricultural lands and pasture lands (except at places where urban infrastructure is found). The large expanses of agricultural lands and forests explain the somewhat lower density of targets per square kilometer found in some sectors of the investigated area. No coherent targets are observed along most of the elevation of the volcanic complexes such as Sierra Nevada and Pico de Orizaba, due to the presence of dense forested covers and snow at the highest elevations, and consequent phase decorrelation.

The reference points to which all the InSAR-derived deformation estimates were referred were set at locations deemed to be geologically stable (see Figure 4), specifically at:

• 19.562 N, −97.244 E within the city of Perote in the Mexican state of Veracruz for data stack #1

• 18.872 N, −97.712 E onto limestone to the east of the town of Tecamachalco for stack #2

• 19.285 N, −98.432 E within the town of Texmelucan (25 km north-west of Puebla) for stack #3

• 18.912 N, −98.428 E within the town of Atlixco (15 km south-west of Puebla) for stacks #4 and #6

• 20.124 N, −98.738 E within the city of Pachuca for stack #5.

The displacement rates of the identified coherent targets reach negative values up to −51 mm/year as estimated along the ENVISAT LOS, indicating deformation occurring in the direction away from the satellite sensor. While modal values for the velocities observed within the whole area are approximately 0 for each of the six processed stacks, evidence of clear negative skewness (i.e., right-modal or left-skewed distribution) is found in the results, due to the occurrence of ground subsidence at several locations, as discussed in detail in Section 4. This is very apparent for the two stacks covering the city of Puebla (i.e., ID #4 and #6 in Table 1), where across each processed frame displacement rates are negatively skewed with average values of −2.7 and −3.0 mm/year. Such values reflect the severe subsidence rates occurring in the city (Section 4.1.1).

The uncertainties (precision) associated with the annual velocity estimates are mainly controlled by the number of small baseline interferograms that were identified and exploited during the processing of each of the six data stacks. Their values can be computed using the empirical relationship proposed in [78]. In particular, the LOS velocity standard errors for the six datasets of results are: 1.5 mm/year for datasets ID#1 and ID#2 given the limited number of input scenes available (i.e., 22 scenes, combined into 56 small baseline interferograms, Table 1), and 1.1–1.2 mm/year for ID#3 to ID#6, thanks to the greater amount of scenes processed and consequent great number of interferograms formed (i.e., 29 to 34 scenes, combined into 84 to 97 small baseline interferograms). It is by accounting for these uncertainties that the InSAR results are interpreted in Section 4.

Cross-comparison of the InSAR results for the area of overlap (~3500 km²) between the two stacks in descending mode (i.e., tracks T212 and T484) shows the agreement of the estimated annual rates for coherent targets that were identified in both InSAR analyses. This outcome is relevant, considering the difference in temporal periods covered by the two input SAR data, i.e., 2003–2007 and 2003–2010, and thus confirms on one side the suitability of the selected reference points as good locations to tie all ground displacement observations to and, on the other side, the persistence of similar annual rates in the period 2007–2010. No adjustment between the datasets of InSAR results from the two descending tracks was therefore deemed necessary prior to their geological interpretation and further analysis. The same evidence was found in the areas of overlap between the processing results from the descending and ascending mode stacks (i.e., tracks T484 and T148, respectively). Further methodological insights into the combination of multiple datasets of SBAS InSAR results for large-area surveys, and calibration to a unique reference datum using GPS measurements is provided, for instance, in [71].

3.3. Estimation of Vertical and Horizontal Displacement Rates

When information from both ascending and descending datasets was available for the same area, for instance for the municipality of Puebla (see Section 4.1.1), the InSAR results providing estimates along the sensor LOS in the two geometries (i.e., V_A and V_D) were combined to derive vertical and horizontal components of the ground displacement vectors. A comprehensive methodological reference to implement such multi-geometry combination is provided, for instance, in [79].

In particular, the combination was possible by assuming absence of north-south displacement components (i.e., V_N = 0), and by considering the east-west, north-south and vertical directional cosines (i.e., the cosines of the angles that the vector forms with the three coordinate axes, namely E, N and U) of the LOS in the ascending and descending modes, i.e., E_A, N_A and U_A, and E_D, N_D and U_D, respectively.

In order to derive the components along the vertical and east-west directions of the ground displacement vectors, i.e., V_U and V_E, the system of two equations:

V_A = U_A*V_U + E_A*V_E + N_A*V_N (1)

V_D = U_D*V_U + E_D*V_E + N_D*V_N (2)

was solved for V_U and V_E under the assumption that V_N = 0, and using the directional cosines of the ENVISAT LOS. For the datasets used in this work, the latter are −0.3791, −0.0945, +0.9205 for the ascending mode, and +0.3791, −0.0945, +0.9205 for the descending mode.

Consequently, the system of two equations above can be simplified as follows:

V_U = 0.543*(V_D + V_A) (3)

V_E = 1.319*(V_D − V_A) (4)

It is worth noting that the processing of descending and ascending mode data stacks leads to the identification of different targets on the ground, depending on their geometrical characteristics and orientation with respect to the satellite LOS. The number of processed scenes for each stack and their baseline distribution are factors that could influence the resulting distribution and coverage of coherent targets across the processed area. To account for this different spatial distribution, the combination of descending and ascending datasets was preceded by the resampling of the point-wise results into a raster with regular grid spacing of 50 m. Subsequently, the combination of descending and ascending results was implemented only for those pixels including InSAR results from both geometries.

In other circumstances (e.g., outside the overlapping area between ascending and descending mode results), only one observation geometry—either ascending or descending—was available. In those cases, in order to solve Equations (1) and (2), the more strict assumption of absence of any horizontal displacement component—either north-south or east-west—was used, i.e., V_N = 0 and V_E = 0. This was, for instance, the case of the Tulancingo valley (see Section 4.1.4) and the sites within the Serdán-Oriental basin (see Section 4.1.3, Section 4.2.2 and Section 4.3), where data along only one geometry was available.

Consequently, under the assumption that V_N = 0 and V_E = 0, and when only velocity information from the ascending mode was available (i.e., V_A), vertical velocities were computed as follows:

V_U = 1.086*V_A (5)

or, analogously, based on the V_D, when descending mode results were available only:

V_U = 1.086*V_D (6)

The same approach was used to covert LOS displacement histories into vertical for targets of interest, for instance for the time series of targets within Los Humeros caldera (see Section 4.3).

It is worth noting that, in zones where significant horizontal displacement is present, the assumption of null horizontal deformation may introduce a bias in the interpretation of the results, with errors depending on the incidence angle and the actual amount of ongoing horizontal deformation (e.g., [79]).

4. Results and Discussion

The region surveyed with InSAR covers a total of 42,200 km². Apparent spatial patterns indicating the occurrence of ground deformation were observed in the areas of Puebla, the valley of Tecamachalco, Tlaxcala, Apizaco, Tulancingo, Las Derrumbadas and Los Humeros (Figure 4).

The identified areas of ground displacement are discussed in the following sections, focusing on three main topics: land subsidence and groundwater management (Section 4.1), identification of ground stability and instability baselines in areas of geothermal exploration (Section 4.2), and analysis of the relationship between surface motion and geothermal energy production (Section 4.3).

4.1. Land Subsidence and Groundwater Management

Many of the observed ground displacement patterns across the surveyed region reflect surface deformation due to groundwater pumping from the shallow and deep aquifers for public, agricultural and industrial use, and consequent water level drop and aquifer depletion, as later discussed in this section for the areas of Puebla, the valley of Tecamachalco, Tlaxcala, Apizaco, Tulancingo and the eastern Mexico City metropolitan area.

Adopted in 1992, the National Waters Law establishes that groundwater pumping concessions in Mexican aquifers should take into account the average annual availability of groundwater, which is estimated by CONAGUA and published in Mexico’s Official Gazette of the Federation (Diario Oficial de la Federación, DOF, e.g., [59,60]), at least every three years. This is of crucial importance since many changes in the annual availability of water for several aquifers in Mexico were observed as a consequence of modifications in the natural water recharge and discharge balance, due to water pumping at many locations.

In 2008 CONAGUA estimated a deficit of a few to up to several tens of millions of cubic meters per year for many of the aquifers within the investigated area. These include the valley of Tecamachalco, the valley of Tulancingo and Libres-Oriental, for a total of 68, 21 and 13 million m³/year of groundwater deficit, respectively. In addition, at the western margin of the region of interest, the aquifers of Cuautitlán-Pachuca, Texcoco, and Chalco-Amecameca also registered deficits of 190, 49 and 16 million m³/year, respectively [59]. Five of these six aquifers were still in deficit as of 2011 (Figure 3). All these locations are where the most pronounced ground displacement patterns are observed based on the InSAR results (Figure 4).

4.1.1. The Valley of Puebla

The most remarkable and apparent feature of identified ground displacement in the InSAR results is found in the relatively-flat young deposits (mainly ravines and alluvial plains) between Sierra Nevada and La Malinche, within the aquifer of the valley of Puebla (see Figure 1, Figure 2, Figure 3 and Figure 4). Peaks of −45 to −50 mm/year LOS displacement rates are estimated along both the ascending and descending mode results within the city of Puebla, mainly in a large sector extending from San Juan Cuautlancingo to the north-west, towards the south through Santiago Momoxpan, San Pedro Cholula, San Baltazar Campeche, and Guadalupe Hidalgo. A notable area undergoing instability is also found in the north-east at the Parque Industrial Puebla 2000, where LOS displacement rates reach −40 mm/year (Figure 5a). The extraction of vertical and horizontal displacement rates through a combination of ascending and descending mode results (see Section 3.3), confirms that the vast majority of the motion occurs along the vertical direction and indicates land subsidence (Figure 5a). Extremely low to null rates are found along the E-W direction (Figure 5b), thus suggesting absence of significant horizontal strain at the margins of and within the subsidence bowls.

Comparison of groundwater levels in 2002 and 2010 for the whole aquifer of the valley of Puebla shows drops of 2–20 m, i.e., about 0.25 to 2.5 m/year, and around 2–10 m within the city (see water level decline information in Figure 5a). The greatest drop (15–20 m) is observed to the northwest of San Martín Texmelucan, while between 5 and 10 m are recorded around the town of Huejotzingo (near Puebla International airport) and in the region around the city of Puebla, and the towns of San Pedro Cholula, Santa María Coronango and Tlaxcalancingo [57]. During the preceding decades, between 1973 and 2002, a groundwater level decline within the city of over 20–40 m was measured, together with a cone of depression of 80 m in a 5 km-wide area located in the southwestern sector of the city of Puebla [58]. It is within the latter sectors that the highest rates of ground displacement are found through the InSAR analysis in 2003–2010. Vertical rates between −30 and −50 mm/year predominate and reach peaks of −53 mm/year (Figure 5a).

The 2003–2010 results agree extremely well (both in terms of spatial pattern and with respect to the observed displacement rates) with those published in the wide-area assessment based on ALOS PALSAR L-band data [12], where a maximum vertical displacement rate of −44 mm/year was estimated in 2007–2011. Around −30 mm/year LOS rates were also estimated for a more recent period (i.e., 2014–2015) through a wide-area SBAS-like analysis of Sentinel-1 imagery [16].

As also observed in [12], higher subsidence rates correlate spatially with industrial land use, especially in the areas of Parque Industrial Puebla 2000 and Industrial Resurrección, where over 20 groundwater wells operate [62]. This corroborates the conclusion that groundwater pumping for industrial activities in those sectors acts as the lead factor that determines ground deformation. Significant subsidence rates that are observed at alluvial deposits and intermediate tuffs in the larger sector encompassing the areas of San Juan Cuautlancingo, Santiago Momoxpan, San Pedro Cholula, San Baltazar Campeche, and Guadalupe Hidalgo, match well with the location of other agglomerates of operating groundwater wells devoted to industrial use and partly to public use. These are mainly located to the west of Río Atoyac (which flows from 40 km to the north of Tlaxcala towards the south up to the city of Puebla) and concentrate along the tributary Río El Zapatero.

4.1.2. Tlaxcala and Apizaco

Three smaller zones of ground movement are found in the cities of Tlaxcala and Apizaco (Figure 5c), located within the aquifer of Alto Atoyac, which extends a total of 2032 km², north of that of the valley of Puebla (Figure 3), and is fed by Río Atoyac.

To better understand the InSAR results, the geological setting and aquifer characteristics need to be accounted for. The oldest stratigraphic units in this area are Tertiary sedimentary rocks, in particular, clastic deposits of continental lacustrine origin (sandstone, limestone, and shale), onto which Miocene andesite and intermediate tuffs, and Pliocene acid tuffs are found [80]. Alluvial and fluvial deposits within channels and valleys host the unconfined aquifer, the main source to address water needs of the region. The deeper portion of the aquifer sits within a sequence of highly fractured extrusive igneous rocks, composed mainly of basalts, tuffs, and andesite. The barriers to horizontal and vertical flow are mainly the same volcanic rocks, where fracturing disappears at depths.

A total of 890 extraction sites were inventoried in 2011 within the Alto Atoyac aquifer, of which 734 wells (600 active and 134 inactive). The total extracted volume was over 150 million m³/year, 53% of which for public-urban use, 41% for irrigation and the remainder for services, livestock and domestic use [80]. Although the aquifer system records an overall water balance showing availability, for instance of more than 40 million m³/year in 2011 (see also Figure 3), comparison of groundwater levels measured in 1997 and 2011 revealed drops of 0.25 to 3 m, i.e., a maximum of 0.2 m/year, in some sectors of the aquifer. Peaks of 0.15–0.2 m/year groundwater level drop were found at Apizaco and Tlaxco, where pumping activity concentrates, and only 0.03 m/year at Tlaxcala.

In Tlaxcala, spatial patterns of InSAR-derived ground displacement seem to follow the river valley of Río Zahuapan at the meander, and vertical rates range from −5 to −17 mm/year (Figure 5c). This sector involves the Xicotencatl National Park designated area, including the old town of Tlaxcala and La Juventud recreation and sports park. Analysis of water levels [80] showed that these are very shallow between Tlaxcala and Zacatelco (at 5–25 m of depth from the ground surface), thus suggesting that compaction in Tlaxcala likely takes place in the upper layer, i.e., the unconfined aquifer consisting of fluvial deposits. In the areas showing the greatest subsidence rates, a combination of public and industrial use wells operate [62]. Similarly, in Apizaco the greatest subsidence rates (−6 to −8 mm/year) occur at alluvial deposits in the southern sector of the urban area, where a few groundwater wells devoted to industrial and public use operate.

4.1.3. The Valley of Tecamachalco

Vertical displacement rates of −10 to −15 mm/year, with peaks of −22 mm/year, are clearly depicted by the InSAR results within the aquifer of the valley of Tecamachalco (Figure 5d).

The aquifer system, as identified by CONAGUA, extends 3600 km² and encompasses three intermontane stepped valleys, namely the Esperanza valley with elevations around 2400 m a.s.l., Palmar de Bravo at 2200 m a.s.l. and Tepeaca-Tecamachalco-Tehuacán at 2000 m a.s.l., predominantly consisting of calcareous and sedimentary deposits (see also Section 2.4). The upper layer includes unconsolidated Quaternary deposits, while closer to the mountain ranges fractured rocks compose the aquifer [67]. Over 1400 water pumping wells in total are found within the three valleys (Figure 3) [62], where water extraction mainly for irrigation occurs since 1975 (see Section 2.4).

Subsidence in this region chiefly occurs in areas covered by Quaternary alluvial deposits, with observed spots of stability at volcanic cones and edifices, and no clear signs of deformation onto travertine and limestone. To the west of the town of Tecamachalco, the highest vertical rates occur within the villages of San Salvador Huixcolotla (−17 mm/year), San Juan Acozac, Los Reyes de Juárez and Actipan de Morelos (−22 mm/year), San Hipólito Xochiltenango and San Pablo Actipan, south of Tepeaca (−15 mm/year). Within the Palmar de Bravo valley, the highest rates are observed at San Miguel Xaltepec and Palmarito Tochapan (−19 mm/year), and partly at Palmar de Bravo (−12 mm/year), upper through the valley towards the southeast (Figure 5d). Other sectors showing displacement but at much lower rates (−5 mm/year) are southwest of Esperanza town, at the margins of a narrower sector of the Esperanza valley where a number of extraction wells operate.

No significant ground deformation rates are found to the south of the town of Tecamachalco, despite the moderate number of wells found, probably reflecting the less intense water pumping and slower declines in the static levels. A detailed hydrogeological study of water levels in this aquifer is however currently unavailable in the literature or in the monitoring reports by CONAGUA, so no comparison of the observed patterns of ground displacement with the measured water level declines is viable.

The valley of Tecamachalco falls within the wide-area analysis based on ALOS PALSAR imagery acquired in 2007–2011 [12], which showed peaks of around −20 mm/year vertical rates at San Pablo Actipan, and around −15 mm/year at San Hipólito Xochiltenango. These are consistent with the rates estimated based on the ENVISAT analysis presented in this paper. Less pronounced in the ALOS study were the displacement rates observed in San Juan Acozac, Los Reyes de Juárez and Actipan de Morelos (i.e., ~−10 mm/year), vs. the ENVISAT rates as high as −22 mm/year (Figure 5d). A more recent regional study carried out using Sentinel-1 data [16], also estimated around −20 mm/year LOS rates for the area of the Tecamachalco valley in 2014–2015, confirming the persistence of such rates also for a more recent period.

4.1.4. The Valley of Tulancingo

A pronounced displacement pattern is also identified in the south-western and western sectors of the municipality of Tulancingo de Bravo (over 150,000 inhabitants), with rates as high as −50 mm/year along the ENVISAT ascending LOS (Figure 4b). The pattern is characterized by rates increasing from east to west, namely from the western sector of the city of Tulancingo towards the margins of the urbanized area and transition into the valley (Figure 5e). Assuming that pure vertical motion occurs, the observed LOS estimates equal vertical rates of up to −55 mm/year.

Looking at the ENVISAT descending mode results, displacement rates of −20 to −25 mm/year along the LOS and at some locations exceeding −30 mm/year are also found (Figure 4a), corresponding to vertical rates of up to −33 mm/year. These mostly agree with estimations based on the ascending geometry at the same locations, though the ascending geometry provided many more coherent points towards the center of the valley, where rates are higher. It is to be noted that the sharp cut in the InSAR results in the descending mode at this location (Figure 4a) is due to the interruption of many of the 32 ENVISAT strips composing data stack #3, and the unavailability of a sufficient number of SAR images to perform a robust multi-temporal image processing for the area to the north. This, therefore, made unfeasible the combination of the ascending and descending mode results to derive vertical and horizontal components for this town. The estimation of vertical rates in Figure 5e is therefore based on the ascending mode data only (see Section 3.3).

Similarly to what discussed above for other subsiding locations, local geology at Tulancingo and the characteristics of the valley of Tulancingo aquifer provide further insights into the interpretation of the observed ground displacements. The municipality of Tulancingo is built within an intermontane alluvial plain, consisting mainly of 3–15 m thick Quaternary fluvial and lacustrine deposits (Figure 2) and bounded to the east and west by two systems of normal faults. The latter delimit the valley and act as barriers to horizontal water flow. The Rio Grande de Tulancingo river crosses the valley from south to north. Alluvial deposits sit onto 50–60 m of basalts and sandy-clayey pyroclastic deposits, which in turn overlay 100 m thick volcanic rocks of the Atotonilco Formation, and the Grupo Pachuca Formation [81].

The aquifer system of the valley of Tulancingo consists of two aquifers. The first is a shallow unconsolidated layer including the alluvium, basalts and pyroclastic material and conglomerates with a minimum thickness of around 200 m. The second deeper aquifer, mostly unconfined (or semi-confined only at places) and highly fractured, consists mainly of volcanic rocks, where groundwater extraction concentrates. Within the aquifer boundary (1054 km²), 313 groundwater extraction sites were inventoried in 2006 (of which 248 active), including 272 wells (see also Figure 3). The extracted total volume was 64 million m³/year for domestic, industrial, irrigation and other use, and the estimated deficit of groundwater resources was 7 million m³/year [81].

There is no record of pumping sites in the central and eastern sectors of the city [62], where the InSAR results show an apparent scenario of ground stability. Instead, over 15 of extraction wells for industrial and public use concentrate on the western, north-western and south-western sectors, at locations matching with the observed greatest rates of ground subsidence. Moving towards the valley to the west, the destination use of the wells is mainly for agriculture, and considerable subsidence rates persist, reaching up to −50 mm/year and even lower vertical rates.

Within this region of the aquifer system, the depth of static level from the ground surface was at 120 to 80 m in 2006, shallower in the valley of the Rio Grande de Tulancingo river. Comparison of static levels between 2000 and 2006 showed drops as high as 6 m within the valley [81], mainly in the northern sector, between the villages of Metepec, Colonia 28 de Mayo, and Alcholoya (Figure 5e), where the InSAR results identify an area with around −10 mm/year vertical displacement rates. Drops in the static level of up to 3–4 m between 2000 and 2006 were estimated along the fluvial and lacustrine deposits, where the highest rates of subsidence (−8 mm/year) are found and a dense concentration of groundwater extraction wells for agricultural use operate. Drops of 1–2 m in the static level were observed in the western sectors of Tulancingo.

In the whole area, the distribution of groundwater extraction wells, the pattern of static level drops and local geology suggest that ground deformation occurs as a consequence of water pumping and compaction of the aquifer system.

4.1.5. Eastern Mexico City Metropolitan Area

At the western boundary of the processed ascending data stacks (Figure 4b), another area showing significant ground displacement rates can be identified, with peaks of −44 mm/year vertical rates. This is over the municipalities of Ciudad Nezahualcóyotl, Ixtapaluca, and Valle de Chalco Solidaridad, located within the eastern Mexico City metropolitan area.

In this area, a number of InSAR based studies discussed the notable deformation patterns and rates observed, their relationship with intensive groundwater pumping and drawdown of the piezometric levels (e.g., up to 8–10 m drop recorded between 2006–2011 in this area, see Figure 5f) [82] and induced compaction of the thick layer of clay-rich deposits of former lakes Texcoco and Xaltocan (e.g., [8,9,12,16,83,84,85,86,87]).

Figure 3 shows that the aquifers encompassed by this subsiding area were all classified as overexploited in 2011, as they also were in the preceding decade. A thorough analysis of this area is beyond the scope of this paper, though the detail of the ground displacement rates as estimated via InSAR processing of ENVISAT imagery is provided in Figure 5f.

4.2. Environmental Baseline Mapping

Moderately dense datasets of InSAR results are found within the permit areas for geothermal exploration of Acoculco, Las Derrumbadas, and Cerro Pinto (see Figure 2 for location, and Figure 4 for the InSAR overview). The ground stability information retrieved contributes to appraising the environmental baseline setting for these sites, in view of their potential future exploitation.

4.2.1. Acoculco

The geothermal exploration area of Acoculco (Figure 6a–b) is covered by InSAR results from both ascending and descending passes, though during the image processing the chiefly forested and rural landscape impeded the identification of coherent radar targets in non-urban sectors, resulting in a relatively sparse coverage of results (Figure 6c–d). For this reason, differently from the approach used for the analysis of the valley of Puebla, the combination of ascending and descending results to derive vertical and horizontal rates in Acoculco was considered not suitable, and the two datasets of InSAR results were analyzed individually.

Both the viewing geometries provide an overall stable scenario, with null to extremely slow rates of ground displacement (i.e., ~−1.1 mm/year along the LOS on average within the exploration site) occurring in 2003–2010. The observed rates are thus in the range of the 1.1–1.2 mm/year uncertainty characterizing the InSAR datasets covering this area (see Section 3.2). No apparent regional pattern of ground deformation can be identified from the InSAR analysis, and the groundwater monitoring data from CONAGUA do not record any deficit over the last years for the aquifer of Tecolutla [60], i.e., the main aquifer across which the exploration area extends (see also Figure 3).

The results along the descending geometry (Figure 6c) identify only localized motion of a few mm/year occurring at locations at a significant distance from the areas with most apparent hydrothermal alteration, which are chiefly at Alcaparrosa, and at Los Azufres-El Potrero Colorado (Figure 7a). The latter area is where the first exploratory well (i.e., EAC-1, Figure 7b) was drilled in 1994 by CFE at up to a total depth of 2 km [50]. A second well was then drilled with 1.9 km total depth (i.e., EAC-2) in 2008 [88], and is located ~500 m east of EAC-1.

As observed during a field campaign carried out in May–June 2018, in the area of Los Azufres-El Potrero Colorado swamp there is a number of clear surface manifestations of geothermal activity, such as argillitic alteration at the surface, several active gas vents (carbon dioxide and hydrogen sulfide discharge), bubbling and acid-sulfate springs (Figure 7c–e). Springs and bubbling due to CO₂ emissions in the water seem to manifest along directions following the main geological structures recognizable in the site, namely long NW-SE trending faults, and a system of shorter SW-NE normal faults determining some ground depressions filled with lacustrine sediments deposited onto ignimbrite.

Despite the absence of a consistent density of InSAR coherent targets in the area of hydrothermal alteration, it can be reasonably concluded that there is no obvious evidence of any wide-scale deformation pattern ongoing in 2003–2010 within the exploration area. Within a wider perspective, general stability is also found over travertine outcrops 6 km south-west of the exploration wells, thus suggesting that no apparent local deformation pattern occurs at that location. As such, this ground stability scenario can be used among reference datasets for future assessment of possible volume changes in the reservoir—if and when exploitation of geothermal resources might be approved and take place at this site.

4.2.2. Las Derrumbadas and Cerro Pinto

The InSAR results for Las Derrumbadas complex and Cerro Pinto permit areas show a predominantly stable ground scenario, with some noticeable areas of land subsidence with around −5 to −10 mm/year vertical rates observed in 2003–2010 (Figure 8).

The most pronounced and well-defined feature is to the east of Totolcingo lake (or El Carmen), and extends from the town of Zacatepec and towards the south-west nearly up to San Salvador El Seco, for a total of around ~125 km², with peaks of −13 mm/year.

Localized subsidence patterns characterized by −6 mm/year rates can also be distinguished along the eastern and southern shore of Alchichica lagoon and around Quechulac lagoon. Evidence of targets moving at −4 mm/year is also found across the alluvial deposits and partly in the breccias and debris avalanches deposits north-east of Las Derrumbadas domes (see also Figure 2), and in the plain between Quechulac lagoon and the town of Guadalupe Victoria.

The InSAR results can be interpreted by considering the hydrogeological context of the region. Las Derrumbadas complex and Cerro Pinto are located, indeed, within the aquifer of Libres-Oriental (see also Figure 3), which extends 3500 km², from Los Humeros caldera to the north, towards the south up to the valley of Tecamachalco, and to the west towards La Malinche. The basin is endorheic, and surface water flow is limited (see Section 2.2). The aquifer mainly consists of extrusive igneous rocks, such as pyroclasts and lavas, except in the central and southern sector—where its upper section includes unconsolidated alluvial deposits—and in the mountain massifs, where it is constituted by fractured lava effusions. In the areas of Totolcingo and Tepeyahualco (or El Salado) lakes, the aquifer is semiconfined by clayey deposits of low permeability and thickness of various tens of meters, whilst in other zones, it can be locally semiconfined by fine-grained alluvial deposits.

Water level depths are around 2 m in the area of Totolcingo lake, while near Tepeyahualco lake water rises at the surface in the lagoons of Alchichica and Quechulac, and increases up to 40–100 m in the foothills of the mountains surrounding the valley [49]. Groundwater flows from the higher elevations towards the lacustrine area, where a pattern of groundwater flow typical of closed basins can be recognized. Water inputs come from the area of Huamantla (to the west of the aquifer) and water discharge occurs in the area of Libres and Oriental. Similarly, groundwater flow around Tepeyahualco does not converge into the lake but flows towards the area where pumping occurs to the west of the latter, and then towards the state of Veracruz.

Starting in the 1960s, a number of wells were drilled, extraction of groundwater (80% of which for agricultural use) increased significantly with up to 103 million m³ extracted in 2002, and water levels started to drop. As a consequence, some shallow lakes such as that of Tepeyahualco nearly dried out. Taking into account the distribution of groundwater wells across the region and the extraction activities, it can be concluded that the localized patterns of land subsidence observed in this area are mostly explained by compaction due to aquifer depletion. This interpretation also finds further supporting elements in the estimated water deficit of 13 million m³/year recorded by CONAGUA in 2008 [59]—and, therefore, referring to the years covered by the satellite data—vs. the availability of 17 million m³/year in the preceding years, as estimated in 2002 [49]. Future studies of surface deformation at Las Derrumbadas and Cerro Pinto will need to account for the dynamics of land subsidence observed across the Libres-Oriental aquifer, to accurately identify any potential displacement patterns due to different anthropogenic processes (e.g., geothermal resources exploitation and energy production, should a concession be granted).

4.3. Surface Deformation and Geothermal Exploitation at Los Humeros

InSAR results for Los Humeros geothermal field reveal the presence of land deformation within the caldera, whereas no large scale patterns of ground deformation are found across the wider region of interest. The highest subsidence rates detected within Los Humeros caldera are −6 to −8 mm/year with peaks of ~−11 mm/year and concentrate within Los Potreros collapse. Subsidence occurs mainly within a sector bounded by Los Potreros caldera rim to the east, Los Humeros and Maxtaloya faults to the west, and Los Humeros caldera rim to the south (Figure 9a,b). This sector is mainly covered by Cuicuiltic tuff (Holocene) and Maxtaloya trachyandesites (Upper Pleistocene), and partly olivine basaltic lavas west of Los Humeros fault.

Evidence of ground deformation from the InSAR analysis occurs at elevations of 2500–2800 m a.s.l., and indicates an apparent, very pronounced variation of ground deformation rates from −2 to −9 mm/year moving from the western rim of Los Humeros caldera and towards the east up to Los Humeros fault (Figure 9b,c). A more gentle, but of equal magnitude, variation of rates is also observed when moving from the eastern rim of Los Humeros caldera, towards Los Potreros caldera rim and further to the west up to Los Humeros fault. The subsidence pattern is thus pointing at the area with the greatest concentration of production wells (examples of which are shown in Figure 10).

To assess the reliability of our results, a comparison was conducted against InSAR results published in recent papers, either focused specifically on Los Humeros geothermal field [23] or including Los Humeros within the processed areas, albeit with a different study purpose [16].

Starting from the latter [16], a clear pattern of land subsidence can be recognized in the area of Los Humeros in the Sentinel-1 processing results (2014–2015) obtained with a modified SBAS implementation. The maximum LOS rate was of around −10 to −30 mm/year.

LOS displacement rates of up to −8 mm/year in 2003–2007 were also recently identified in [23] in the northern sector of the geothermal field, by processing a small set of 13 ENVISAT scenes with the Persistent Scatterers InSAR over a 30 by 35 km subset area. It is to be noted that, differently from the PS-InSAR approach used in [23], the implementation of an SBAS InSAR approach in this paper allowed the use of the full stack of 22 ENVISAT scenes (i.e., ID #1 in Table 1), without exclusion of any available image in the monitoring period. This was possible thanks to the multi-master approach typical of an SBAS implementation [5], and the creation of a redundant network of 56 small baseline interferograms (instead of 12 used in [23] with a single-master approach). This resulted not only in the successful identification and removal of atmospheric phase delays from the differential interferograms but also in the robust estimation of the full deformation histories and annual rates for the identified coherent targets.

The cause for the observed deformation was attributed by [23] to field operations (rather than to volcano-tectonic processes), owing to the maximum subsidence rates located around the production wells, the estimated source depth at 1650 m, and the short wavelength deformation pattern. The observed ground subsidence could, therefore, be explained by volume reduction in the andesite reservoir due to the extraction of geothermal fluids. Fast-forward modeling by [23] also allowed the estimation of a reservoir compaction volume ~50 times smaller than the net production volume, i.e., a much larger ratio than that observed at other geothermal fields such as Cerro Prieto. The relatively low rates of subsidence thus suggested that it is likely that pressure support of the reservoir takes place, thanks not only to injection at two wells (average yearly injection rates were ~0.3–0.5 million tons per year, i.e., an order of magnitude smaller than the average yearly production of ~5.2 to 5.7 million tons per year) but also to deep recharge.

Figure 11 compares the InSAR-derived deformation history for a set of coherent targets located in the area with the greatest number of production wells, and the geothermal energy production data for Los Humeros field as available in the literature [44]. The plot shows that in 2003–2007 the average steam production rates were between 500–600 t/h (i.e., a total of over ~5 million tons of steam produced each year), increasing from the rate of ~200 t/h in 2001. The increase in production observed starting in early 2002 was attributed by [22] to the occurrence of a 3.2 magnitude earthquake in January 2002, and the consequent expansion of the pores in the reservoir system, which facilitated fluid flow, thus increasing the production of steam.

From the deformation history in Figure 11, it can be observed that a +300 t/h boost in production rates occurred in early 2003 and a less strong boost of around +100 t/h was recorded at the end of 2004. While fluctuations of ~5 mm in the deformation history of the coherent targets can be associated with the uncertainty of the estimated time series, the overall trend observed across the 2003–2007 period clearly reveals a match between the two production boosts and periods of accelerated displacement. The more pronounced acceleration occurred throughout the year 2003, then deformation decelerated in the following years as production rates stabilized around 550 t/h.

The published literature reports about induced seismicity at low magnitude (generally <2.0) occurring at the local scale, with greatest density of epicenters in the area of the production wells, and to the west of the injection wells, as recorded by the 6 triaxial digital seismographs of the permanent telemetric seismic network installed by CFE in 1997 [90]. Water injection into the reservoir results in the fracturing of the rock skeleton, which in turn determines local seismicity. Over 230 earthquakes were recorded in 1997–2008 at 2000–2700 m depths, especially over the traces of La Cuesta, Los Humeros, Loma Blanca and Las Papas faults, where temperature is higher (~310 °C) [22]. On the other hand, the analysis of regional seismicity using data from the national seismic service (i.e., Servicio Sismológico Nacional, SSN, [32]) managed by the Institute of Geophysics of the National Autonomous University of Mexico (UNAM), confirms that no major events were recorded in the area of Los Humeros between 2003 and 2007, i.e., the period covered by ENVISAT data (see earthquake records in Figure 2). Figure 11 also shows that the subsidence process was persistent throughout the period 2003–2007, with no evidence of any other sudden displacement that could relate to time-limited, specific events of deformation. Therefore any hypothesis that the ground deformation observed within Los Humeros caldera in 2003–2007 could be due to a natural phenomenon such as an earthquake in the regional context of the TMVB, cannot be supported with these data.

Nevertheless, regional tectonics and the occurrence of high magnitude earthquakes over the last few years (i.e., after the period observed with ENVISAT data) contributed to damaging some of the housing and public buildings within the geothermal field, for instance in the village of Humeros. Two of the most recent and significant earthquakes that occurred in the geothermal field include the Mw 2.0 earthquake recorded at 1.6 km depth on 16th August 2015 and the already recalled Mw 4.2 earthquake recorded on 8th February 2016 [25,26,39]. Evidence of rockfalls triggered along Los Humeros fault is visible just outside the southern outskirts of Humeros. More recently, the 7.1 Mw Puebla earthquake struck on 19 September 2017.

Photographs taken during the field visits in May 2018 (Figure 12) show evidence of cracks opened in the façade and tower-bell of the church that was built in the 2000s. The diagonal orientation of the cracks on the masonry and brickwork and the presence of X-cracking on masonry piers suggest that these are not due to ground subsidence and settlement, but rather to earthquake shaking, and combination of vertical and horizontal loads.

Looking at the hydrogeological setting, the territory where Los Humeros caldera complex sits, spans across three aquifers, namely Libres-Oriental (3500 km²) within which lies the southern portion of the caldera, Tecolutla (7584 km²) within which lies the northern portion, and Perote-Zalayeta (1016 km²) covering its eastern flank. The first of these aquifers is the only one recording a groundwater deficit as of 2008 [59], but mainly attributed to intensive activities of pumping for agricultural use within the Serdán-Oriental basin (see Section 4.2.2). The monitoring report by CONAGUA for the Libres-Oriental aquifer, on the other hand, highlights that in the northern half of the aquifer, the presence of faults and fractures, the caldera of Los Humeros and other geological structures probably controls the groundwater flow, generating a deep hydrothermal system [49].

To the east of the caldera, within the Perote-Zalayeta aquifer (see also Figure 3), groundwater flows towards lower elevations at the center of the aquifer. In 2013, a total of 192 extraction sites were inventoried (of which 183 wells), pumping a total of nearly 33 million m³/year, mainly for agricultural use, and some as source of potable water and domestic use for the towns and villages of the valley [91]. In 1996–2013 localized drops of the water levels as high as 2–5 m were observed (corresponding to a rate of around 0.1–0.3 m/year), especially in the central and southern sectors of the valley, where the most of the groundwater pumping occurs [91]. In these areas, however, no alterations in the groundwater flow or formation of depletion cones were observed as a consequence of the pumping activity. On the other hand, the CONAGUA report for the Tecolutla aquifer [92] does not provide additional key information about significant groundwater level variations for the area north of Los Humeros caldera.

The above analysis of the ENVISAT InSAR results in the context of published literature and the wider geological context, finally suggest a spatial and temporal association of ground displacement rates and the location and distribution of geothermal energy infrastructure and production activities. Further research studies are ongoing to better constrain the relation with ground deformation processes, also in connection with regional tectonics and seismicity.

It is with this goal in mind that in 2016 the Universidad Michoacana de San Nicolás de Hidalgo (UMSNH) has started dedicated monitoring campaigns including GNSS surveying to estimate direction and magnitude of ground displacements occurring at the regional level. Together with multi-disciplinary studies carried out in the framework of the CONACYT-H2020 project GEMex [93], such investigation will allow improved characterization of structural and hydraulic properties of the reservoir, which could potentially act as essential information to better manage the geothermal resource in future.

5. Conclusions

This wide-area survey provides the first overview of ground displacement hotspots across the eastern Trans-Mexican Volcanic Belt (TMVB). An assorted representation of use-cases is provided with regard to land subsidence due to groundwater pumping for public, industrial and agricultural use, plus elements for the assessment of baselines and ground deformation scenarios in areas of geothermal exploration and energy production, respectively. The findings discussed in this paper corroborate earlier studies of ground stability carried out in some sectors of the study area (e.g., [12]). These results then expand the knowledge base by going into details with local-scale analysis of displacements at several locations including Puebla, Tulancingo, Tlaxcala, Apizaco and the Tecamachalco valley, plus at four sites of geothermal exploration and exploitation, only one of which was previously investigated with a multi-temporal InSAR method [23].

In this sense, the knowledge gathered thanks to this regional survey and hotspot identification can provide supporting information for land use management and planning in urban areas of the eastern TMVB, hazard and risk assessments for urban infrastructure (e.g., housing, transport infrastructure), as well as management of groundwater and energy resources.

In densely urbanized areas such as Puebla, Tulancingo, Tlaxcala, and Apizaco, direct impacts of ground subsidence on urban infrastructure are expected to be concentrated at the margins of the subsiding areas. Fissures and surface faulting, indeed, usually occur in areas of significant differential motion (e.g., [83,84,85]), where the horizontal gradients of subsidence are much more pronounced. On the other hand, the InSAR-derived findings for rural areas devoted to agriculture can enhance the understanding of aquifer characteristics (e.g., storage and compressibility, [94,95]) and deforming behavior in response to pumping activities. This would contribute to a more informed resource management in the context of water inputs and the overall aquifer balance.

Surface deformation analysis in the geothermal exploration sites of Acoculco, Las Derrumbadas, and Cerro Pinto will serve as a baseline reference against which any future studies of ground stability in the areas should be compared in order to identify any potential effects of energy exploitation. The case of Los Humeros—where operations started in the early 1990s, but no earlier SAR data were available to study the pre-exploitation scenario—in a way demonstrates the importance that such baseline knowledge could have in the future, should operations start at those three sites.

Analysis of InSAR-derived ground displacement rates and patterns at Los Humeros confirms that there is a need for further modeling and investigations to better understand deformation mechanisms and local structures at the site, assess the relationship with local seismicity, and surface and bedrock geology, and carry out a detailed correlation analysis of production at each well. In such context, it is clear the potential benefit arising from further satellite monitoring and multi-disciplinary investigations at the site, especially accounting for the ongoing site developments as part of Los Humeros III project [33,44]. Towards an informed geothermal energy exploitation based on space observations, future use of Sentinel-1 SAR imagery acquired at the site since 2014, together with the GNSS surveying initiated by UMSNH in 2016, may help to further develop the initial assessment and extend the monitoring to the 2010s and up to present.

Author Contributions

Conceptualization, F.C. and D.T.; methodology, F.C.; software, F.C.; validation, F.C., D.T., V.H.G.-M., J.A.M.-J. and A.J.-H.; formal analysis, F.C. and D.T.; writing—original draft preparation, F.C. and D.T.; writing—review and editing, V.H.G.-M., J.A.M.-J., O.H.G.-H. and A.J.-H.

Funding

V.H.G.-M., J.A.M.-J., O.H.G.-H. and A.J.-H. received funding to carry out fieldwork activities from the CONACYT project GEMex: Cooperation in Geothermal energy research Europe-Mexico for development of Enhanced Geothermal Systems and Superhot Geothermal Systems, grant agreement n.2015-04-268074.

Acknowledgments

ENVISAT Advanced SAR data were made available by ESA; NASA’s 1 arc-sec SRTM global DEM by USGS via the LP DAAC Global Data Explorer; Copernicus Sentinel-2 data from the Copernicus Open Access Hub. ENVISAT data processing was carried out with the Parallel-SBAS method, developed by CNR-IREA, within ESA’s Grid Processing On Demand (G-POD), and using ESA’s ESRIN CE 01 SL6 64bits.

Conflicts of Interest

The authors declare no conflict of interest. The funders of the GEMex project had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Figures and Table

Figure 1. Location of the area of interest within the eastern Trans-Mexican Volcanic Belt (TMVB) in Mexico (inset), and topography based on 30 m ground resolution SRTM digital surface model [30], with the indication of volcanoes, main mountain ranges, cities, and urban areas.

Figure 2. Simplified geology of the area of interest (1:250,000 scale, adapted from [31]), with indication of earthquakes recorded by the national seismic service in 2003–2010 [32] and the sites of geothermal exploration and exploitation: Las Derrumbadas, Cerro Pinto, Acoculco, and Los Humeros (boundaries of geothermal sites adapted from [27,28,29,33]).

Figure 3. Groundwater well distribution from the CONAGUA’s Public Registry of Water Rights (REPDA) database [62] and aquifer groundwater availability/deficit classification as of 2011 [63].

Figure 4. Results of the InSAR analysis for the area of interest, showing ground displacement velocity estimated along the (a) descending and (b) ascending mode Line-Of-Sight (LOS) of the ENVISAT Advanced SAR (ASAR) sensor, overlapped onto 10 m resolution Copernicus Sentinel-2 data 2019.

Figure 5. InSAR-derived displacement velocity in (a–b) the city of Puebla, (c) Tlaxcala and Apizaco, (d) the valley of Tecamachalco, (e) the town of Tulancingo, and (f) the eastern Mexico City metropolitan area, overlapped onto 10 m resolution Copernicus Sentinel-2 data 2019. (a) and (c–f) show vertical rates, while (b) shows E-W horizontal rates. Static level changes in (a–c,e,f) are those recorded in the aquifers of the valley of Puebla in 2002–2010, Alto Atoyac in 1997–2011, the valley of Tulancingo in 2000–2006, and Texcoco in 2006–2011, as documented by the National Water Commission (adapted from [57,80,81,82]). Aquifer names are indicated in blue.

Figure 6. (a) Key lineaments of the Acoculco caldera and identified areas of hydrothermal alteration, overlapped onto SRTM digital surface model [30], and (b) simplified geology of the Acoculco exploration site. (c–d) InSAR-derived displacement velocity in the Acoculco geothermal exploration site in 2003–2010, estimated along the (c) descending and (d) ascending Line-Of-Sight (LOS) of ENVISAT, overlapped onto 10 m resolution Copernicus Sentinel-2 data 2019. Data sources (adapted from): lineaments [50], geology and hydrothermal mapping [53,89], geothermal permit [27], exploration wells [50,88].

Figure 7. (a–e) Surface manifestations of hydrothermal activity in the area of Los Azufres-El Potrero Colorado pond within the Acoculco geothermal field. One of the two exploration wells drilled by CFE (i.e., EAC-1, see location in Figure 6a) is shown in (b). Photos from the May-June 2018 field survey carried out at Acoculco.

Figure 8. InSAR-derived displacement velocity observed in the geothermal exploration sites of Las Derrumbadas and Cerro Pinto in 2003–2007, in relation to groundwater well distribution from CONAGUA’s Public Registry of Water Rights (REPDA) database [62]. Data are overlapped onto 10 m resolution Copernicus Sentinel-2 data 2019. Boundary of geothermal permits adapted from [28,29].

Figure 9. (a) Simplified geological map of Los Humeros geothermal field (modified from [40,41]) with indication of the boundary of the concession area (boundary adapted from [33]), and (b) InSAR-derived vertical displacement velocity in 2003–2007, overlapped onto 10 m resolution Copernicus Sentinel-2 data 2019. (c) Vertical displacement velocity observed for coherent targets within a buffer distance of 600 m from cross-section A-A’. Elevation data derived from the 30 m ground resolution SRTM digital surface model [30].

Figure 10. Energy production infrastructure at Los Humeros geothermal field in the context of (a,b) land cover and use within the caldera complex, with (c) detail of a production well. Photos from the May–June 2018 field survey carried out at Los Humeros.

Figure 11. Comparison of InSAR-derived ground displacement in 2003–2007 for a set of coherent targets located in the area with the greatest number of production wells within the Los Humeros geothermal field, and the historic steam production data adapted from [44]. Periods of moderate increase in steam production are indicated in grey.

Figure 12. Cracks in the (a,c) bell-tower and (b) façade of the town church in Humeros, caused by earthquake shaking. Photos from the May–June 2018 field survey carried out at Los Humeros.

Table 1

Satellite ENVISAT ASAR IS2 data stacks, and results from SBAS InSAR analysis. Notation: ID, identifier of the data stack, Asc, ascending orbits, Desc, descending orbits, T, track. LOS velocity information reports minimum and maximum velocity observed within each frame.

ID	No. Scenes	Dates	Orbit	Track	Area [km2]	Identified Targets	Density [Targets/km2]	LOS Velocity [mm/year]
#1	22	2003–2007	Desc	T212	9,500	104,805	11	−11.7 – +8.1
#2	22	2003–2007	Desc	T212	11,300	239,428	21	−20.3 – +8.2
#3	32	2003–2010	Desc	T484	10,200	109,647	11	−30.6 – +7.4
#4	31	2003–2010	Desc	T484	11,500	118,720	10	−45.3 – +7.9
#5	29	2003–2010	Asc	T148	11,000	144,548	13	−51.0 – +10.2
#6	34	2003–2010	Asc	T148	11,100	125,757	11	−49.0 – +10.0