Study Area

The Huehuetán River basin is located in the southern part of the state of Chiapas, a region known as the Coast of Chiapas. It covers an area of 319.27 km² from the upper reaches to the Huehuetan hydrometric station located at coordinates 15°05’56” N and 92°24’02” W (Figure 1). The basin was delimited using the Soil and Water Assessment Tool (SWAT) model (^{Arnold et al., 1998}) and a Digital Elevation Model (DEM) with a 10-m spatial resolution and the drainage network from the Instituto Nacional de Estadística y Geografìa (Mexico’s National Institute of Statistics and Geography) (^{INEGI, 2012}).

Figure 1 Location of the Huehuetán River basin and subbasins. 

This basin is part of the Sierra Madre de Chiapas. It is characterized by a rugged relief in the upper portion, with an altitudinal elevation ranging from 19 to 2 690 meters a.s.l. (Figure 2a); 83.15% of its area has hillslopes greater than 10% (Figure 2b). These variations play differential roles during heavy short-lasting rains, producing large volumes of runoff and sediments from landslides due to the saturation of soil pores and erosion in hillslopes and channels (^{Mora-Ortiz and Rojas-González, 2012}).

Figure 2 Spatial distribution of relief in the Huehuetán River basin. (a) Elevations; (b) Hillslope. 

Mean annual precipitation in the study area was obtained from the Normal Weather Conditions recorded by the National Meteorological Service for the period 1981-2010 (^{SMN, 2022}); this varies between 2 772 mm and 4 654 mm, mainly between May and October. Spatially, its volume is concentrated in the central part because the Sierra Madre de Chiapas is a barrier that makes clouds discharge water on the hillslope facing the Pacific Ocean (Figure 3a). Based on INEGI Series II, the main soil groups are Lithosols, with 36.63%, followed by Andosols, with 28.71%; the remaining 34.66% is distributed between Cambisol and Luvisol soils (^{INEGI, 2006}).

Figure 3 Spatial distribution in the Huehuetán River basin. (a) Mean annual precipitation; (b) geological units. 

The geological features are associated with igneous and sedimentary rocks (Figure 3b), with the former occupying just over 85.05% of the surface area (^{SGM, 1982}). Sedimentary rocks are located mainly in lower areas, particularly on the Huehuetán River channel and margins, associated with surface runoff from upper areas and the detachment, transport, and deposition of sediments in areas with more gentle hillslopes. The interaction between precipitation, topography, soil type, geology, and plant cover directly influences the spatial distribution of landslides (^{Lin and Wang, 2018}).

SINMAP Model Conceptual Framework

The SINMAP (Stability Index Mapping) model combines a uniform-regime hydrological model sensitive to relief variations with an infinite slope stability model (^{Pack et al., 1998}, ²⁰⁰⁵). This model requires a DEM, mechanical and hydrogeological properties of the study area, and an inventory of landslides (^{Tarolli and Tarbuton, 2006}). With this information, SINMAP classifies land stability based on the spatial variation of the hillslope and the watershed (^{Michel et al., 2014}).

For the infinite slope stability model (Figure 4), SINMAP is based on the Mohr-Coulomb law indicating that, at the time of the failure, the shear stress (τ) generated by the descending component of soil weight is equal to the resistance caused by soil cohesion (c) and the resistance to friction due to effective normal stress across the failure plane (^{Michel et al., 2014}), Equation 1.

Source: adapted from ^{Michel et al. (2014)}.

Figure 4 Representation of the infinite slope model. P is the weight of the sliding ground block. 

Where: σ It is normal stress (N m^-2), u is the pore pressure that opposes the normal stress (N m^-2), and ∅ is the angle of internal friction of soil (degrees).

However, ^{Selby (1993)} modified Equation 1 to consider the presence of roots and the level of the groundwater table (Equation 2):

Where: ρ_s is wet soil density (kg m^-3), z is soil depth (m), θ is hillslope (degrees), ρ_w is bulk water density (kg m^-3), h is groundwater table depth (m), c_r is root cohesion (N m^-2), c_s is soil cohesion (N m^-2), and g is gravitational acceleration (m s^-2).

The Safety Factor (SF) is estimated by dividing the right side of Equation 2, representing the soil structural forces (stabilizing forces), by the left side, which corresponds to the destabilizing forces (Equation 3).

