Determination of Susceptibility to the Generation of Discontinuities Related to Land Subsidence Using the Frequency Ratio Method in the City of Aguascalientes, Mexico

Abstract

Land subsidence in the Aguascalientes Valley, documented since the 1980s, has developed a large number of discontinuities that damage infrastructure. There is currently no methodology to accurately predict the site and time at which a discontinuity will occur, making it difficult to make decisions in urban planning or risk management. However, it is possible to determine the susceptibility of an area to the generation of fractures based on the factors associated with their formation. This study presents a zoning method based on the ground failure susceptibility index (GFSI) in the city of Aguascalientes, using the frequency ratio (FR) method and employing the depth of the basement, the subsidence rate, the subsidence gradient, and the groundwater level drawdown as variables. The zoning method included three categories of land subsidence susceptibility to fracturing, moderate, high, and very high, which were divided using the first (3.76) and second (4.24) quartiles of the GFSI. The zoning method was created with the discontinuities reported in 2010 and was validated with data from 2022. The results obtained show that 11.19% of the discontinuities developed between 2010 and 2022 were located in a zone of moderate susceptibility, 41.97% were located in a zone of high susceptibility, and 46.87% were located in a zone classified as having very high susceptibility.

1. Introduction

Land subsidence is a growing concern for construction and infrastructure in countries around the world, including Thailand, Italy, China, the United States, Indonesia, India, Iran, Spain, Japan, and Mexico [1]. In Mexico, several cities, including San Luis Potosi, Queretaro, Celaya, Salamanca, Morelia, Tepic, Zamora, Puebla, Toluca, Jocotepec, Ciudad Guzmán, and Aguascalientes, have experienced this problem [2,3,4,5].

Groundwater extraction in the Aguascalientes Valley has caused fractures and active faults in an area of 900 km² as reported in 2012 [6]. By 2022, using velocity maps derived from InSAR images, it was found that subsidence associated with groundwater extraction covered 878 km², indicating that the affected areas had not increased. This regional subsidence can also lead to differential settlements and terrain fractures. In Aguascalientes, reports of active surface discontinuities date back to 1981 [7], suggesting that land subsidence due to the overexploitation of the aquifer began earlier than previously thought.

The nomenclature used in this work to differentiate between the two types of surface discontinuities in the terrain observed in the Aguascalientes Valley is as follows [8]. The term ‘surface fault’ is used when two blocks formed by the same substratum or by different substrata move parallel to the rupture plane. The term ‘earth fissure’ refers to discontinuities whose blocks, generated through the breaking of the substrate, floor, or rock due to tensional stressors, only present displacements perpendicular to the fracture plane. The term ‘ground failure’ can be used to refer to both of the types described above.

Various studies have proposed methods to develop susceptibility maps for land subsidence using different techniques, or land subsidence hazard maps [9], which is a problem affecting construction and infrastructure in many parts of the world [10,11,12,13,14,15,16,17,18]. However, few have tackled the problem of generating or reactivating the faults and fractures associated with subsidence. Some studies have used numerical models to investigate this problem, such as those of Ye et al. (2018), who defined the trace where a fracture can be generated due to high topography in the underlying rock floor of the deformable fills of the aquifer [8], and He et al. (2020), who reproduced the same high topography configuration of the bedrock in a scaled model to investigate under what conditions and how fractures associated with this bedrock configuration occur [19]. Zang et al. (2021) proposed a method for mapping the susceptibility to earth fissuring associated with tectonic forces and anthropogenic activities by integrating the analytical hierarchical process (AHP) method, the area under the curve (AUC), and the certainty factor model (CFM) [20]; this method identifies areas that are more likely to experience earth fissures based on factors such as land use, soil type, and topography. Risk management against land-subsidence-induced fractures is a crucial challenge. However, the current methods for determining susceptibility rely on regional geological information and observations, such as the area affected by subsidence, active faults, and geological formations.

While these factors are useful, they do not directly address the mechanisms that generate fractures. To improve fracture prediction, other factors, such as the shape of the bedrock, subsidence gradient, water table drawdown, and aquifer thickness, must be considered [6,21,22,23,24,25]. Although recent studies, such as that of Zang et al. (2021), have proposed methods to map the susceptibility to fractures, their results have not yet been validated. Therefore, future studies should evaluate the validity of these approaches by assessing the occurrence of new fractures in areas identified as having a high susceptibility to fracturing.

