Introduction
The subsidence of sedimentary basins is a process that affects many cities worldwide (e.g., Amelung et al., 1999; Buckley et al., 2003; Chai et al., 2004; Dehghani et al., 2009; Ding et al., 2004; Le Mouélic et al., 2005; Lu and Liao, 2008; Marfai and King, 2007; Stramondo et al., 2008; Strozzi et al., 2001; Vilardo et al., 2009; Watson et al., 2002). In general, differential ground subsidence in urban areas is related to the withdrawal of groundwater to meet human, industrial and agricultural needs (Poland, 1984). The over-exploitation of the aquifers decreases the pore pressure and induces the compaction of the sediments as a result of the inelastic deformation of the aquitard (Amelung et al., 1999).
The firsts cases of differential subsidence in Mexico due to water extraction were observed in Mexico City in the early 1920s (Gayol, 1925). Subsidence increased rapidly from 1930 to 1960 due to the rapid growth of the population (Ortega-Guerrero et al., 1993). Carrillo (1948) was the first to relate directly the subsidence with the over-exploitation of the aquifers in Mexico City. In the 1980s, differential subsidence began to be observed and documented in several cities in central Mexico, such as: Aguascalientes, Morelia, Celaya, Irapuato, Salamanca, Querétaro, Toluca and San Luis Potosí (e.g., Aguirre-Díaz et al., 2000; Ávila-Olivera, 2008; Ávila-Olivera et al., 2010; Ávila-Olivera and Garduño-Monroy, 2006; Calderhead et al., 2011; Carreón-Freyre and Cerca, 2006; Castañeda et al., 1993, 1995; Chávez-Alegría, 2008; Cigna et al., 2011; 2012; Esquivel, 2009; Farina et al., 2007, 2008; Garduño-Monroy et al., 2001; López-Doncel et al., 2006; Mejía-Gómez and Sandoval-Minero, 2004; Pacheco-Martínez; 2007; Pacheco-Martínez and Arzate-Flores, 2007; Pacheco-Martínez et al., 2006; Schroeder-Aguirre, 2010; Trujillo-Candelaria, 1989, 2009; Vargas, 1999; Zermeño de León et al., 2004, 2005;).
More recently, interferometric techniques were used to study the rapid subsidence of Mexico City and its suburbs (e.g., Strozzi and Wegmüller, 1999; Strozzi et al., 2001; Carreón-Freyre et al., 2006; Cabral-Cano et al., 2008; López-Quiroz, 2008; López-Quiroz et al. 2009; Osmanolu et al., 2011). The common denominator shared by these cities is their location on lacustrine and fluvial-lacustrine basins filled with sediments that are heterogeneous in both composition and structure.
In the particular case of the city of Morelia, discussed in this paper, subsidence of the ground was first recognized in 1983, coinciding with the rapid increase of the urban population (Arreygue-Rocha et al., 2002; 2005). The differential subsidence was manifested on the surface as soil fractures damaging houses and buildings in the urban area (e.g., Garduño-Monroy et al., 2001; Ávila-Olivera and Garduño-Monroy, 2006; Ávila-Olivera, 2008). The Central Camionera and La Colina faults were the first physiographic features identified as resulting from the differential subsidence of the ground (Figure 1). The rates of subsidence reported from in situ measurements range from 4 to 6 cm/yr (Garduño-Monroy et al., 2001). To our knowledge, no routine and repeated measurement of ground deformation are done in the city.
Several faults affecting the city of Morelia have been mapped: La Paloma, La Colina, Central Camionera, Cuautla, Torremolinos, El Realito, La Soledad, Chapultepec and Ventura Puente (Ávila-Olivera, 2008; Ávila-Olivera and Garduño-Monroy, 2008; Cabral-Cano et al., 2010; Cigna et al., 2012) (Figure 1). The orientation of these faults is predominantly NE-SW and E-W. Several authors proposed that these faults are parallel to the mapped regional tectonic faults (Ávila-Olivera, 2008; Ávila -Olivera and Garduño-Monroy, 2006; Farina et al., 2007; 2008; Garduño-Monroy et al., 2001). Furthermore, Garduño-Monroy et al., (2001) suggested that La Colina and La Paloma faults are the surface expression of motion on basement faults of tectonic origin.
The early studies that measured subsidence in the city of Morelia using InSAR were based on interferograms produced using only images pairs (Table 1). More recently, Cabral-Cano et al. (2010) and Cigna et al. (2012) used the method of persistent scatterers (PSI) (Ferreti et al., 2000; 2001) to obtain time series of the subsidence at specific locations. Chaussard et al., (2014) used the Small Baselines (SBAS) methodology on ALOS images to estimate average subsidence in the city. Here, we extend the work of Cigna et al. (2012) using the method of SBAS (Berardino et al., 2002; Schmidt and Bürgmann, 2003; Usai, 2003; Cavalié et al., 2007; López-Quiroz et al., 2009) to obtain time series of the subsidence across several cross sections and specific locations. The purpose of applying this technique was to increase the signal coherence and to minimize its temporal and spatial decorrelation (Berardino et al., 2002; Yan et al., 2009). The use of this method allowed us to obtain reliable time series that reflect the evolution of the ground deformation as a function of time, induced by subsidence of the soil in Morelia.
Author | Method | Interferograms spanning time | Subsidence [mm/yr] |
---|---|---|---|
Farina et al. (2007) | DInSAR1 |
July 2003 - November 2004 December 2004- December 2005 |
10-35 |
Ávila-Olivera (2008) | DInSAR1 |
July 2003 - November 2004 December 2004- December 2005 |
10 |
Farina et al. (2008) | DInSAR1 |
July 2003 - November 2004 December 2004- December 2005 |
10-35 |
Ávila-Olivera et al. (2010) | DInSAR1 | July 2003 - November 2004 | 35 |
Cabral-Cano et al. (2010) | PSI2 | July 2003 - October 2009 | 50 |
Cigna et al. (2011) | PSI2 | July 2003 - October 2009 | 50 |
Cigna et al. (2012) | DInSAR1 | July 2003 - May 2009 | 70-80 |
PSI2 | 40-50 |
1 DInSAR - Differential Interferometry
2 PSI: Permanent Scatterer Interferometry
InSAR processing
Twenty-eight ENVISAT-ASAR images, acquired from May 2003 to September 2010, were processed using the software JPL/CalTech Repeat Orbit Interferometry Package (ROI_PAC) (Rosen et al., 2004). ROI_PAC was used to construct and to unwrap the interferograms. All scenes are descending-orbit images (track 69, frame 3213). The orbital fringes were removed using the satellite orbital corrections provided by the Department of Earth Observation and Space Systems (DEOS) of Delft University of Technology and by the European Space Agency (ESA). The topographic contribution was corrected using a 1” arc, digital elevation model (DEM) resampled from the 3” arc DEM of the Shuttle Radar Topography Mission (SRTM) (Farr and Kobrick, 2000). Multilooking by a factor of five was applied in the direction of azimuth. This process leads to a ground pixel resolution of ~20 x 20 m. In order to increase the signal to noise ratios, a nonlinear adaptive spatial filtering (Goldstein and Werner, 1998) was applied to all interferograms. As described below in more detail, the interferograms were unwrapped using the branch-cut algorithm (Goldstein et al., 1988).
In order to avoid loss of coherence in the interferograms, their baselines were limited to 1) spatial perpendicular baselines of less than 400 m; and 2) temporal baselines of less than 420 days. As a result, 65 interferograms were selected, constituting two separate groups (Figure 2). The first one spans from July 2003 to January 2007, and the second from December 2008 to September 2010. Unfortunately there are no InSAR acquisitions in the region from February 2007 to November 2008. Usai (2003) proposed to overcome large gaps in the interferometric data by creating synthetic interferograms that predict a particular deformation behavior. However, the period missing is very long to select a priori a reliable behavior. Therefore, the two observational periods are treated as independent groups.
Atmospheric phase delays due to the vertical stratification of the atmosphere (Hanssen, 2001) and to the orbital inaccuracies, for which we estimate the best fitting ‘twisted plane’ (Cavalié et al., 2007), were calculated simultaneously following López-Quiroz et al. (2009). The corrected interferograms were co-registered prior to the inversion process. The phase delay time series were obtained using a least-squares inversion (e.g.,Cavalié et al., 2007). The mean square errors (RMS) of all interferograms and pixels were calculated in order to detect unwrapping errors and inconsistencies in the interferometric dataset (Appendix A).
Finally, average subsidence rates were calculated as a function of time for all pixels meeting the following criteria: 1) To have the complete number of acquisitions and; 2) To show an RMS error of less than 0.5 rad. All analyses were performed in radar geometry and later geo-referenced to the DEM.
Description of the inversion process
As suggested by Schmidt and Bürgmann (2003), the inversion process was performed without smoothing in order to avoid aliasing effects. The names of the interferograms are identified by the year, month and day of the images used to calculate them. Interferograms 20040313-20040417, 20061223-20070127, 20100220-20100327, 20100327-20100501, 20100327-20100605, and 20100501-20100605 were discarded because they did not meet the quality criteria described before. Furthermore, all interferograms sharing the common image 20040731 (20040207-20040731, 20040417-20040731, 20040731-20050122, 20040731-2050820) were not included in the analysis due to the strong atmospheric effects observed in that particular scene, which greatly increased the RMS errors.
Fifty-five interferograms were used in the analysis (Tables 2 and 3). The RMS values of the interferograms used to evaluate the time series are described in equations 2 and 5 in Appendix A and graphically shown on Figure 3.
ID1 | Interferogram1 | Bt2[days] | BD3[m] | ID1 | Interferogram1 | Bt2[days] | BD3[m] |
---|---|---|---|---|---|---|---|
1 | 20030712-20030816 | 35 | -83 | 16 | 20041218-20051203 | 350 | -63 |
2* | 20030712-20041113 | 490 | -28 | 17 | 20041218-20060211 | 420 | -141 |
3* | 20030816-20041113 | 455 | 55 | 18 | 20050122-20050611 | 140 | -335 |
4* | 20040207-20040313 | 35 | 457 | 19 | 20050122-20050820 | 210 | -16 |
5 | 20040207-20040417 | 70 | 335 | 20 | 20050607-20050716 | 70 | -23 |
6 | 20040207-20050122 | 350 | -167 | 21 | 20050611-20050820 | 70 | 319 |
7 | 20040313-20050507 | 420 | 188 | 22 | 20050611-20051203 | 175 | -182 |
81 | 20040313-20050716 | 490 | 165 | 23 | 20050611-20060211 | 245 | -260 |
9 | 20040417-20050507 | 385 | 310 | 24 | 20050611-20060527 | 350 | -335 |
101 | 20040417-20050716 | 455 | 287 | 25 | 20051203-20060211 | 70 | -78 |
11 | 20041113-20041218 | 35 | 355 | 26 | 20051203-20060527 | 175 | -153 |
12 | 20041113-20051203 | 385 | 292 | 27 | 20051203-20070127 | 420 | 166 |
131 | 20041113-20060211 | 455 | 214 | 28 | 20060211-20060527 | 105 | -75 |
14 | 20041218-20050122 | 35 | 454 | 29 | 20060211-20070127 | 350 | 244 |
15 | 20041218-20050611 | 175 | 119 | 30 | 20060527-20070127 | 245 | 319 |
1 Interferogram number. The interferogram was built using the satellite image pairs expressed as year, month and day. The asterisk indicates the interferogram did not meet the prescribed quality criteria.
2 Bt is the temporal baseline between subsequent radar images
3 BD is the spatial perpendicular baseline between subsequent radar images
ID1 | Interferogram1 | Bt2[days] | BD3[m] | ID1 | Interferogram1 | Bt2[days] | BD3[m] |
---|---|---|---|---|---|---|---|
1 | 20081227-20091003 | 280 | -1 | 14 | 20100327-20100710 | 105 | -150 |
2 | 20091003-20100220 | 140 | -90 | 15 | 20100327-20100814 | 140 | -266 |
3 | 20091003-20100327 | 175 | 224 | 16 | 20100327-20100918 | 175 | -79 |
4 | 20091003-20100501 | 210 | 207 | 17 | 20100501-20100710 | 70 | -133 |
5 | 20091003-20100605 | 245 | 178 | 18 | 20100501-20100814 | 105 | -249 |
6 | 20091003-20100710 | 280 | 74 | 19 | 20100501-20100918 | 140 | -62 |
7 | 20091003-20100814 | 315 | -42 | 20 | 20100605-20100710 | 35 | -104 |
8 | 20091003-20100918 | 350 | 145 | 21 | 20100605-20100814 | 70 | -220 |
9 | 20100220-20100501 | 70 | 297 | 22 | 20100605-20100918 | 105 | -33 |
10 | 20100220-20100605 | 105 | 268 | 23 | 20100710-20100814 | 35 | -116 |
11 | 20100220-20100710 | 140 | 164 | 24 | 20100710-20100918 | 70 | 71 |
12 | 20100220-20100814 | 175 | 48 | 25 | 20100814-20100918 | 35 | 187 |
13 | 20100220-20100918 | 210 | 235 |
1 Interferogram number. The interferogram was built using the satellite image pairs expressed as year, month and day.
2 Bt is the temporal baseline between subsequent radar images
3 BD is the spatial perpendicular baseline between subsequent radar images
Discussion of results
Regional ground subsidence in Morelia
The results of the analysis of the two independent time series show that, on a regional scale, most of the urban sprawl of Morelia is not affected by significant ground subsidence (Figures 4a and 4b). During the seven years studied, the average surface deformation was ±0.5 cm/yr. These values are within the error of the methodology. The areas showing subsidence occur on the footwall of some of the mapped faults. The subsidence, however, is not continuous along the whole length of the faults. Subsidence is concentrated in the western part of the city and on specific segments of La Paloma, Central Camionera and La Colina faults. The areas showing ground deformation are consistently the same in the two time series spanning this ten-year period, indicating that the subsidence is consistent in time and areal extent (Figures 4a and 4b). However, our results show a small general increase in the subsidence rate of the affected areas from December 2008 to September 2010, relative to the one observed from July 2003 to January 2007 (Figures 4a and 4b). Cigna et al. (2012) observed also this increase in subsidence in the latter years.
Although La Paloma fault represents one of the starker lithological contacts in the city, between the Miocene ignimbrite sequence of the Sierra Mil Cumbres and the basin fill formed by Quaternary sedimentary deposits and cemented tuffs (Figure 1), subsidence is observed mostly in the western end of the fault (Figures 4a and 4b). In this area of La Paloma fault, the zones that were mostly affected by a rapid rate of subsidence are residential areas in the southwestern part of the city: INFONAVIT Villa Universidad, Valle Quieto, Arboledas, and Residencial del Sur.
In the central part of the city of Morelia, the Chapultepec and Cuautla faults show no evidence of ground deformation during the seven years of observation (Figures 1, 4a and 4b). Both time series indicate that no discernable subsidence takes place near these faults during the observation period (Figures 4a and 4b). Farther to the north, the Central Camionera fault reflects rapid subsidence in the central part, gradually tapering off towards the west. The eastern Central Camionera fault shows no evidence of subsidence. The relatively rapid subsidence observed in the central segment of the Central Camionera fault affects the neighborhoods of El Porvenir, Ampliacion Porvenir and Las Flores.
La Colina fault, to the northwest of the city, reflects the geological boundary between the Pliocene basalts and andesites and the sediment fill of Quaternary age (Figure 1). The footwall of La Colina fault reflects subsidence along its whole length (Figures 4a and 4b). Comparing the two periods analyzed, the largest increase in the rate of subsidence took place along the eastern end of La Colina fault from December 2008 to September 2010 (Figure 4b). Subsidence on La Colina fault affects the housing units of INFONAVIT Adolfo López Mateos, Las Águilas, Agua Clara, and the neighborhoods La Quemada, Irrigacion and Vicente Guerrero.
Cigna et al. (2012) did not report a correlation between the ground subsidence in the city and the rate of water extraction and the decrease of the static level of the water wells. In the same manner, a comparison of the ground subsidence with the thickness of the Quaternary deposits showed no relation to the areas of high subsidence rates. Our results confirm the observations made by Cigna et al. (2012), suggesting that the subsidence is not directly and exclusively controlled by the thickness of the sediments or by the water extraction rates. The variables controlling the subsidence of the ground appear to be the complex structural distribution of the volcanic deposits forming the Morelia basin and the preexisting basement topography.
Subsidence time series and cross sections in the areas of higher subsidence
As mentioned above, we use the method of Small Baselines (SBAS) (Berardino et al., 2002; Schmidt and Bürgmann, 2003; Usai, 2003; Cavalié et al., 2007; López-Quiroz et al., 2009) in order to extend the work of Cigna et al. (2012). The use of this method allowed us to obtain reliable time series that reflect the evolution of the deformation induced by the subsidence in Morelia. Based on this approach, time series of the ground deformation were estimated across several cross sections of the faults. Also, rates of subsidence as a function of time were obtained for several locations in the city (Figures 4a, 4b and 5).
In order to construct the time series and the cross sections, only pixels that met the criteria of having a RMS < 0.5 rad were selected (see Equation 3 in Appendix A). In addition, the selected pixels had to be present in all interferograms. These criteria guarantee a reliable depiction of the behavior and velocity of the ground deformation (e.g.,López-Quiroz, 2008). All the cross sections and time series calculated are shown in the Electronic Supplement to this paper. In the discussion, we select some relevant examples to illustrate and discuss the observations.
Deformation in the vicinity of La Paloma fault
Three cross sections were calculated across the La Paloma fault (Sections CC’, DD’ and EE’ on Figures 4a and 4b). Also, subsidence time histories were estimated for representative locations near the fault (points D, E, F, and G on Figure 5). Points D and E show that the hanging wall of La Paloma fault, on the Mil Cumbres Sierra, is stable in the seven-year period observed by the data. In these two points the rate of deformation is less than 1 mm/yr (see Electronic Supplement). Immediately to the north of point E, point G shows a rate of subsidence of approximately 0.7 cm/yr (Figure 6). This rate of subsidence is maintained consistently throughout the seven years covered in this study. Point F, in the western end of La Paloma fault shows a rate of subsidence of 1.8 cm/yr in the period from July 2003 to January 2007 (Figure 6); the subsidence rate later increased to 2.6 cm/yr, from December 2008 to September 2010 (Figure 6). In all cases, the subsidence of the ground behaves quasi linearly as a function of time, displaying no significant seasonal variations.
Cross sections across La Paloma fault show very different subsidence histories. Sections CC’ and DD’, located in the central part of La Paloma fault (Figure 4a), show an average rate of subsidence of approximately 1 cm/yr; a value similar to the subsidence rate observed in location G (see Figure 7 and Electronic Supplement). In the case of cross section EE’, located in the western end of La Paloma fault, no discernible subsidence is present (Figure 7). The results of these three cross sections drawn across La Paloma fault depict the variability of ground subsidence in Morelia in different segments of the same fault.
Subsidence near Central Camionera fault
Subsidence time histories were estimated for several locations and for six cross sections in the vicinity of the Central Camionera fault. As in the case of La Paloma fault, the cross sections indicate great variability of the ground subsidence rates along the strike of the fault (see the Electronic Supplement). Point H lies on the hanging wall of the Central Camionera fault and shows no appreciable subsidence during the time covered (Figure 8). The area called Prados Verdes (point I on Figure 5) is on the footwall of the fault. This area shows the largest ground subsidence rate observed in the city (Figure 8). During the first period of observation, point I shows a relatively low rate of subsidence (Figure 8). During this time, the water well located in the vicinity of this area was not operational due to the installation of new casing. The bore diameter was increased from 96 to 102 m (Cigna et al., 2012). The drastic increase of subsidence, going from 1.2 to 4.4 cm/yr during the second observational period (Figure 8), is presumably related to the increase of the water extraction rate and a further reduction of the static level of the well (Cigna et al., 2012). To the west of the Prados Verdes area, point K, also on the footwall of the Central Camionera fault, shows a stable and linear rate of subsidence of approximately 1 cm/yr.
Six cross sections were calculated across the Central Camionera fault. In the central part of the fault, on cross section II’, there is a clear change in the subsidence rate on both sides of the fault. Along this cross section, the variation in the rates of subsidence, observed also in point I, is clearly evident (Figure 9). This rapid change in the rate of subsidence is probably due also to the increase in production of Prados Verdes well. In the eastern part of the fault (Section HH’) there is no discernable ground deformation during the seven years of the study (Figure 9). A similar result is observed also in cross section LL’ (see the Electronic Supplement). The gaps observed in cross sections HH´and II´ are due to points on the ground where coherence was lost and where the RMS values were above the prescribed levels of quality.
Subsidence in the footwall of La Colina fault
In the area near La Colina fault, several locations were selected to calculate the variation of subsidence as a function of time. Points N and P, on the hanging of the fault, show a very similar subsidence behavior. In the first three years of observation, the subsidence in both cases ranged from 1.7 cm/yr in Point N, to 1.8 cm/yr in Point P (Figure 10). From 2009 to 2010, the subsidence rate increased to almost 3 cm/yr in point P and to 2 cm/yr in Point N, again showing an increase in the subsidence rate of the ground in the second period studied. In all cases, however, the behavior is linear, showing no temporal or seasonal variations. Point O is also in the footwall of La Colina Fault. However, this location is farther away from the fault trace than points N and P. The average rate of subsidence is slower than in the previous two locations, showing an average subsidence of between 0.6 and 0.8 cm/yr (Figure 10).
Two cross sections were calculated crossing La Colina fault (cross sections NN’ and OO´ (Figure 11). Cross sections NN´ and OO´ show subsidence of the ground to the north of the surface expression of La Colina fault. In both cases, a broad area subsides at a maximum rate of approximately 2 to 2.5 cm/yr. The deformation of the ground stops to the north of the city, where the Pliocene basalt and andesite deposits delineate the northern boundary of the Morelia basin (Figure 11). The deformation rate on these two cross sections remains constant during the two periods studied. In both cases, some coherence is lost in the central part of the cross sections.
Summary and conclusions
A new analysis of ground subsidence in the city of Morelia, in central Mexico, is presented using the small baselines methodology to estimate interferograms. The purpose of this work is to extend the work done by Cigna et al. (2011); (2012) in mapping the ground subsidence in the city. The identification of the areas of the city where subsidence occurs and the estimation of subsidence rates in various locations is important to mitigate and predict the locations of the city that will be more prone to damage and ground faulting in the future.
From a regional point of view, most of the city of Morelia shows no significant subsidence. In fact, there is no deformation observed on many of the mapped faults, such as the Chapultepec, Cuautla and Torremolinos. The subsidence evidenced by the InSAR observations is concentrated in the western half of the city and along three major faults: La Paloma, Central Camionera and La Colina. However, the ground deformation is not continuous and homogeneous along these faults. Subsidence takes place only on specific locations along the footwall of these three faults. La Colina fault is the one that exhibits a more continuous deformation of the ground along strike.
The largest subsidence rate of 4 cm/yr is observed to the north of the central part of the Central Camionera fault. This region showed a rapid increase in subsidence from 2009 to 2010, probably due to the increased flow of Prados Verdes well. In general, however, there is no direct correlation between the rates of water extraction and the areas of larger subsidence. In the same manner, as observed by Cigna et al. (2012), there is no clear correlation between the thickness of the Quaternary sediments in different parts of the Morelia basin and the rate of subsidence.
Garduño Monroy et al. (2001) suggested that the ground deformation in the city of Morelia is the surface expression of motion on subsurface basement faults. These authors suggest that these buried faults are continuously deforming by aseismic creep. These authors attribute the presence of surface deformation of the ground as a reflection of the motion on basement faults beneath the sediments, caused by creep at depth. However, the rates of subsidence observed on the ground are of several cm/yr. The rates of subsidence observed in Morelia are almost an order of magnitude larger than creep deformation measured on active faults. For example, fault creep observed on San Andreas is in the order of mm/yr (e.g.,Schulz et al., 1982; Lienkaemper et al., 2001; Lyons and Sandwell, 2003). Thus slow aseismic creep motion on subsurface faults would not explain the large subsidence rates observed in the Morelia basin.
The ground subsidence observed in the city of Morelia is highly variable in areal extent and concentrated in specific areas. The lack of direct correlation between subsidence of the ground and water extraction and the thickness of the sediments suggests that the deformation of the ground is due to the preexisting topography of the basement that controlled the complex depositional history of volcanic and lacustrine deposits in the Morelia basin. This poses a challenge to manage the exploitation of the aquifer in order to minimize and control the ground deformation.