In Equation 3, SF = 1 indicates a state of equilibrium; if SF < 1, a hillslope failure occurs; if SF > 1, the hillslope is stable. This value does not indicate absolute stability or instability. Stability increases directly with SF values (^{Selby, 1993}). On the other hand, hillslope stability/instability is closely related to the hydrological conditions of the study area; therefore, SINMAP requires a hydrological model to estimate soil moisture and its relationship with hillslope.

According to Figure 5, α is the area of flow contribution or accumulation (m²), b is contour length (m), and q is the uniform recharge rate (m d^-1). Therefore, defining moisture as the portion of saturated soil based on uniform recharge conditions, ^{O’Loughlin (1986)} proposed a relationship between water entering the soil as uniform or constant recharge and water leaving the soil by saturation through Equations 4 and 5.

Source: ^{Michel et al. (2014)}.

Figure 5 Representation of the hydrological model used by SINMAP. 

Where: W is soil moisture (m m^-1), T is soil transmissivity (m² d^-1), and K_s is the saturated hydraulic conductivity (m d^-1) considered homogeneous through the soil profile.

Based on the mechanical and hydrogeological properties, the SINMAP model calculates the Stability Index (SI), defined as the probability that a site will be stable assuming a uniform distribution of the parameters within the established thresholds (^{Pack et al., 1998}, ²⁰⁰⁵). The SI value normally ranges from zero (most unstable) to one (slightly unstable); however, in places where parameter values allow stability, SI values are greater than one, leading to stability (^{Michel et al., 2014}). Table 1 shows the SINMAP model stability classes.

The SINMAP model calculates the spatial distribution of SI and SSI; for the former, it is based on the calculation of the SF for each DEM pixel. SF is defined as the relationship between stabilizing (soil cohesion and friction) and destabilizing (force of gravity, mainly) forces (^{Pack et al., 2005}). Thus, Equation 6 is obtained by combining the infinite slope stability model and the uniform-regime hydrological model, corresponding to the SF used in the SINMAP model (^{Pack et al., 2005}; ^{Pradhan and Kim, 2015}).

Where: c_α is the dimensionless cohesion parameter (ca=cr+csρs∙g∙z∙cos⁡θ), and r is the ratio of water density to soil density (r=ρwρs).

Due to the spatial and temporal variability of rainfall and soil type (^{Pack et al., 1998}, ²⁰⁰⁵; ^{Pradhan and Kim, 2015}), the SINMAP model works with minimum and maximum limits of the parameters; assuming the following:

Estimation of indexes

The DEM was hydrologically corrected using the ArcGIS^® Fill tool, which allows for eliminating model depressions so that each cell has a defined flow direction, eventually reaching lower elevations when moving from one cell to another (^{Barnes et al., 2014}; ^{Survila et al., 2016}; ^{Wang and Liu, 2006}); then the hillslope, flow direction, and contribution area were obtained. The parameters for each calibration region of the SINMAP model are shown in Table 2, obtained from field data collection by ^{Velescu (2009)} and laboratory tests. Plant cover was obtained from the association between calibration regions and land use of INEGI’s Series VI (^{INEGI, 2017}). The limit values of the T/R parameter associated with soil depth were determined with the average values of the landslide depth identified in the field, which varied between 6 and 14 m.

The two calibration regions defined according to the basin geology allowed the application of the SINMAP model to estimate SI and SSI. Subsequently, a reclassification was performed to determine the spatial variation of the index classes in the delimited subbasins.

Landslide inventory

The landslide inventory was carried out through field trips from 2005 to 2021, recording the location of each landslide and, subsequently, its spatial characterization based on Landsat-7 satellite images, Sentinel-2A, and historical images from Google Earth^®. Satellite images were obtained from the United States Geological Survey (USGS) official website, https://earthexplorer.usgs.gov/, with the processing levels 1TP for Landsat and 2A for Sentinel. During the first years of the inventory, landslides were recorded through their spatial location with GPS equipment supported by satellite imagery; their spatial distribution in the basin was known. This methodology was similar to the approach used by ^{Alanis-Anaya et al. (2017)}, although these researchers used SPOT images and orthophotos for improved detail.