This paper presents a method for generating susceptibility maps for ground cracking associated with groundwater extraction. The method uses variables that have been identified in the literature as contributing to ground discontinuities and applies the frequency ratio method in the analysis and correlation. The variables used in the method are specified in the literature. To validate the susceptibility map, variable values from 2007 to 2011 and an inventory of the discontinuities reported in 2010 were used. The results indicate that the fractures generated between 2010 and 2022 are primarily located in areas with high and very high susceptibility.

Description of the Study Area

The study area is located within the Aguascalientes Valley, where the aquifer system is composed of poorly consolidated alluvial and lacustrine sedimentary materials [6]. It is a predominantly urban area within the city of Aguascalientes that covers a total area of 78.33 km². The climatological and hydrogeological conditions of the Aguascalientes Valley are favorable for the development of discontinuities due to subsidence [26] since (a) the aquifer system is mainly composed of poorly consolidated sedimentary fillings; (b) the valley has a semi-arid climate and no surface water deposits, so the aquifer undergoes long periods of groundwater pumping; and (c) differential subsidence occurs due to lateral variations in the thickness of the fillings that form the aquifer. Figure 1 shows the study area within the context of the Aguascalientes Valley with the terrain discontinuities.

According to INEGI [27], the Aguascalientes Valley ranks second in the land subsidence rate among the areas experiencing subsidence in México and third in surface area affected by subsidence, making it one of the most studied places in Mexico for subsidence and cracking [6,7,28,29,30,31,32,33,34].

The generation of new cracks in the Aguascalientes Valley has significantly increased in the last 10 years. In 2013, 207 discontinuities with a total length of 321,159 m were recorded [6]; currently, 802 discontinuities with a total length of 361,938 m are registered [35]. Additionally, different types of discontinuities have been documented in the Aguascalientes Valley, such as faults due to tertiary tectonism that are currently inactive, faults generated by differential land subsidence, faults reactivated by land subsidence, and fractures reactivated by paleochannels [6,29]. For the development of the fracture susceptibility map presented in this study, we only considered discontinuities generated by differential land subsidence and reactive faults due to land subsidence. Other types of discontinuities, such as faults due to tertiary tectonism and fractures reactivated by paleochannels, were not included.

Initially, fault and fracture mapping was only carried out for the city of Aguascalientes [36,37,38]. Currently, mapping covers the entire Aguascalientes Valley and is publicly available online [35]. For the susceptibility analysis presented in this paper, we considered the fault and fracture mapping reported in SOPMA [38] up to 2010, and we used the updated mapping presented in SIFAGG 2022 [35] to validate our results.

2. Materials and Methods

To evaluate the susceptibility of generating new faults and fractures in the study area (Figure 1b), we applied the frequency ratio (FR) method, a statistical technique that uses geographic information. To implement this technique, it is convenient to use geographic information systems. The FR method involves quantifying the magnitude and frequency of the variables associated with the occurrence of the studied phenomenon. The variables are worked in layers in raster images, where each pixel contains information about its spatial location (longitude and latitude) and the value of the variable at that location. We determined the number of classes into which each variable was divided by testing different numbers of classes and attempting to find a trend between the classes and RF. For this study, we divided the variables into three to four classes.

In this study, we used the FR technique to develop an index that identified the combination of variables and their magnitudes that were present in fault and fracture zones. A high FR value indicated the presence of faults and fractures, resulting in a high index value. There were also zones that did not exhibit discontinuities but had a high index value because they fell within classes with high frequency ratios. This suggested that the variables associated with the generation of faults existed but that the faults themselves had not yet developed. These zones were therefore considered prone to cracking.

The data correlation analysis was carried out using the free software Qgis (version 3.18.2). The final result was a susceptibility index which, in this work, was called ground failure susceptibility index (GFSI), which was determined by applying Equations (1) and (2) of the FR technique [39].

The variables used in this study, which, according to the literature, are involved in the process of generating cracking, were (a) the depth of the rock layer underlying the consolidating sediments, (b) the subsidence accumulated during the analysis period, (c) the horizontal subsidence gradient, and (d) the lowering of the water table. Changes in lithology, which is also an important variable for the formation of discontinuities, were not considered in this research, as there was only one type of material within the study area.

The variables used were based on data from 2007 and 2011 and were chosen due to the availability of information on subsidence and water table lowering. This analysis period was selected to produce a susceptibility map that could be verified against records of new fractures developed from 2010 to the present.

$$F{r}_i=\frac{Npix/Npt}{{{\displaystyle \sum }}^{}Npix/{{\displaystyle \sum }}^{}Npt}$$

(1)

