1. Introduction
Ground subsidence can be caused by several geological factors, climatic processes and anthropogenic sources, or by the mixing of the above phenomena. The subsidence is frequently linked to intense faulting and opening of fissures in urban areas, generating a significant geologic hazard that needs to be accurately assessed and monitored [
1,
2,
3]. Both standard InSAR (Interferometric Synthetic Aperture Radar), based on the analysis of a pair of SAR images, and Advanced InSAR (A-InSAR), based on the analysis of series of SAR images, have been applied to assess natural hazards such as processes of slow and/or fast surface movements [
4]. InSAR measurements are particularly capable of mapping ground deformation with a very high spatial resolution over a large area, with high precision and a moderate-to-zero cost [
5,
6,
7,
8]. Among others, SAR processing techniques have been successfully used to study vertical surface movements caused by anthropogenic ground subsidence, like mining [
9,
10] and the withdrawal of subsoil fluids in combination with loading [
11,
12,
13,
14,
15,
16]. Urban areas located within confined basins in active tectonic and volcanic environments are especially subject to the dangerous effects of the ground subsidence, faulting and fissuring, due to the combination of neotectonics, seismicity, thick unconsolidated sedimentary deposits and anthropogenic activities [
17,
18,
19].
One of the intense surface deformations linked with ground subsidence in urban areas occurred in the Mexican area of Ciudad Guzmán (CG) (Jalisco state) [
20]. On 21 September 2012, in the center of CG several fissures opened, causing deformations of roads and serious damages to houses and facilities [
20,
21]. This event was very fast, with no precursors and not related to significant seismic activity. These characteristics, the coexistence of active faults, thick recent sediments and human activities make the CG 2012 fissuring event an interesting case study for the understanding of the mechanisms and evolution of ground deformation in urban areas.
This work presents the results of a detailed field study and InSAR analysis of the CG surface deformations, focusing on the 2003–2012 period. Through the interferometric processing of SAR images of CG and surrounding area with standard and advanced methods, our study attempts to quantify and characterize the behavior of the subsidence deformation processes. Remote sensing data, accurate field mapping of brittle deformations, and numerical modeling of the subsidence process are used in order to propose a genetic model of the observed deformations. The final aim of our paper is to contribute to a better knowledge of ground subsidence in urban areas located within active tectonic and volcanic environments, showing that the integration of InSAR remote sensing, structural fieldwork and numerical modeling is an effective methodological approach for the study of these hazardous geological processes.
2. Geological Setting and Ground Subsidence of Ciudad Guzmán
The Trans-Mexican volcanic Belt (TMVB) is a 1200-km-long active continental volcanic arc originated by the subduction of the Cocos and Rivera plates along the Middle American Trench [
22] (
Figure 1). Some of the main Mexican cities in the TMVB, such as Mexico City, Querétaro, Morelia, Toluca, Guadalayara, Puebla and CG experience surface subsidence and ground fissuring [
17,
18,
23]. These cities are located in lacustrine basins bounded by fault scarps and volcanic reliefs. The sedimentary fill of the basins is highly heterogeneous in composition and texture, ranging from fine lacustrine sediments to coarse alluvial, colluvial and volcaniclastic deposits [
24]. The neotectonic activity, seismic shacking, compaction of unconsolidated recent sediments and human activities, like water withdrawal and ground loading within these basins, have been considered as the principal causes of the occurrence of ground subsidence and subsequent soil fissuring in the urban areas of the TMVB [
17,
18,
23]. CG, with a population of about 100,000 inhabitants, is located in the Mexican state of Jalisco, at 1500 m a.s.l., inside the Colima Rift, which is the southern branch of the Colima-Tepic-Chapala triple junction in the western sector of the TMVB [
17,
25,
26,
27,
28,
29,
30,
31] (
Figure 1).
Figure 1.
Main tectonic and volcanotectonic structures in the Northern Colima Graben (NCG). The red box in the upper left panel locates the Colima and CG area, the blue box in the main panel locates
Figure 2 while the white box locate
Figure 3. CG: Ciudad Guzmán, NC: Nevado de Colima volcano, FC: Fuego de Colima volcano, LZ: Laguna Zapoltlan (modified from [
32]).
Figure 1.
Main tectonic and volcanotectonic structures in the Northern Colima Graben (NCG). The red box in the upper left panel locates the Colima and CG area, the blue box in the main panel locates
Figure 2 while the white box locate
Figure 3. CG: Ciudad Guzmán, NC: Nevado de Colima volcano, FC: Fuego de Colima volcano, LZ: Laguna Zapoltlan (modified from [
32]).
The Colima Rift consists of three structural segments, the Northern Colima Graben (NCG), the Central Colima Graben, and the Southern Colima Graben. The NCG, where CG is located, is flanked by reliefs consisting of Late Miocene-Pleistocene volcanic deposits, Jurassic-Eocene sedimentary and intrusive rocks. The depression is floored by Pliocene-Holocene lacustrine sediments, alluvium, colluvium and volcaniclastic deposits of the nearby Colima Volcanic Complex (CVC) [
26,
27] (
Figure 2). The NGC is 20 km wide and 60 km long, and flanked by sharp and parallel NNE-SSW-trending active faults. Bounding faults in the NCG dip 70° toward the graben axis, the relief of their fault scarps is up to 1–2 km, and the mean displacement rate is up to 1–3 mm/year [
26,
27,
31]. The kinematics of these faults is normal with a minor right-lateral strike-slip component of motion, consistent with a minimum principal stress oriented from E-W to NW-SE [
26,
27,
28,
31].
Figure 2.
(
a) Geological map of Ciudad Guzmán surroundings with the location of the 2012 fissures (purple dots, see
Figure 3, here located with light blue dashed box). (Ba) Basanites—Quaternary deposits; (VS) Volcanic Sediments; (BT) Brown Tuffs—Plio-Quaternary deposits; (At) Andesitic Tuffs—Pliocene; (Rb) Red beds—Late Cretaceous; (Li) Limestone; and (Da) Dacites—Early Cretaceous. (
b) Geological cross-section of the A-A’ trace in
Figure 2a, with the location of the 2012 fissures (modified from [
32]). (
c) Geometry of the finite element model, taken from the black rectangle in
Figure 2b. The model has been made with MSC Marc 2013 software [
33].
Figure 2.
(
a) Geological map of Ciudad Guzmán surroundings with the location of the 2012 fissures (purple dots, see
Figure 3, here located with light blue dashed box). (Ba) Basanites—Quaternary deposits; (VS) Volcanic Sediments; (BT) Brown Tuffs—Plio-Quaternary deposits; (At) Andesitic Tuffs—Pliocene; (Rb) Red beds—Late Cretaceous; (Li) Limestone; and (Da) Dacites—Early Cretaceous. (
b) Geological cross-section of the A-A’ trace in
Figure 2a, with the location of the 2012 fissures (modified from [
32]). (
c) Geometry of the finite element model, taken from the black rectangle in
Figure 2b. The model has been made with MSC Marc 2013 software [
33].
The urban area of CG is bordered to the east by a Cretaceous continental sedimentary sequence exposed along the main fault scarps (red beds), and lies on the graben fill sequence composed of volcaniclastic deposits with intercalation of alluvial and lacustrine sediments. The thickness of the mainly unconsolidated graben fill sequence in the CG area ranges from 300 to 1200 m [
32] (
Figure 2).
This area is exposed to hazardous natural events such as landslides, volcanic eruptions and earthquakes. The seismic hazard is linked to the main tectonic structures of the volcanic arc and subduction zone, capable of strong earthquakes, and to minor local active faults generating moderate earthquakes [
34,
35,
36]. Earthquakes have hit the NCG several times, in 1911, 1931, 1932, 1941 and 1973 [
35]. The latest seismic events that affected CG occurred in 1985 (M8.0), 1995 (M7.4) and 2003 (M7.5) [
36]. Some of these earthquakes (e.g., 1985) generated ground fissuring in the urban area of the town [
37]. In addition, ground fissures opened in CG without any significant seismic shake (e.g
. 1993). Indeed, on 21 September 2012, the urban area of CG experienced intense ground fissuring, causing deformations of roads and serious damages to some houses. The 2012 fissures opened in the same location reported for a number of NE-SW striking superficial cracks opened in 1985 and 1993. These fissures have been associated to the faults cropping out along the border of the graben, which are buried by recent sediments under the urbanized area [
21]. The withdrawal of water from the ground has been considered the main driving factor of the surface subsidence and the opening of cracks in the city [
20,
32,
38]. The upper sediment sequence under the urban area consists of heterogeneous sediments with relatively low cohesion, that can be affected by water table dynamics and piping [
32].
4. InSAR Analysis: Methods and Results
We analyzed ground deformation in the CG area by applying the classical and multi-temporal InSAR techniques to the 2003–2010 ENVISAT and 2012 RADARSAT-2 data (
Table 2). In order to investigate the ground deformation during 2003–2010 period (prior to the rupture), we used 40 ascending and 41 descending ENVISAT images acquired by the ASAR (Advanced Synthetic Aperture Radar) sensor during March 2003–October 2009 and December 2003–August 2010 periods, respectively. These two sets of images were processed with the Multi-Baseline method implemented in GAMMA software [
39]. This method computes deformation time series and residual topographic heights using the Singular Value Decomposition (SVD) Least-Squares inversion technique [
40,
41]. All interferograms were computed with a multi-look factor of 1 × 5 in range and azimuth directions, respectively, leading to a 20 × 20 m pixel on the ground. After the computation of the interferograms, only coherent points,
i.e., pixels characterized by signal-to-noise-ratio of the interferometric phase higher than 1.3 were selected. This resulted in a stack of 42 ascending and 58 descending point-wise interferograms used for the SVD inversion. The velocity fields computed from ascending and descending data are shown in
Figure 5. Both velocity maps show a Line Of Sight (LOS) distance increase (surface movement away from the satellite,
i.e., subsidence) in the area NW of CG. This area is confined by the location of the 2012 surface fissures (
Figure 2,
Figure 3 and
Figure 5). The other side of the ruptures alignment (Eastern side of CG) is stable. It is worth noting that a step gradient in the surface velocity is located exactly where the surface cracks appear. The deformation rates measured by ENVISAT InSAR dataset in the subsiding area reach values of up to 25 mm/year for ascending data and of about 17 mm/year for descending data, and both deformation fields show similar spatial pattern.
Table 2.
The ENVISAT 2003–2010 dataset (left); and list of interferograms created using four 2012 RADARSAT-2 images (6 March; 26 June; 6 September; 11 December) (right).
Table 2.
The ENVISAT 2003–2010 dataset (left); and list of interferograms created using four 2012 RADARSAT-2 images (6 March; 26 June; 6 September; 11 December) (right).
ENVISAT | RADARSAT-2 |
---|
Ascending Orbit—Figure 5a | Descending Orbit—Figure 5b | Interferogram | Figure |
---|
20030117 | 20040416 | 20050506 | 20060804 | 20030124 | 20050826 | 20061020 | 20061020 | 20120306–20120626 | 8A |
20030328 | 20040521 | 20050610 | 20061117 | 20031205 | 20050930 | 20061124 | 20061124 | 20120626–20120906 | 8B |
20030502 | 20040730 | 20050715 | 20061222 | 20040319 | 20051104 | 20061229 | 20061229 | 20120906–20120306 | 9 |
20030815 | 20040903 | 20050819 | 20070126 | 20040423 | 20051209 | 20070309 | 20070309 | | |
20030919 | 20041008 | 20050923 | 20071207 | 20040528 | 20060113 | 20070413 | 20070413 | | |
20031024 | 20041112 | 20051028 | 20080912 | 20040806 | 20060217 | 20070518 | 20070518 | | |
20031128 | 20041217 | 20051202 | 20081121 | 20041119 | 20060428 | 20071005 | 20071005 | | |
20040102 | 20050121 | 20060106 | 20090306 | 20050513 | 20060602 | 20071109 | 20071109 | | |
20040206 | 20050225 | 20060317 | 20091211 | 20050617 | 20060707 | 20071214 | 20071214 | | |
20040312 | 20050401 | 20060526 | 20100813 | 20050722 | 20060915 | 20090313 | 20090313 | | |
| | | | | | | 20100924 | | |
The applied processing methodology allowed estimation of time deformation histories for each of the coherent points identified in the SAR scene.
Figure 6 shows an example of such time series for three points characterized by three different behaviors. Point A in
Figure 6 has almost flat horizontal trend,
i.e., it is representative of a stable area. Point B experiences a subsidence of about –15 mm/year; for this point very small oscillations are visible. Point C shows a higher rate of subsidence, of about –22 mm/year, with similar oscillations as point B, but more pronounced. The observed oscillations seem to be regularly distributed in time and are likely due to seasonal recharge and discharge of the aquifers.
Figure 5.
Surface deformation velocities (deformation rate) estimated by the Multi-Baseline method, applied to (a) ascending and (b) descending ENVISAT data. Known faults and the cracks observed on the field are superimposed.
Figure 5.
Surface deformation velocities (deformation rate) estimated by the Multi-Baseline method, applied to (a) ascending and (b) descending ENVISAT data. Known faults and the cracks observed on the field are superimposed.
Figure 6.
Example of three deformation histories estimated by the Multi-Baseline approach on ENVISAT ascending dataset. Curves A, B, and C are related to the points depicted in the inset (black circles).
Figure 6.
Example of three deformation histories estimated by the Multi-Baseline approach on ENVISAT ascending dataset. Curves A, B, and C are related to the points depicted in the inset (black circles).
The vertical and east-west horizontal components of the deformation have been calculated combining the ascending and descending ENVISAT data (
Figure 7). It is well know that space-borne SAR, flying on a quasi-polar orbit, are very weakly sensitive to the north-south component of the deformation (about 8% for ENVISAT) [
42]. Thus, this component is always neglected during displacement estimation. Only a very low percentage of the north-south movement can be detected. These maps confirm that on the Eastern side of CG (
i.e., to the east of the observed fissures), the ground is almost stable, with horizontal and vertical deformation rates close to zero. Moving towards west, a sudden increase of the deformation is observed, crossing the fissure alignment. The subsidence rate reaches values of about 25 mm/year (
Figure 7a), together with eastward oriented horizontal movements with a magnitude of about 10–15 mm/year,
i.e., one half of the subsidence rate (
Figure 7b).
Figure 7.
(a) Vertical deformation rate (negative values mean subsidence) and (b) east-west deformation rate (positive values mean eastward movement) computed from ascending and descending ENVISAT data.
Figure 7.
(a) Vertical deformation rate (negative values mean subsidence) and (b) east-west deformation rate (positive values mean eastward movement) computed from ascending and descending ENVISAT data.
Besides the multi-temporal analysis of ENVISAT images, additional SAR data were exploited to further investigate the pre-event deformation prior to the ruptures occurrence. Three descending SAR images from the Canadian RADARSAT-2 mission (6 March, 28 June and 6 September, 2012,
Table 2) allowed capturing and measuring the displacement between March and September, 2012, just before the event. We calculated two interferograms by paring 6 March–26 June, and 26 June–6 September images. The resulting deformation maps are reported in
Figure 8. Both maps show subsidence located in the same area detected by multi-temporal analysis of ENVISAT data. Also in these two cases, the deformation is confined within the alignment marked by the fissures, with a clear increase in deformation with a spatial extent of 4 km
2. Indeed the 26 June–6 September interferogram shows a larger displaced area with respect to the 6 March–26 June map, and in addition, subsidence increase toward NNW with values reaching 0.5–1.0 cm. Such discrepancy in the observed displacements between the two images is probably caused by the seasonal recharge of the groundwater, which is clearly visible in the slight oscillations of the displacement of point C in
Figure 6.
The rupture event was also analyzed by exploiting RADARSAT-2 images. In this case, we used the additional SAR image acquired on 11 December 2012 (
Table 2). This data, combined with the 6 September, image allowed estimation of the co-event deformation (
Figure 9). The deformed areas (
Figure 9) are very similar to the one already observed in the previous interferometric analyses. In addition, a more pronounced subsidence of 1.5 cm is visible in the proximity of the surface ruptures. Furthermore, the loss of coherence immediately at NW of the fissures is consistent with the very large deformations observed during the field surveys (
Table 1). It is worth noting that the image pair used to estimate the co-event deformation has a time lapse of about three months. Therefore, the resulting displacement cannot be considered strictly a co-event one, but additional post-event deformation is present. Indeed, this can be observed looking at the northern section of the subsidence where the values are similar to the pre-event interferogram (
Figure 8b),
i.e., more regular in time.
Figure 8.
Pre event deformation estimated by exploiting 2012 RADARSAT-2 images: (a) 6 March–26 June differential interferogram, and (b) 26 June–6 September differential interferogram.
Figure 8.
Pre event deformation estimated by exploiting 2012 RADARSAT-2 images: (a) 6 March–26 June differential interferogram, and (b) 26 June–6 September differential interferogram.
Figure 9.
Co-event deformation by means of DInSAR applied to RADARSAT-2 data (6 September–11 December 2012). The maximum measured deformation in LOS is about 1.5 cm, located close to the surface ruptures and in the northern area of CG. The loss of coherence, visible in proximity of the fissures, can be due to a large deformation.
Figure 9.
Co-event deformation by means of DInSAR applied to RADARSAT-2 data (6 September–11 December 2012). The maximum measured deformation in LOS is about 1.5 cm, located close to the surface ruptures and in the northern area of CG. The loss of coherence, visible in proximity of the fissures, can be due to a large deformation.
5. Discussion
ENVISAT time series from 2003 to 2010 (
Figure 5 and
Figure 7) clearly show that the subsidence rate affecting the Eastern side of the town is approximately constant in time, with ground velocities increasing from 0 mm/year east of the fissures alignment, up to 25 mm/year northwest of CG. The linearity of the detected deformation likely reflects a consistency in the cause of subsidence. A seasonal oscillation can be noted for points B and C in
Figure 6, related to the dry and rainy seasons (October to January is the dry season while most of precipitations are concentrated from June to September).
Artificial loading (due to urbanization), natural consolidation, earthquakes, or groundwater withdrawal in general may cause the subsidence phenomenon. Such events may act separately or in combination.
It is unlikely that artificial loading produced by urbanization caused the observed subsidence. Urbanization generally results in subsidence after the load has been applied for at least two decades [
43]. A comparison between satellite images of the urban area of CG shows that the urban layout has not changed from 2005 to 2013, suggesting that there is no direct connection with the detected subsidence and fissuring. Furthermore, natural processes such as isostatic sediment loading or consolidation of weak quaternary deposits, which presents typical rates of few millimeters per year [
44], cannot explain the detected rapid subsidence rates.
The hydraulic erosion or piping may act as the main process driving the fissuring of September 2012 in CG [
37], but a careful analysis of daily rainfall estimation by Tropical Rainfall Measuring Mission (TRMM) over the CG area shows no heavy precipitations or storms prior or during the subsidence event (
Figure 10). Because of these reasons, we exclude hydraulic erosion as the main process generating ground deformations.
Figure 10.
Rainfall estimates from TRMM (Tropical Rainfall Measuring Mission—NASA) 3B42 Daily ver.007 data for 2012.
Figure 10.
Rainfall estimates from TRMM (Tropical Rainfall Measuring Mission—NASA) 3B42 Daily ver.007 data for 2012.
The data discussed in
Section 3 and
Section 4 show that ground subsidence extends over a large area and the fissures are aligned NE–SW, suggesting a geometric control of the deep buried structures and the dislocated bedrock. In fact, the NCG bounding faults, prolonged southward under the urbanized area (
Figure 2), clearly overlap with the fissures opened on September 2012. These faults probably act as a structural control on the detected subsidence, even if no tectonic deformation occurred, as showed by the lack of any seismic shaking during the 2012 ground deformations. Furthermore, similar observations on the role that existing faults play on subsidence partitioning have already been observed in Mexico [
17,
23].
Land subsidence caused by compaction of aquifer systems is a worldwide problem in agricultural and urban areas heavily dependent on groundwater supplies. In some places, this kind of land subsidence is associated with structural faults, generating fissures and surface faults due to vertical differential compaction of lacustrine and/or fluvio-lacustrine sediments overlying the faults [
45,
46]. This phenomenon is called Subsidence-Creep-Fault Process (SCFP) [
18,
46] and is characterized by different stages. The groundwater extraction from the pore spaces in unconsolidated sandy to gravelly aquifers causes a lowering of pore-water pressure. This results in an increase of the effective stress in the high-permeability low-compressibility coarse-grained aquifers and a time-dependent pore-pressure reduction in the low-permeability high-compressibility fine-grained aquitards [
47]. The reduction in pore-water pressure produces an increase in the overburden stress, causing the immediate compaction of the soil. If the stress increment due to the groundwater depletion is larger than the preconsolidation stress (
i.e., stress ever experienced by these sediments) the deformation is irreversible as it is caused by the non-reversible grain rearrangement of the sediments. Generally, the compaction of thecoarse-grained soils, which constitute the aquifer, is negligible. If the aquifer has silt and clay beds (aquitards) within, the lowered pore-water pressure in the sand and gravel causes the slow drainage of water from the pore spaces in the silt and clay beds, allowing the fine-grained particles to compress or compact. The overall aquifer volume change is due mainly to the compaction of fine-grained sediments. Reaching pore-water pressure equilibrium between aquifers and aquitards may take months or years, and thick clay layers may take hundreds of years to reach equilibrium. Thus, the resulting compaction may continue long after groundwater withdrawals are brought back into equilibrium with groundwater recharge, or cease completely.
In the presence of bedrock dislocated by faults and overlaid by soft sediments, the ground compaction may produce differential settlements depending on the different thickness of sediments, causing tension and fissuring along the surface projection of the buried fault planes. If the process extends over time, one of the borders of the soil fissure sinks (surface fault), generating a scarp at the surface projection of the buried fault planes.
The SCFP is compatible with the subsidence pattern observed in CG by SAR data and with the development of the fissures of September 2012. The observed high gradients of ground velocities occurring on a narrow zone (<1 km wide) suggest that there is a sharp transition between laterally adjacent lithological units juxtaposed by faulting, which present different thickness, compressibility and permeability. In fact, the aquifer of CG is formed by alternating levels of sandy to gravelly layers (acquifers) and silty to clayey levels (aquitards) [
32,
37]. Furthermore, an important drop of piezometric water level has been reported in CG and in other Mexican cities [
45,
48]. In particular, CG has experienced a reduction in the phreatic level of about 69 cm per year [
37].
The connection between SCFP and the development of fissures at ground has been verified with a numerical model of the ground subsidence using the finite element method. A 2D numerical model crossing CG from NW to SE has been developed (
Figure 2c). The fully coupled approach proposed by Biot [
49] has been considered in the model in order to assess the deformation of porous media that results from fluid withdrawals. The model consists of two materials: the rocky basement, assumed as infinitely stiff and impermeable; and the aquifer body, supposed homogeneous and isotropic. An elastic constitutive model has been assumed for the aquifer body, with the geotechnical parameters listed in
Table 3 [
50].
Table 3.
Geotechnical parameters for the finite element model of
Figure 2c
Table 3.
Geotechnical parameters for the finite element model of Figure 2c
Material | ρ (kg/m3) | E (Pa) | ν | k(m2) | e |
---|
Basement | 2400 | - | - | - | - |
Quaternary soil | 1520 | 6e7 | 0.26 | 1e8 | 0.4 |
The model is vertically fixed at the bottom and horizontally sideways. The only force acting is the gravity load. This simplified model is useful for the analysis of the stresses and deformations into the soil mass that caused the surface cracks. Of course, a multi-layered aquifer system is the best solution in order to best fit the observed displacements, but this kind of modeling is not possible at this time because of the lack of data regarding the complex aquifer system, the exact geometry of the buried faults and the amount of the extracted water.
The modeling has been performed in two steps. At first stage, a physically valid distribution of in-situ stress state is obtained by applying the gravity load. The groundwater table is horizontal at the ground level and pore pressures present a hydrostatic stress distribution. In the second stage, the water table is lowered in order to simulate the displacements observed by the ENVISAT data during 2003–2010. Since the information regarding the amount of the extracted water during this period is missing, we assumed a mean groundwater loss of 67 cm/year, for a total reduction of about five meters in six years. The selected rate corresponds to the annual phreatic surface reduction observed until 2004 [
37]. This assumption is valid because generally the water demand tends to increase increases following the population growth. Certainly, there is the possibility that between 26 June–2012 and 6 September 2012 a strong variation of the water table has happened, causing the deformation pattern in
Figure 8b. Anyway this aspect enforces the hypothesis of a connection between groundwater changes and the ground fissures. A first comparison between observed and calculated LOS ENVISAT displacement profiles at 2010 (
Figure 11a,b) shows a fairly good agreement, despite the simplicity of the adopted model. The displacements are overestimated to the eastern part.
Such discrepancy is probably due to the presence of different lithology, which has different subsidence potentials [
50,
51], or because the groundwater depletion is uneven, thus it does not affect the east side of CG. Moreover, the presence of buried faults, which act as a barrier to groundwater flow, may produce differential subsidence [
45], but such feature is not introduced in our model because of the lack of field data.
Figure 11.
Comparison between calculated and observed LOS displacements considering the cumulated displacements of (a) ascending and (b) descending ENVISAT data. (c) Vertical and (d) horizontal displacement contours, together with the resultant displacement vectors. (e) Horizontal and vertical displacement profiles at the ground level. (f) Horizontal and vertical gradient at the ground level. (g) Contour of the horizontal strains after the lowering of the water table.
Figure 11.
Comparison between calculated and observed LOS displacements considering the cumulated displacements of (a) ascending and (b) descending ENVISAT data. (c) Vertical and (d) horizontal displacement contours, together with the resultant displacement vectors. (e) Horizontal and vertical displacement profiles at the ground level. (f) Horizontal and vertical gradient at the ground level. (g) Contour of the horizontal strains after the lowering of the water table.
In order to explain the features of the observed displacement and the development of ground fissures, the results of the numerical model are shown in terms of horizontal and vertical displacement profiles and contours, and horizontal strains.
Generally, a uniform lowering of the water table, as we modeled in this case, produces only vertical displacements when the thickness of the compacting stratum is constant. In presence of sudden discontinuities and/or gradual changes in sediment thickness, the magnitude of vertical displacements depends strictly on the thickness of the compacting soil, giving rise to the development of a complex pattern. For the CG model, the vertical displacements increase towards NW, following the higher thickness of the compacting stratum (
Figure 11c) as it can be clearly seen by observing the vertical ground displacement profile in
Figure 11e.
Moreover, a change in the thickness produces some local horizontal movements. The soil mass in the ticker part, subjected to higher vertical deformation, pulls the soil in the thinner part towards the deeper side of the profile; this effect causes the horizontal movement of the soil particles. In fact, the modeled thickness heterogeneity causes horizontal displacements close to the surface. More in detail, SE directed horizontal movements develop between 0 and 1250 meters (
Figure 11d,e). These displacements are in agreement with the observed horizontal displacements in
Figure 7b and are mainly driven by the gradual increment of the sediment thickness towards SE. The difference in magnitude between the observed and computed horizontal displacements (
Figure 11e and
Figure 7b) depends strictly on the geometry of the bedrock, which probably is more complex than the one we modeled, and also on the uneven groundwater depletion. Such horizontal displacements are distributed over a large area thus they do not develop high horizontal strains (
Figure 11e). Moreover, larger values of NW directed horizontal displacements develop close to the surface projection of the buried fault at about 1200 meters (
Figure 11d,e).
The general displacement is well represented in
Figure 11c,d. At the left of the model (
i.e., at 0 m), the vectors are oriented downward because of the imposed constraint on the horizontal displacements. This is an artifact of the model, which produces zero horizontal displacements. Proceeding to the right, the vectors are slightly rotated towards SE, until they rotate towards NW close to the buried fault. This displacement pattern justifies the development of ground fissures. In fact, the modeled horizontal and vertical displacement profiles (
Figure 11e) show high gradients in proximity of the cracks (
Figure 11f), where the thickness of the compressible stratum changes because of the buried fault. In particular, maximum vertical displacements in this area are about four times higher than horizontal displacements, suggesting a dominant role of the vertical displacements in the development of fractures. This result is consistent with the geometry of the observed fissures, which show a main vertical offset.
Because the horizontal particle movement is not uneven along the surface, it causes horizontal strain, and the possible generation of ground fissures. In fact, in the area between 0 and 1000 meters, the horizontal strains are generally positive, indicating compression. Contrariwise, maximum tensile horizontal strains are localized in a narrow zone at the surface (at about 1200 meters in
Figure 10f), coinciding with the position of the observed ground cracks. This suggests that ground fissures are mainly driven by vertical displacements, but develop in narrow zones characterized by high values of tensile strains. Such phenomenon could be enhanced by the presence of faults with different mechanical and hydraulic properties with respect to the neighboring soils, and by the presence of a superficial brittle unsaturated zone [
52,
53].