A total of 281 landslide sites were identified for the period 2005 to 2021, mostly located in the upper part of the basin. The main hydrometeorological events that led to landslides included Hurricane Stan in October 2005, with intense rainfall that reached a 100-year return period over six days (^{Arellano-Monterrosas, 2010}; ^{Murcia and Macías, 2009}). Figure 6 shows a landslide of a surface of approximately 1 950 m² and an average depth of 14 m.

Figure 6 Landslide identified in the upper part of the Huehuetán river basin. 

Evaluation of the SINMAP Model Performance

The accuracy of the map of unstable zones produced with the SINMAP model was evaluated using the Receiver Operating Characteristic (ROC) Curve (^{Pradhan and Kim, 2015}), an effective method for representing the quality of probability detection and prognosis systems (^{Svets, 1988}).

The ROC curve method provides a curve given by a binary classification confusion matrix according to four possible outcomes (Table 3). These outcomes are derived by comparing the model results with the landslide inventory, examining the rates of true positives and false negatives (^{Moresi et al., 2020}).

The area under the ROC curve (AUC) represents the accuracy evaluation index of the model, ranging from 0.5 (random prediction, represented by a diagonal straight line) to 1.0 (perfect prediction) and can be used for model comparison (^{Hosmer and Lemeshow, 2005}), This classification is shown in Table 4.

According to ^{Pradhan and Kim (2015)} regarding physical models, the stability index (SI) is reclassified into two classes (stable, SI > 1.0; unstable, SI ≤ 1.0). It is then superimposed with the landslide inventory to calculate the number of land units with the presence or absence of landslides; afterward, the contingency table is derived for the specific cut when crossing susceptibility classes and landslide presence/absence.

Stability index

Based on the spatial distribution of landslides, all were found to occur in areas with intrusive and extrusive igneous rocks; in contrast, areas with sedimentary rocks have no records of landslides (Figure 7). A high correlation was observed between areas with steep hillslopes and the presence of landslides, consistent with ^{Fabbri et al. (2003)}, who determined that topography is one of the key elements defining the location of landslides.

Figure 7 Landslide inventory in the Huehuetán River basin. 

The Stability Index (SI) generated by the SINMAP model classified the Huehuetán River basin into six classes (Figure 8). Each class represents the range of SI values, indicating the probability that each pixel evaluated will get an SF > 1 (^{Michel et al., 2014}). Table 5 shows that the SINMAP model classified 24.36% of the area (77.76 km²) as unstable and 75.64% (241.51 km²) as stable. Previously, ^{Arellano-Monterrosas (2012)} determined the vulnerability to landslides in terms of maximum precipitation over 24 hours, soil porosity, hillslope, current land use, and population density; as a result, a surface area of 35.42% of the basin was classified as having high to very high vulnerability. On the other hand, ^{Bischetti and Chiaradia (2010)} conducted a study in forest areas and found that 17.3% and 82.7% of the study area was unstable, respectively, emphasizing the effect of soil and root cohesion on landslides.

Figure 8 Spatial variation of the Stability Index obtained with SINMAP. 

A total of 229 landslides occurred in the area classified as unstable by SINMAP (Figure 8), and 52 landslides in the area classified as stable. This finding suggests an excessive number of landslides in the stable area; however, their area is significantly smaller, 0.2 km² (20 ha), while the area with landslides in the unstable area was 1.19 km² (119.22 ha). The main difference is that the areas with higher stability determined the presence of smaller landslides than those in unstable areas, where landslides with areas greater than one hectare were recorded.

Given the conditions in the basin related to topography and the occurrence of heavy rains, precipitation was one of the main factors identified as a trigger in most landslides. This finding is consistent with ^{Yilmaz and Keskin (2009)} and ^{Nery and Vieira (2015)}, who indicate that pore pressure and air displacement through the infiltrated water increase during heavy rains of long duration, which increase soil weight and reduce shear stress, leading to landslides. Besides, ^{Zhuang et al. (2017)} analyzed precipitation scenarios and concluded that soil stability in hillslope areas is sensitive to prolonged periods of heavy rainfall. Additionally, ^{Zúñiga and Magaña (2018)} found a trend of an increasing number of events of 100 mm of precipitation accumulated over three days on the Chiapas coast in recent decades, implying a higher incidence of landslides in the region.

The susceptibility to landslides is associated with steep hillslopes, relief, and soil characteristics, magnified by extreme hydrometeorological events that cause groundwater accumulation in soil, triggering mass movements. ^{Sengupta et al. (2010)} found that landslides occur from 250 mm of rain accumulated over 15 days; in the Huehuetan river basin, ^{Murcia and Macias (2009)} reported 242 mm of rainfall on 04 October, which produced landslides in the upper part.