where Fr_i is the frequency ratio for the i-th variable, Npix is the number of pixels in the class, ∑Npix is the total number of pixels of a variable in the entire study area, Npt is the number of points defining the trace of fractures in the class, and ∑Npt is the total number of points defining the trace of all fractures in the study area.

To obtain the GFSI, it was first necessary to divide the studied variable into classes; to obtain Npix and ∑Npix; then, to obtain Npt and ∑Npt. Next, Equation (1) was used to obtain the frequency ratio for the variable under study. This procedure was repeated for each of the variables studied.

The FR for each class was calculated, and this value was assigned to all pixels within that class for each of the four variables. The resulting maps for each variable were then used to calculate the GFSI by adding the frequency ratio value for each variable of each pixel as shown in Equation (2). Figure 2 provides a schematic representation of this procedure.

$$GFSI=F{r}_I+F{r}_2+F{r}_3+F{r}_4$$

(2)

The variables used to determine the GFSI (Figure 3) are listed below, and how they were obtained is explained.

(a)

Bedrock depth. The thickness of the fill materials in the Aguascalientes Valley, which have consolidation potential, was obtained through gravimetric surveys [40]. The values of the depth of bedrock were categorized into 3 classes as shown in Figure 3a.

(b)

The LOS deformation accumulated during the analysis period was calculated from the results of previous investigation [40] in which a subsidence velocity map of the Aguascalientes Valley based on 34 ALOS satellite INSAR images was developed using the SBAS technique. The subsidence values were categorized into 4 classes as shown in Figure 3b. Although the variable associated with the generation of fractures is the vertical deformation, for this work, the LOS deformation was used, which was proportional to the real vertical subsidence.

(c)

Zones with differential subsidence were determined based on the subsidence gradient according to [41,42]. The values of differential subsidence were categorized into 3 classes as shown in Figure 3c.

(d)

The water table drawdown was calculated from the difference between the water table level in 2011 and that in 2007; the data was obtained from the official CONAGUA website [43]. The magnitude of the drawdown was classified into 3 classes as shown in Figure 3d. In total, 101 wells within the Aguascalientes Valley were used, the location of which is shown in Figure 4.

The original technique requires the information of the studied phenomenon to be expressed in points. However, the traces of the terrain discontinuities are linear entities; therefore, in this study, they were discretized into points spaced two meters apart to achieve adequate resolution of the discontinuities without duplicating information.

Since the range of values for each variable was different, the number of classes into which they were divided was also different. The images of the variables over the study area were discretized into pixels with a resolution of 30 by 30 m, which is the resolution of the digital elevation model used.

To eliminate errors associated with class size, we identified outliers for each variable. To do so, we determined the first quartiles (Q1), second quartiles (Q2), and third quartiles (Q3) as well as the interquartile range (IQR = Q3 − Q1). We calculated the lower limit with Q1 − 1.5 × IQR and the upper limit with Q3 + 1.5 × IQR. Values that fell outside these limits were considered outliers and were reclassified as extreme values in their respective class. This approach allowed us to create classes between the lower and upper limits for each variable. The data distributions, lower and upper limits, and outliers are shown in Figure 5.

The GFSI map was calculated using variables from 2007 to 2011 and discontinuities reported in 2010. In order to validate the results, this map was evaluated with discontinuities reported in 2022. Figure 6a shows the reported discontinuities in the city of Aguascalientes over different years, and Figure 6b shows the accumulated length of reported discontinuities in the Aguascalientes Valley.

The susceptibility to cracking was classified into three levels: moderate, high, and very high. Typically, studies on hazard zoning use five categories; however, this research only used three levels because the entire study area was located on sediments with high consolidation potential, making the area susceptible to cracking. The zoning method for susceptibility to cracking in this research coincided with the zoning method for subsidence hazards proposed by Pacheco-Martínez et al. in 2015 [40]. This occurred because the study area fell within the moderate and high zones proposed by Pacheco-Martínez. In addition, their study explains that areas with zero hazards are only found in locations where there are no materials with consolidation potential.

The zoning of susceptibility zones was performed based on quartiles of the GFSI value range. The moderate-susceptibility zone corresponded to pixels with GFSI ≤ Q1; the high-susceptibility zone was defined by pixels with Q1 < GFSI ≤ Q2; and finally, the very-high-susceptibility zone was defined by GFSI > Q2.

3. Results

Table 1. Ground failure factors, classes, and weight.

Variable	Class	Class Interval	Number of Cells per Class	Percent	Number of Failure Points per Class	Percent	Frequency Ratio
Bedrock depth	1	≤129.3	43,381	49.71%	11,832	66.24%	1.33
	2	129.3–261.8	32,054	36.73%	5766	32.28%	0.88
	3	≥261.8	11,825	13.55%	264	1.48%	0.11
Accumulated LOS deformation or Subsidence	1	≤−5.90	8729	10.00%	84	0.47%	1.05
	2	5.90–11.0	38,029	43.58%	7558	42.31%	1.29
	3	11.0–16.2	30,292	34.71%	8020	44.90%	0.97
	4	>16.2	10,210	11.70%	2200	12.32%	0.05
Gradient subsidence	1	≤0.06	46,815	53.65%	8084	45.26%	0.84
	2	0.06–0.12	27,851	31.92%	6208	34.76%	1.09
	3	>0.12	12,594	14.43%	3570	19.99%	1.38
Lowering of water table	1	≤7.76	25,305	29.00%	5561	31.13%	1.07
	2	7.76–10.13	36,969	42.37%	8319	46.57%	1.10
	3	>10.13	24,986	28.63%	3982	22.29%	0.78

shows the limit values that define each class for each of the four variables with the corresponding percentages that define the frequency relationship between the fractures and the classes of each variable. The highest correlation was reported for the range of higher subsidence gradient values, with a frequency relationship value of 1.38, and the lowest correlation corresponded with the class with the highest accumulated deformation, with a frequency relationship of 0.05. Figure 7 shows the frequency relationship graphs for each class per variable studied from which the following observations can be made:

(a)

According to the graph in Figure 7a, the sediment thickness variable showed a behavior that was inversely proportional to the frequency ratio, meaning that zones with a lower thickness of fillings had a higher frequency ratio, i.e., they are areas where a greater number of discontinuities have occurred. This is consistent with the conceptual models of fracture generation reported by [6,21,22,24]. These models show that the generation of discontinuities takes place in zones where the underlying bedrock of the aquifer is shallower. This occurs where there is a lateral change in the depth of the aquifer.

(b)

The graph of the accumulated deformation versus the frequency ratio (Figure 7b) suggests that the zone where the highest number of fractures has occurred was not where the deformation was the highest or the lowest; instead, the highest number of fractures occurred in the zone that corresponded to the intermediate values of accumulated deformation. According to the subsidence conceptual models [6,21,22,24], a greater sediment thickness produces greater deformation, and a smaller sediment thickness produces less deformation. Therefore, the highest number of fractures should occur in the zone of transition from the thinner to the thicker sediment layers, i.e., in the subsidence zone with intermediate values as shown in the graph of Figure 7b.

(c)

The graph in Figure 7c, which shows the frequency relationship between the subsidence gradient and the occurrence of fractures, presents a directly proportional relationship. This means that the highest number of fractures was observed in areas with a higher horizontal subsidence gradient, which is consistent with a study carried out in México City [35] in which the authors used the horizontal subsidence gradient as a parameter to identify the fracturing zones.

(d)

The graph in Figure 7d suggests an inversely proportional relationship between the generation of fractures and the lowering of the water table. Figure 7d shows that the highest number of fractures occurred outside the cone of depression, which is consistent with [44,45]. They noted that the greatest subsidence occurs in areas where the lowering of the water table is highest; that fractures occur at the edges of the subsidence zone; and that, therefore, they occur outside the cone of depression, where the water level drop is lower.

4. Discussion

Figure 8 shows a map of the GFSI that defines the level of susceptibility to the generation of fracturing based on the frequency relationship of the four analyzed variables with the presence of cracking. The resulting GFSI values ranged from 2.6 to 5.12, with higher values representing areas with greater susceptibility to cracking.

Figure 9 shows a histogram of the distribution of values obtained from the GFSI, which presented a range from 1.78 to 5.12. The histogram has a negative or left skewness, which means that the majority of the study area was in a zone of very high susceptibility; meanwhile, some areas had medium values, and there were some outliers up to 2.60. The zoning was done as follows: from 0 to 3.76 (Q1), it was considered to be moderately susceptible; from Q1 to Q2 (4.24), it was considered to be high-risk zoning; and values higher than Q2 were considered to represent very high susceptibility. Within the range of GFSI values, there were values below the lower limit that were calculated with Q1 − 1.5 × (Q3 − Q1), which were considered to be outliers. These values represent the pixels with the least susceptibility to cracking, but, being outliers, they were considered to be within the moderate-susceptibility zone.