The proximity of the Tacana volcano influences the relief of the Huehuetán river basin; as a result, the soil properties in the study area are associated with processes that have occurred over millions of years, giving rise to a steep relief and landslide susceptibility in the upper subbasins.

Currently, 75.64% of the basin area is stable (Table 5). However, hillslope stability is strongly associated with relief; thus, areas of higher stability are located in the lower part, while more unstable areas are found in the upper part due to the rugged relief (Figure 8). An issue in unstable areas is that heavy rainfall over a short period combined with a decreased plant cover leads to higher landslide susceptibility by saturation of soil pores (^{Collins and Znidarcic, 2005}; ^{Guns and Vanacker, 2013}).

The most unstable subbasins are 7 and 8 (Table 6), located in the upper part. The SI value by subbasin is important for prioritization, as it determines the level of soil degradation due to relief conditions and geotechnical properties. In addition, they allow for a prioritization of subbasins according to their susceptibility to landslides.

Soil saturation index

As regards the distribution of soil saturation, saturated areas are located in the lower part, close to stream channels. In the upper part, the association between hillslope, precipitation, plant cover, and soil type leads to increased susceptibility to landslides (Figure 9). Besides, in the high part, no saturated areas are observed on hillslopes, so surface runoff and sediment transport have a greater impact. On the other hand, in the lower part, there are areas where runoff and sediments accumulate due to the slight hillslopes in channels, favoring sediment deposition and causing a significant reduction of the hydraulic capacity of the channel (^{Pérez-Nieto et al., 2012}); as a result, channels tend to overflow due to the rising water level during heavy rains, causing flooding.

Figure 9 Spatial distribution of the Soil Saturation Index obtained with SINMAP. 

The main driver of instability in upper hillslopes is the combination of soil saturation, geology, hillslope, and plant cover (Table 7). According to ^{Keles and Nefeslioglu (2021)}, dense vegetation cover in temperate and tropical forests is a key element for decreasing the spatial distribution of unstable areas and functions as a natural mitigation measure for landslides. Although ^{Bruijnzeel (2004)} mentions that plant cover can prevent erosion and surface landslides, deep landslides (>3 m) are associated with geological and climatic factors. For their part, ^{Alanís-Anaya et al. (2017)} indicate that all areas covered with natural vegetation will be associated with areas of instability at some point from the effects of changes in land use.

Most of the surface area in the basin is classified in the moist class, followed by the low-humidity class, and 19.12% is in the saturated zone. However, the upper subbasins (5, 6, 7, and 8) have a greater influence on the degree of soil saturation when combined with heavy rains, eroded soils, degraded plant cover — areas subjected to changes in land and vegetation use, mainly for seasonal agriculture, coffee plantations, urban growth. and areas devoid of vegetation — (^{Pérez-Nieto et al., 2012}; ^{Revenga et al., 1999}), and areas susceptible to landslides, leading to hillslope instability by increasing the water content in soil pores, reducing shear stress, and ultimately generating landslides.

SINMAP is an adequate model that can be easily implemented in comparison with other methodologies such as uncertainty analysis, neural networks, fuzzy logic, and logistic regression (^{Huabin et al., 2005}). The data required in SINMAP are easy to obtain, as this model is threshold-based (^{Michel et al.2015}); in contrast, the other methods require detailed field information and complex laboratory tests, which are more time- and resource-consuming to yield results similar to those of SINMAP.

Evaluation of the SINMAP Model

To determine the spatial distribution of landslide stability, SI was reclassified based on the SF classification ranges: classes with an SF ≤ 1.0 were considered for unstable areas, and classes with SF > 1.0 for stable areas (^{Michel et al., 2014}). The success rate evaluates the number of landslides used in the model that are successfully captured by the susceptibility map, thus representing a measure of the model efficiency.

The SINMAP map of landslide susceptibility was validated with the landslide inventory, yielding an accuracy of 85.60% using the ROC curve (Figure 10). Based on these results, the SINMAP model can be regarded as an excellent tool for predicting landslide areas in the Huehuetán River basin.

Figure 10 Results of the ROC curve for the map of landslide susceptibility.