The zoning method for susceptibility to cracking, developed with the information on cracking that is available up to the year 2010, is presented in Figure 10a. Figure 10b shows the same susceptibility zoning with the updated discontinuities up to the year 2022. The moderate-susceptibility zone represents 22.62% of the total area; meanwhile, the high-susceptibility zone represents 25.87% of the total area, and the very-high-susceptibility zone represents 51.51% of the total area.

After zoning, the percentage of total discontinuities was counted in each zone, first for the year 2010 in order to have a reference and to compare it with the percentage of cracks in each zone in the year 2022. In addition, the percentage of new cracks generated from 2010 to 2022 was counted for each susceptibility zone (Figure 11). In the year 2010, the percentage of total discontinuities in the moderately susceptible zone was 17.76%; in the high-susceptibility zone, it was 34.74%; and, in the very-high-susceptibility zone, it was 47.50% (Figure 11a). For the year 2022, the percentage of total discontinuities in the moderately susceptible zone was 16.47%, that of the high-susceptibility zone was 36.16%, and that of the very-high-susceptibility zone was 47.37% (Figure 11b). Figure 11c shows that, of the new discontinuities generated from 2010 to 2022, 11.19% were generated in the moderately susceptible zone, 41.97% occurred in the high-susceptibility zone, and 48.87% were generated in the very-high-susceptibility zone.

Figure 10 shows the zoning of susceptibility to cracking, which indicates areas with a higher or lower likelihood of new discontinuities occurring based on the analysis of factors related to their genesis. Unlike the previous work [34], which presented the zoning of hazards to the effects of subsidence and cracking for the year it was conducted, this study warns about the possibility of new cracking occurrences in the future.

Some of the discontinuities used in zoning or validation may not be a result of differential subsidence caused by groundwater lowering or lateral variations in the depth of bedrock as previously explained [6,21,22,24]. Other fractures could be associated with the reactivation of paleochannels [6]. This could introduce errors in susceptibility determination. To improve the reliability of the susceptibility map, it is necessary to identify which fractures are due to this mechanism and to exclude them from the calculation.

Figure 6b shows a decrease in the generation of new cracks since 1995. In that year, the study area had 13,990 linear meters of active fractures. By 2004, the number had increased by 119.26% to 30,675 linear meters of discontinuities, with an annualized rate of new fractures of 13.25%. In 2010, the inventory of discontinuities was 35,940 m, which represents an increase of 17.16% in six years, with an annualized rate of 2.86%. By 2022, the number of faults and fractures reported was 44,722 m, an increase of 24.44% in 12 years, with an annualized rate of 2.04%. All percentages were calculated with respect to the previous record.

5. Conclusions

The application of the frequency ratio method in the development of a susceptibility map for groundwater-extraction-induced fracturing in a zone of the city of Aguascalientes as presented in this work shows that the method can be directly applied and supported through geospatial information on the factors recognized in the scientific literature as involved in the genesis of fracturing. The results of the frequency ratio found between each variable and the occurrence of fracturing are consistent with the conceptual models reported in the literature. In particular, fracturing is generated in the shallowest zones of the aquifer, where there is a lateral change in the depth of the rocky bed underlying the fill with consolidation potential, in zones where the highest values of the horizontal subsidence gradient are present and in zones outside the cone of depression.

The susceptibility zoning method that was carried out was categorized into three levels, moderate, high, and very high; the frequency ratio technique was applied to the variables and the fractures that were recorded up to 2010. The qualitative validation of this zoning method was carried out by contrasting the fractures generated from 2010 to 2022 with the obtained zones, finding that 11.19% of the discontinuities developed in areas of moderate susceptibility, 41.97% developed in areas of high susceptibility, and 46.87% developed in areas of very high susceptibility, validating the zoning method.

It is important to consider that the study area was limited to a circle with a radius of 5 km within the urban area of Aguascalientes due to the availability of information on the variables. To produce a susceptibility map that covers the peripheral non-urbanized areas outside the city, where new buildings and urban infrastructure will be constructed, more work will be required to obtain information on the entire subsidence zone.

The influence of the scale and resolution of the variable maps on the resulting susceptibility map should also be considered. More work is required to define the optimal scale and resolution of the variable maps to produce a more accurate susceptibility map that successfully defines the areas where new fractures will occur in the future. In conclusion, this work will be an important reference for the management and mitigation of the hazards associated with the presence of discontinuities in subsidence zones, and it will be a useful tool for urban development decision makers in designing less vulnerable infrastructure and in reducing the risk of fractures associated with subsidence.