A Multivariate Geomorphometric Approach to Prioritize Erosion-Prone Watersheds

Abstract

Soil erosion is considered one of the main degradation processes in ecosystems located in developing countries. In northern Mexico, one of the most important hydrological regions is the Conchos River Basin (CRB) due to its utilization as a runoff source. However, the CRB is subjected to significant erosion processes due to natural and anthropogenic causes. Thus, classifying the CRB’s watersheds based on their erosion susceptibility is of great importance. This study classified and then prioritized the 31 watersheds composing the CRB. For that, multivariate techniques such as principal component analysis (PCA), group analysis (GA), and the ranking methodology known as compound parameter (Cp) were used. After a correlation analysis, the values of 26 from 33 geomorphometric parameters estimated from each watershed served for the evaluation. The PCA defined linear-type parameters as the main source of variability among the watersheds. The GA and the Cp were effective for grouping the watersheds in five groups, and provided the information for the spatial analysis. The GA methodology best classified the watersheds based on the variance of their parameters. The group with the highest prioritization and erosion susceptibility included watersheds RH24Lf, RH24Lb, RH24Nc, and RH24Jb. These watersheds are potential candidates for the implementation of soil conservation practices.

1. Introduction

Soil erosion is considered one of the most important degradation processes in the world [1,2]. The soil resource is limited and its wide use is of utmost importance; it sustains biogeochemical processes and is the habitat for a great diversity of microorganisms [3]. Sustained soil development, conservation, and restoration is one of today’s main challenges for humankind.

Hydric erosive processes affect the fertile soil layer, which is a key factor in the primary production of ecosystems [4]. The production of goods and satisfiers for the population such as wood, food, fiber, fodder, water, and recreational areas, among others, in addition to industrial expansion and the need for infrastructure facilities, have increased land-use/land-cover changes, increasing the pressure over the soil [5]. This has caused experts to pay more attention to the growing trend of soil erosion and the importance of water and soil conservation for achieving sustainable development.

Integrated watershed management is an alternative for soil management [6,7,8]. Watersheds are one of the spatial units that are used for the planning and management of soil resources [9]. Management implies the characterization of the ecosystems inside the watershed and the understanding of the relationships between uplands, lowlands, land use/land cover, geomorphic processes, slope, and soil [10,11]. In watershed management, erosion control is one of the main components [3]. Thus, the hydrological planning and monitoring of a watershed is important for the development of environmental policies [12]. In this sense, the quantification of the watersheds’ characteristics is fundamental to understanding their dynamics and degradation levels. This knowledge serves to define and implement strategies to prevent soil erosive processes and promote the conservation and restoration of watersheds [13].

Morphometry is used in the analysis of the watershed configuration [14]. Such methodology was developed by Horton [15] and then modified by Strahler [16], and provides information on the behavior of the basin [17]. It is an important tool for identifying and prioritizing eroded watersheds [18].

Nevertheless, monitoring soil erosion in situ is costly and time consuming in large watersheds. Thus, the analysis of geomorphometry is often carried out based on geographic information systems (GIS) [19,20,21,22,23]. On a spatial scale, geomorphometric parameters, e.g., the Gravelius compactness coefficient [24] and elongation ratio [25], among others, are important to know the hydrological configuration of watersheds. The relationships among these parameters are useful for developing hydrological models, which allow prioritizing watersheds based on their condition, such as erosion susceptibility. To determine the aforementioned relationships, statistical methods, such as multivariate techniques, have been widely used worldwide [26,27]. For instance, Gavit et al. [28] characterized 13 geomorphometric parameters in 11 watersheds located in the Godavari river in Maharashta, India. Youssef et al. [29] estimated the erosion risk by using remote sensing technology, GIS, and geomorphometric parameters in 11 watersheds located in Sinai, Egypt. Makwana and Tiwari [30] used seven geomorphometric parameters to characterize 19 watersheds in the region of Gujarat, India. Sharma et al. [31] applied the multivariate technique of principal component analysis (PCA) to 13 geomorphometric parameters from eight watersheds located in the Madhya district, India. Meshram and Sharma [32] and Farhan et al. [33] used PCA to analyze the geomorphometric parameters of a group of watersheds located in the Shakkar Catchment River in India and Jordan.

A large number of shape, relief, and hydrological parameters are associated with watershed geomorphometry [34,35]. The statistical techniques of PCA and group analysis (GA) [36,37,38,39,40,41], as well as multivariate analysis of variance (MANOVA) [41] and the ranking methodology known as compound parameter (Cp) [42] have been widely used in recent years for the analysis of environmental data from watersheds. These techniques assist with analyzing the spatial variability of watersheds, their structure, as well as the relationships existing among them.

The most important basin in the state of Chihuahua, as a runoff source [43], is the Conchos River Basin (CRB). Yet, this basin has experienced water stress conditions during the past years. According to Ordoñez [44], approximately 8000 km² (11.8%) of the basin high lands present strong erosion problems, which could impact waterflow and water quality. In these high lands, deforestation and land-use/land-cover changes had contributed to reduce the amount of infiltrated water, impacting on groundwater availability [43]. In the low lands of the basin, agriculture consumes 90% of the water harvested in the basin. Other consumers include the industrial and domestic sectors [43,45]. In addition, the international water trade between Mexico and the U.S. from 1944 commands Mexico to deliver annually from this basin a total of 432 × 10⁶ m³ of water to the U.S. [46]. Therefore, specific knowledge about water management and the status of the basin’s soil erosion is required to implement strategies to solve water-related problems and to promote the sustainable development of the region.

The objective of this study was: (1) to describe the behavior of the 31 watersheds located along the CRB in the state of Chihuahua, Mexico, based on the values of their geomorphometric parameters; (2) to spatially classify the 31 watersheds into groups; and (3) to prioritize the groups of basins according to their erosion susceptibility. For the grouping, multivariate techniques and the compound parameter (Cp) were used and their efficacy compared.

2. Materials and Methods

2.1. Study Area

The CRB is located in the state of Chihuahua and Durango, Mexico, and is part of the 24th Rio Bravo-Conchos Hydrological Region [47] (Figure 1). The basin has an area of 67,800 km² [48], with an altitudinal distribution ranging from 841 m to 2348 m [49]. It presents a diversity of climates ranging from temperate in the upper, semi-arid in the middle, and arid in the lower parts of the basin [50]. The physiography of the upper basin belongs to the mountainous zone made up of temperate forests dominated by species of pines and oaks. The middle part of the Altiplano or central valleys is made up of transition zones where oaks and bushes are present. Regarding the lower part, it belongs to the arid region and is made up of shrublands and grasslands [51]. The basin has a precipitation regime from June to September, with July and August being the wettest months, and fluctuating from 200 mm to 700 mm [48].

2.2. Data

Data of the CRB was obtained from the online GIS source of CONABIO [52]. Likewise, the data of the 31 watersheds composing the CRB (Figure 2) were obtained from the Watershed Water Flows Simulator [53]. The relief and hydrology type parameters were estimated by processing the necessary data from a Digital Terrain Model (DTM), with a resolution of 15 × 15 m, downloaded from INEGI [54]. The values of the basic parameters from the watersheds were obtained by using the Hydrology tool [55] of ArcMap© 10.3 (ESRI, Redlands, CA, USA; https://wwwesri.com/en-us/home). The values of the shape, relief, and linear type parameters were calculated from the equations listed in Table 1.

Basic parameters include the area (A), perimeter (P), watershed length (Lb²), stream order (u), main channel length (Lc), all channel lengths (Lu), and contour length (Li). The area and perimeter are the most important parameters of the watersheds to understand their hydrological design and reflect the volume of water that can be discharged in a rainfall event [61].

Shape parameters include the Gravelius compactness coefficient (Cc), elongation ratio (Re), shape factor (Rf), elongation index (Ia), unit shape factor (Ru), and circularity ratio (Rc). The Gravelius compactness coefficient is the relationship between the perimeter of the watershed and the perimeter of a circle of area equal to that of the watershed [24]. Low values of the Gravelius compactness coefficient, shape factor, and elongation index indicate a more elongated watershed, while high values correspond to a less elongated watershed. Watersheds with less elongated shapes have a shorter maximum flow duration, while elongated watersheds correspond to watersheds with longer flow duration [24]. The elongation ratio is the diameter of a circle with an area equal to that of the watershed divided by the length of the watershed [32]. Those watersheds with values close to or greater than the unit correspond to flat watersheds, while values that move away from the unit are watersheds with pronounced relief [25]. The shape factor is the relationship between the area and the length of the watershed [15]. The elongation index is also a relationship between the length and width of the watershed [15]. The unit shape factor is the relationship between the length and the area of the watershed. Values less than 2 indicate that they have weak flood discharge periods, while values greater than 2 indicate that their flood discharge is strong [30]. The circularity ratio is defined as the relationship between the area and the perimeter of the watershed. If the circularity of the watershed is low, then discharge will be slow, and the susceptibility to erosion will be low [61].

Relief parameters include the mean watershed slope (J), massivity coefficient (tgα), relief ratio (Rh), relative relief (Rr), and orographic coefficient (Co). The mean watershed slope indicates the degree of the terrain roughness. As the slope increases, the watershed will be prone to erosion. The massivity coefficient represents the relationship between the mean watershed height and its surface area, which is expressed as a percentage. Small values of the massivity coefficient correspond to very mountainous watersheds, while large values correspond to watersheds from moderately mountainous to flat. The relief ratio is directly related to the slope of the currents and the watershed surface, being an indicator of hydrological processes and erosion. The relief ratio, similar to the relative relief, has a direct correlation with the watershed erosion processes [25].

Linear parameters include drainage density (Dd), mean slope of the main channel (j), mean distance (Am), sinuosity of the main channel (Scp), Kirpich concentration time (Tc), average peak flow (Qp), texture ratio (T), rivers frequency (Fu), resistance number (Rn), general flow length (Lo), and drainage intensity (Di). High drainage density values indicate a high current density, and therefore a rapid response to precipitation events [62]. The mean slope of the main channel indicates the slope of the longest channel in the watershed. The high values of this parameter indicate that sediment flow and entrainment will quickly exit the watershed [15]. The sinuosity of the main channel indicates the velocity of flow movement in the channels. The lowest values of sinuosity correspond to channels where the flow travels at greater speed, whereas the channels with the highest values of sinuosity indicate the slow movement of the flow [57]. The Kirpich concentration time is the time when a drop of water falls at the furthest point until it reaches the exit point [58]. Average peak flow is defined as the mean maximum amount of water passing through a specific section [59]. The texture ratio corresponds to the relationship between the total number of streams and the watershed perimeter. Rivers frequency represents the total number of streams of all orders per unit area [15]. The resistance number is used to measure the flood potential of rivers. It has a direct relationship with erosion, where increasing its value would represent an increment in erosion susceptibility [15].

2.3. Watershed’s Description and Classification

Prior to the watersheds’ classification, their geomorphometric parameters were correlated [6,63]. Correlation indicates when part of the information contained in a set of geomorphometric parameters is also contained in the remaining ones [32]. That served to reduce the number of parameters included in the subsequent analysis.

To describe the variability of the geomorphometric parameters, principal component analysis (PCA) was performed. The technique of PCA reduces the data dimensionality, simplifies the dataset, and makes it easier to explain through a set of new principal components (PCs) [64,65]. The first principal component is the linear combination of the original geomorphometric parameters that contributes the most to the total variance; the second principal component, not correlated to the first, contributes the most to the residual variance, and so on until the total variance is analyzed. The PCs are orthogonal variables that could be obtained by multiplying the original variables, which are correlated, with coefficients similar to the ones described in Equation (1):

(1)$z_{ij}=a_{i1}{x}_{1j}+a_{i2}{x}_{2j}+\dots +a_{im}{x}_{mj}$

where z represents the coefficient of the component, a represents the weight of the component, x represents the measured value of the variable, i corresponds to the component number, j represents the sample number, and m represents the total number of variables.

The PCA was applied to the values of geomorphometric parameters to calculate the correlation matrix and to obtain the PCs [66]. Both the correlation analysis and the PCA were performed in the SAS© 9.1.3 software (The SAS Institute Inc., Cary, NC, USA, https://www.sas.com/en_us/home.html). Then, the PCs were mapped for interpretation.

For the classification, the first method used for the watersheds was a hierarchical group analysis, which was based on the Ward criterion [67]. The Ward criterion was applied to the GA by using the square Euclidean distance to explore the clustering of the 31 watersheds. The definition of the watershed groups was performed based on the coefficient of determination (R²) [68]. Finally, a multivariate analysis of variance (MANOVA) was used to analyze whether significant multivariate differences exist between the groups based on the values of their geomorphometric parameters [69].

The second classification method considered in this work was the compound parameter (Cp). Previous research has employed this approach for sustainable watershed planning and management [42]. Linear and shape parameters have been commonly used for this method, whereas the relief and basic parameters were additionally included in this study. Linear geomorphometric parameters have a direct relationship with erosion susceptibility, where high values are more likely to result where high erosion probabilities are present [61,70]. Thus, for watershed classification, the highest value of linear parameters was ranked as 1, the second highest was ranked as 2, and so on. Conversely, shape parameters have an inverse relationship with erosion [61,70], and low values are related to high susceptibility to erosion and vice versa. Then, the lowest value of the shape parameters was ranked as 1, the next lowest value was ranked as 2, and so on. Regarding the relief and basic parameters, the highest value was ranked as 1, the second highest value was ranked as 2, and so on. After this procedure was completed, the ranked values from each watershed were summed and then averaged to produce the Cp of each watershed. This average represents the collective impact of all the parameters, and is calculated according to Equation (2) [42]:

(2)$C_P=\frac{1}{n}{\displaystyle \sum }_{i=1}^{n}R$

where Cp is the compound parameter, R is the ranked value of each parameter, and n is the number of parameters.

Based on the Cp, the highest priority was assigned to the watersheds with the lowest Cp value, the second priority was assigned to the next higher Cp value, and so on. Then, the Cp was classified into five categories or groups of erosion susceptibility. The categories were defined as: Very High (Group 5), High (Group 4), Moderate (Group 3), Low (Group 2), and Very Low priority (Group 1), similar to classifications made in previous studies [71].

2.4. Comparison of the Classification Methods

To compare the grouping methods, the following procedure was followed.

First, an ANOVA was carried out for each geomorphometric parameter (independent variable), and separated for each grouping method. The source of variation or class was considered to be the group. Such analyses served to determine possible significant differences among the groups within each grouping method. After the ANOVA analyses were completed, the grouping method that achieved the highest number of significant p-values (α < 0.05) was considered the most effective for grouping the watersheds of the CRB.

A list of abbreviations can be found in Appendix B (Table A3).

3. Results

3.1. Characterization of Conchos River Basin Watersheds

The values of the basic, shape, relief, and linear-type parameters from the 31 watersheds used in this study are shown in Table A1. Watershed RH24Ia has the smallest area and perimeter with 78.62 km² and 50.59 km, respectively. Meanwhile, the watershed with the largest area and perimeter is RH24Lb, with 5428.29 km² and 640.77 km, respectively. The watershed with the highest stream order is RH24Kb, while watershed RH24Lb has the longest main channel. First-order streams do not have tributaries, and their flow depends on the secondary surface contributions that converge to them [61]. The watershed with the highest number of channels of order one is RH24Ia, whereas the watershed with the lowest number of channels is RH24Jb. The watershed lengths vary from 11.17 km to 133 km. Watershed RH24Kb presents the largest elongation ratio value, indicating that it is the flattest, while watershed RH24Ib has the lowest value for this parameter, indicating that it has the steepest slope [25]. The values of sinuosity of the main channel vary from 0.05 km to 4.1 km. The watershed with the shortest Kirpich concentration time is RH24Kb, while watershed RH24Lb has the longest. The lowest average peak flow value corresponds to watershed RH24Ia, whereas watershed RH24Lb presents the highest. The texture ratio values are between 4.98–75.1, which are considered as moderate to high values; the low values correspond to watershed RH24Na, while the highest value belongs to watershed RH24Lg, so the former is not susceptible to erosion, while the latter is. The watershed with the lowest resistance value is RH24Na, while the watershed with the highest value is RH24Lg.

3.2. Correlations and Principal Component Analyses

The data matrix of the 31 watersheds and the 33 parameters, which included the basic shape, linear, and relief type parameters, were analyzed through a correlation analysis. A set of parameters showing high correlations were identified. From each pair of highly correlated parameters, only one parameter was chosen, and the rest were eliminated to reduce the data dimensionality. After the reduction, the final number of geomorphometric parameters was 26, as listed in Table A2.

The PCA was performed on the 26 parameters and showed that the first five principal components were the most important for explaining the data variance (Figure 2). The most important parameter was selected according to its contribution to the principal component, as it is shown by the values of the eigenvectors (in bold) of the correlation matrix (Table A2). The first five principal components accounted for 88.44% of the total variance of the dataset. The linear parameters (hydrology) are the ones mainly explaining the CRB behavior (Table 2).

3.3. Watershed’s Classification Based on Group Analysis

Figure 3 shows the dendrogram, resulting from the group analysis (GA) of the 31 watersheds. Five groups were identified based on their basic, relief, shape, and hydrology-type parameters, and considering a value of R² = 0.84. Group 1, with seven watersheds, has the largest amount of low values for the shape, relief, and linear-type parameters, as can be verified in Table 3. Group 2, with eight watersheds, presents the lowest average values of drainage density, sinuosity of the main channel, and general flow length. Group 3, also with eight watersheds, showed the highest values of elongation ratio, drainage density, mean slope of the main channel, and general flow length. Group 4, with four watersheds, has the highest values in maximum and minimum height, mean slope of the watershed, mean distance, unit shape factor, resistance number, and drainage intensity. In this group, the lowest values correspond to the elongation ratio. Group 5, also with four watersheds, presents the largest amount of high values for the basic, shape, relief, and linear-type parameters. The multivariate analysis of variance (MANOVA) showed significant differences among the groups of watersheds, showing a value of Wilks’ lambda equal to 0.0025, with a value of p < 0.0001.

The geospatial distribution of the groups is shown in Figure 4. Group 1 shows a homogeneous pattern in its distribution, which is concentrated in the central part of the study area. In contrast, Group 2 shows a dispersed distribution, mainly at the edges of the CRB. Group 3 is distributed in the northern, central, and southern parts, and is represented by small clusters of two or three watersheds. Group 4 corresponds to watersheds spatially dispersed over the basin. The watersheds of these groups are isolated from the other watersheds of the same group. Group 5 shows watersheds clustered in the southern part of the basin, except for the watershed RH24Jb, which is located in the northern region. The GA showed a clustered geospatial pattern for Groups 1, 3, and 5, who share characteristics in their parameters and space.

3.4. Watershed’s Classification Based on the Compound Parameter (Cp)

Considering the 26 geomorphometric parameters selected after the correlation analysis was performed, the value of the compound parameter (Cp) was calculated for the 31 watersheds of the CRB (Figure 5a). The watershed RH24Lg received the highest priority (1), followed by the watershed RH24Le (2). The watershed with the lowest priority (31) was watershed RH24Kg. A high priority is an indicator of a high degree of erosion susceptibility in the watershed.

The resulting Cp map (Figure 5a) was reclassified into the following five categories: Very High, High, Moderate, Low, and Very Low (Figure 5b). The spatial distribution of the groups was reclassified by natural breaks [71]. The watersheds classified as Very High show a homogeneous pattern in their distribution in the southwestern part of the watershed. Meanwhile, the watershed RH24Jb is isolated in the northwestern part. The High, Moderate, and Low classes show a dispersed distribution, with at least two of their watersheds clustered in space. The Very Low class shows a homogeneous distribution in space in the southeastern part of the study area, with only one dispersed watershed (RH24Hf).

3.5. Comparison of the Classification Methods

Regarding the GA classification method, the results from the ANOVA analyses performed on 14 geomorphometric parameters, out of 26, detected significant differences among the groups defined by the method. In the case of the Cp classification method, the ANOVA analysis of only two parameters detected significant differences among the groups defined by this classification method (Table 4).

4. Discussion

The prioritization of watersheds, based on susceptibility to erosion, has been carried out in different regions of the world [72,73], using different prioritization methods [74]. This study contributes to the lack of knowledge regarding the susceptibility to erosion in northern Mexico. This was assessed by implementing two methods for prioritization based on the analysis of a set of 33 parameters, which differ from other studies [5,7,75]. The inclusion of several geomorphometric parameters and their relationships within several connected watersheds enriched the study of their erosion susceptibility. In this sense, multivariate techniques have proved to be appropriate methods for establishing priorities, reducing the dimensionality of the dataset by losing the least amount of information [76].

This study integrated a multivariate analysis of several geomorphometric parameters that served to identify those watersheds, which may be prone to erosion. That was possible by evaluating the behavior of such geomorphometric parameters and representing them in a geospatial basis [77,78]. Their relationships provided significant information about the main sources of variability among the studied watersheds [74]. Previous research studies have reported that topography, geomorphology, and land use/land cover are the most important factors in the watershed susceptibility to erosion [79,80,81,82]. In this study, the factors with the greatest influence on the hydrological behavior of watersheds and their erosion susceptibility were the average peak flow and the all channel lengths, as it has also been found in previous studies [6,83].

The PCA is considered a statistical exploratory technique, whose results have helped explain the distribution of environmental attributes [69]. Results from the PCA were useful to identify the sources of variance, which were mainly represented by the dominant parameters influencing the data structure. Then, the basin’s hydrological configuration was explained by those geomorphometric parameters explaining the greatest portion of the variance among the watersheds. The PCA results from this study are consistent with the observations made by Meshram and Sharma [32] and Farhan et al. [33].

From the PCA analysis, PC1 and PC2 are mainly influenced by linear geomorphometric parameters. Some of the linear parameters with an influence on PC1 are the average peak flow. This is shown in Figure 3, where the lowest PC1 coefficients correspond to the watersheds with the lowest mean slope values of the small channels. Regarding PC2, drainage density is one of the linear parameters with an influence. Watersheds with low drainage density indicate the presence of permeable surface material, good vegetation cover, and low relief [84,85]. The map of PC2 (Figure 3) was highly influenced by drainage density, since the watersheds with low values of this parameter are located in the south–central part of the study area and grouped in such a map.

PC3 and PC4 are influenced by linear parameters such as mean distance and shape parameters such as elongation ratio, as well as topographic parameters such as minimum height and mean slope. These factors are associated with the main channel, relief, and slopes, among others. In Figure 3, the watersheds with the greatest heights and slopes correspond to the watersheds located in a mountainous zone, while the watersheds with the lowest elevations and slopes correspond to the arid and semi-arid zones of the state of Chihuahua [86]. Regarding PC5, it is mainly influenced by the elongation index (shape parameter) and the general flow length (linear parameter). The high values of the elongation index correspond to enlarged watersheds, which are related to high drainage densities. A watershed with a high drainage density implies a quick hydrological response to rainfall events, while non-enlarged watersheds correspond to fan-shaped watersheds, which are characterized by short channels [87].

Group analysis (GA) was one of the methodologies used in this study to group and then prioritize the watersheds. It was useful to relate watersheds that share the same characteristics based on their geomorphometric parameters. The groups delineated by the analysis have a unique combination in terms of their geomorphometric attributes [40]. The groups of watersheds follow a territorial pattern. Group 1 includes watersheds located in the zones with the least slopes, where the predominant economic activity is agriculture. Group 2 and 3 belong to watersheds located in transition zones because of their moderate slopes. Groups 4 and 5 are watersheds with rugged topography, with vegetation of shrublands and oak forests.

The compound parameter (Cp) was the second methodology employed to prioritize the watersheds. The value of Cp was calculated for each of the 31 watersheds composing the CRB (Figure 5a). Based on the value of Cp, watershed RH24Hf received the highest priority (1), followed by the watershed RH24Le (2). By comparing the results from GA and Cp, Group 4 was identically integrated by the two methodologies. This group is characterized by watersheds having the highest average values of maximum and minimum height, elongation ratio, elongation index, mean watershed slope, slope of the main channel and unit shape factor. The high values of these parameters correspond to watersheds with a high erosion susceptibility [88]. Conversely, Group 5 was formed by watersheds having the highest values of main channel length, watershed length, contour length, all channel lengths, Gravelius compactness coefficient, shape factor, massivity coefficient, mean distance, Kirpich concentration time, sinuosity of the main channel, average peak flow, texture ratio, river frequency and resistance number. This coincide with high values of Cp, which correspond to watersheds distributed in the southwestern zone of the study area and may have a low erosion susceptibility [61].

The two prioritization schemes used in this study gave similar results according to the spatial distribution of watershed groups. The prioritization of watersheds, obtained through GA and Cp, highlighted those watersheds with potential for the implementation of soil and water conservation practices. Based on the ANOVA analyses performed to statistically compare the GA and Cp methodologies, the former resulted in more effectively classifying the watersheds, since it permitted better differentiating the watershed groups.

Results from the GA show that erosion susceptibility is strongly related to linear parameters (surface hydrology) for southwestern watersheds, where steep slopes of both the watershed and the main channel influence soil erosion [89]. Watersheds RH24Lg, RH24Le, RH24Lf, RH24Mc, RH24Lb, RH24Nc, and RH24Ne have the steepest slopes, making them more prone to erosion [29].

One of the advantages of using the watershed as a territorial unit is the analysis of multiple geomorphometric parameters, which are related to the watershed’s hydrological configuration, topography, and shape. Most of the watershed surface attributes depend on local topographic conditions [5]. In this study, the Basin’s altitudinal gradient, a surface attribute, assists in exhibiting the contrasts among watersheds groups, while showing a homogeneous geographic distribution within them. The linear and shape-type parameters are important because of their influence on soil erosion.

The description and spatial grouping of the 31 watersheds through their 26 parameters using multivariate techniques proved to be useful to understand the main factors that control the variance in the CRB. Prioritization through the two types of grouping was also effective in detecting those watersheds susceptible to erosion. The proposed methodology for prioritizing watersheds on a geospatial basis is a feasible approach for identifying watersheds that are susceptible to erosion. However, prioritization with parameters that are based on shape, linear, and relief of the watersheds may not be sufficient. Thus, the incorporation of information regarding the activities on the territory of the CRB would help to improve the efficacy of the classification of watersheds based on their erosion susceptibility. Therefore, future research could include socioeconomic attributes that contribute to soil loss, such as agriculture [39]. Despite the limitations of this study, the contribution of this work represents an advance in the identification of the watersheds that are most susceptible to erosion in the CRB. This in turn contributes to land-use planning, which may help mitigate soil degradation processes.

5. Conclusions

The application of GA and Cp methodologies allowed integrating a large set of geomorphometric parameters, which served to classify watersheds according to their characteristics.

GA more effectively clustered the watersheds of the Conchos River Basin than Cp, since the groups formed by GA were more differentiated based on the analysis of the watersheds’ geomorphometric parameters. The results of GA show that watersheds RH24Lf, RH24Lb, RH24Nc, and RH24Jb might be subjected to strong erosion processes, and are potential candidates to be subjected to soil conservation practices.

The present study demonstrates the usefulness of integrating GIS and multivariate techniques to prioritize watersheds based on their erosion susceptibility. Such an integration approach showed the spatial relationships of the different geomorphometric parameters analyzed. Although the present study permitted a definition of watershed groups according to the values of their geomorphometric parameters and their relation with erosion susceptibility, the integration of additional variables in the analysis may provide a more insightful classification and thus a more reliable watershed prioritization. Such variables could include land use/land cover, soil type, lithology, geomorphology, and socioeconomic activities, among others.

Author Contributions

Conceptualization, A.P.-A. and F.V.-G.; Data curation, G.V.-Q.; Formal analysis, J.A.P.-A., A.P.-A., M.C.V.-A., A.E.R.-R., and M.M.-S.; Investigation, M.C.V.-A.; Methodology, M.M.-S.; Project administration, A.P.-A.; Software, J.A.P.-A. and G.V.-Q.; Supervision, F.V.-G.; Writing—original draft, J.A.P.-A.; Writing—review and editing, M.C.V.-A., A.E.R.-R., and F.V.-G.

Funding

This research received no external funding.

Acknowledgments

We deeply thank Consejo Nacional de Ciencia y Tecnología (CONACYT-Mexico) for supporting the doctoral graduate studies of the first author of this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Table A1

Values of 33 geomorphometric parameters used for the analysis of the 31 watersheds composing the Conchos River Basin.

Wid	A	P	Lb 2	Lc	Li	Lu	Nu	No 1	Hmin	Hmax	Hmed	Cc	Re	Rf	Ia	J	tgα
RH24Hf	183.19	71.99	12.18	12.66	123.87	505.43	649	336	940	1100	780	3.00	1.25	14.47	0.87	6.76	0.19
RH24Ia	78.62	50.60	24.19	11.17	69.21	189.88	252	111	920	1040	800	3.22	0.41	7.04	1.59	8.80	0.09
RH24Ib	1209.37	259.75	140.20	34.15	2261.86	3915.19	5039	2367	1470	2140	800	4.21	0.28	35.41	0.96	18.70	0.82
RH24Ja	2344.14	297.44	132.14	61.84	3197.22	6430.37	7716	3854	1470	2140	800	3.47	0.41	37.91	1.63	13.64	1.59
RH24Jb	4527.91	373.79	121.30	56.88	5415.42	14,899.51	28,074	13,371	1580	2260	900	3.13	0.63	79.60	0.71	11.96	2.87
RH24Jc	1880.25	290.99	40.60	82.93	1877.94	5304.81	8784	4197	1610	2240	980	3.79	1.21	22.67	3.66	9.99	1.17
RH24Ka	2345.09	283.55	104.99	45.81	3530.78	6993.14	8771	4431	1620	2240	1000	3.30	0.52	51.19	0.89	15.06	1.45
RH24Kb	3447.12	546.62	5.47	106.05	2200.40	7188.10	10,601	5123	1645	2210	1080	5.25	12.12	32.51	3.26	6.38	2.10
RH24Kc	2928.21	433.41	141.09	83.92	3015.00	6069.72	9153	4357	1700	2320	1080	4.52	0.43	34.89	2.41	10.30	1.72
RH24Kd	1076.82	190.50	64.88	37.64	2111.47	2426.44	4431	2100	2020	2640	1400	3.28	0.57	28.61	1.32	19.61	0.53
RH24Ke	380.77	120.27	33.60	22.50	872.80	1023.82	1812	893	2000	2460	1540	3.48	0.66	16.92	1.33	22.92	0.19
RH24Kf	920.06	232.60	93.65	64.15	1030.60	1673.31	2311	1123	1690	2260	1120	4.33	0.37	14.34	4.47	11.20	0.54
RH24Kg	1167.54	193.88	43.29	30.35	717.78	1864.60	2331	1120	1540	1800	1280	3.20	0.89	38.47	0.79	6.15	0.76
RH24La	953.97	169.49	65.07	31.48	781.19	2467.01	3747	1843	1530	1760	1300	3.10	0.54	30.30	1.04	8.19	0.62
RH24Lb	5428.30	640.78	283.10	133.17	10,026.32	13,139.53	19,455	9786	2070	2820	1320	4.91	0.29	40.76	3.27	18.47	2.62
RH24Lc	1760.18	228.02	93.55	51.28	3884.50	2710.82	2872	1428	2320	2820	1820	3.07	0.51	34.32	1.49	22.07	0.76
RH24Ld	2017.67	324.63	16.10	75.21	4868.58	3104.04	3952	1912	2400	2980	1820	4.08	3.15	26.83	2.80	24.13	0.84
RH24Le	2267.52	286.04	145.51	73.49	5722.48	5511.96	12,338	6122	2200	2780	1620	3.39	0.37	30.85	2.38	25.24	1.03
RH24Lf	4296.09	482.75	158.28	97.01	9577.97	9829.22	22,334	11,021	2150	2860	1440	4.16	0.47	44.28	2.19	22.29	2.00
RH24Lg	2618.69	322.00	156.71	85.60	7070.69	6394.05	17,221	8428	2420	3280	1560	3.55	0.37	30.59	2.80	27.00	1.08
RH24Lh	1457.89	220.72	111.39	65.27	946.94	2997.72	5510	2692	1840	2360	1320	3.26	0.39	22.34	2.92	6.50	0.79
RH24Ma	3778.37	465.35	215.32	130.55	1552.45	5018.08	7351	3374	1680	2140	1220	4.27	0.32	28.94	4.51	4.11	2.25
RH24Mb	3808.53	347.58	126.67	92.62	3022.41	8366.91	15,137	7443	1970	2560	1380	3.18	0.55	41.12	2.25	7.94	1.93
RH24Mc	2854.62	435.90	117.91	78.48	5141.29	7178.68	12,804	6264	2330	3020	1640	4.60	0.51	36.37	2.16	18.01	1.23
RH24Md	2485.07	400.59	20.76	122.53	1646.85	3661.01	5738	2822	2070	2880	1260	4.53	2.71	20.28	6.04	6.63	1.20
RH24Me	804.95	211.92	95.40	68.67	1032.98	1770.22	2756	1355	2170	2820	1520	4.21	0.34	11.72	5.86	12.83	0.37
RH24Na	1430.41	232.19	104.82	71.05	574.84	1594.10	2164	1023	1580	2040	1120	3.46	0.41	20.13	3.53	4.02	0.91
RH24Nb	1115.90	186.37	72.02	51.87	1298.61	2417.57	3364	1666	1750	2260	1240	3.15	0.52	21.51	2.41	11.64	0.64
RH24Nc	4655.43	494.37	207.40	110.90	7684.79	10,752.04	23,253	11,600	2020	2780	1260	4.09	0.37	41.98	2.64	16.51	2.30
RH24Nd	2761.21	363.46	166.11	118.17	4289.87	5719.63	8250	4146	2080	2800	1360	3.90	0.36	23.37	5.06	15.54	1.33
RH24Ne	2077.04	315.23	128.60	93.97	3020.71	5385.38	9996	4964	1930	2500	1360	3.90	0.40	22.10	4.25	14.54	1.08
Wid	Co	Dd	j	a	TcK	Scp	Qp	T	Ru	Fu	Rn	Rh	Rr	Lo	Rc	Di
RH24Hf	183.19	2.76	2.63	0.90	3.12	0.96	108.12	9.01	0.94	3.54	3034.96	86.90	15.28	1.38	0.44	1.28
RH24Ia	78.62	2.42	0.99	2.73	7.69	2.17	69.53	4.98	1.26	3.21	2511.63	93.10	20.55	1.21	0.39	1.33
RH24Ib	1209.37	3.24	0.96	4.03	30.21	4.11	289.58	19.40	0.98	4.17	6927.97	62.67	8.24	1.62	0.23	1.29
RH24Ja	2344.14	2.74	1.01	2.73	28.21	2.14	409.07	25.94	1.28	3.29	5870.37	34.61	7.19	1.37	0.33	1.20
RH24Jb	4527.91	3.29	1.12	1.80	25.41	2.13	576.83	75.11	0.85	6.20	7436.74	39.73	6.05	1.65	0.41	1.88
RH24Jc	1880.25	2.82	3.10	0.94	7.39	0.49	364.60	30.19	1.91	4.67	6319.78	27.01	7.70	1.41	0.28	1.66
RH24Ka	2345.09	2.98	1.18	2.17	22.28	2.29	409.16	30.93	0.95	3.74	6679.75	48.89	7.90	1.49	0.37	1.25
RH24Kb	3447.12	2.09	20.67	0.09	0.76	0.05	500.29	19.39	1.81	3.08	4608.40	20.84	4.04	1.04	0.14	1.47
RH24Kc	2928.21	2.07	0.88	2.61	31.35	1.68	459.45	21.12	1.55	3.13	4808.99	27.64	5.35	1.04	0.20	1.51
RH24Kd	1076.82	2.25	1.91	1.98	12.78	1.72	272.55	23.26	1.15	4.11	5948.79	70.15	13.86	1.13	0.37	1.83
RH24Ke	380.77	2.69	2.74	1.72	6.70	1.49	158.41	15.07	1.15	4.76	6614.51	109.32	20.45	1.34	0.33	1.77
RH24Kf	920.06	1.82	1.22	3.09	20.17	1.46	251.06	9.94	2.11	2.51	4110.25	35.23	9.72	0.91	0.21	1.38
RH24Kg	1167.54	1.60	1.20	1.27	11.19	1.43	284.31	12.02	0.89	2.00	2874.65	59.31	9.28	0.80	0.39	1.25
RH24La	953.97	2.59	0.71	2.11	18.79	2.07	255.85	22.11	1.02	3.93	4551.47	55.91	10.38	1.29	0.42	1.52
RH24Lb	5428.30	2.42	0.53	3.84	65.12	2.13	634.11	30.36	1.81	3.58	6825.98	21.18	4.40	1.21	0.17	1.48
RH24Lc	1760.18	1.54	1.07	2.23	21.19	1.82	352.25	12.60	1.22	1.63	4343.02	54.99	12.37	0.77	0.43	1.06
RH24Ld	2017.67	1.54	7.20	0.36	2.62	0.21	378.27	12.17	1.67	1.96	4584.51	39.62	9.18	0.77	0.24	1.27
RH24Le	2267.52	2.43	0.80	3.06	33.33	1.98	402.04	43.13	1.54	5.44	6757.71	37.83	9.72	1.22	0.35	2.24
RH24Lf	4296.09	2.29	0.90	2.41	33.98	1.63	561.22	46.26	1.48	5.20	6543.52	29.48	5.92	1.14	0.23	2.27
RH24Lg	2618.69	2.44	1.10	3.06	31.20	1.83	433.42	53.48	1.67	6.58	8008.77	38.32	10.19	1.22	0.32	2.69
RH24Lh	1457.89	2.06	0.93	2.92	25.53	1.71	319.25	24.96	1.71	3.78	4852.66	36.16	10.69	1.03	0.38	1.84
RH24Ma	3778.37	1.33	0.43	3.50	57.31	1.65	524.84	15.80	2.12	1.95	2842.15	16.39	4.60	0.66	0.22	1.46
RH24Mb	3808.53	2.20	0.93	2.05	28.21	1.37	527.02	43.55	1.50	3.97	5624.03	27.64	7.37	1.10	0.40	1.81
RH24Mc	2854.62	2.51	1.17	2.21	24.45	1.50	453.38	29.37	1.47	4.49	7594.57	38.48	6.93	1.26	0.19	1.78
RH24Md	2485.07	1.47	7.80	0.42	3.09	0.17	421.73	14.32	2.46	2.31	4242.82	23.50	7.19	0.74	0.19	1.57
RH24Me	804.95	2.20	1.36	3.36	19.59	1.39	234.14	13.00	2.42	3.42	6201.64	41.06	13.31	1.10	0.23	1.56
RH24Na	1430.41	1.11	0.88	2.77	24.95	1.48	316.10	9.32	1.88	1.51	2273.45	28.71	8.79	0.56	0.33	1.36
RH24Nb	1115.90	2.17	1.42	2.16	15.54	1.39	277.67	18.05	1.55	3.01	4896.23	43.57	12.13	1.08	0.40	1.39
RH24Nc	4655.43	2.31	0.73	3.04	45.23	1.87	585.26	47.04	1.63	4.99	6420.61	25.07	5.62	1.15	0.24	2.16
RH24Nd	2761.21	2.07	0.87	3.16	35.74	1.41	445.58	22.70	2.25	2.99	5799.98	23.70	7.70	1.04	0.26	1.44
RH24Ne	2077.04	2.59	0.89	2.82	29.09	1.37	384.04	31.71	2.06	4.81	6482.05	26.60	7.93	1.30	0.26	1.86

W_id = Watershed identification, A = Area, P = Perimeter, Lb² = Length of watershed, Lc = Main channel length, Li = Length of contour lines, Lu = Length of channels, Nu = Number of channels, No₁ = Total number of channels of order 1 of the watershed, Hmin = Minimum height, Hmax = Maximum height, Hmed = Medium height, Cc = Gravelius compactness coefficient, Re = Elongation elation, Rf = Form factor, Ia = Elongation index, J = Medium slope of the watershed, tgα = Massive coefficient, Co = Orographic coefficient, Dd = Drainage density, j = Mean slope of main channel, a = Medium distance, TcK = Kirpich concentration time, Scp = Sinuosity of main channel, Qp = Average peak flow, T = Texture ratio, Ru = Unitary shape factor, Fu = River frequency, Rn = Resistance number, Rh = Relief ratio, Rr = Relative relief, Lo = Length of general flow, Rc = Circularity ratio, Di = Drainage intensity.

Table A2

Eigenvectors of the correlation matrix.

GP	PC1	PC2	PC3	PC4	PC5
Lc	0.2678	0.0223	−0.2753	−0.1330	0.0095
Lb 2	0.2345	−0.2721	−0.0327	−0.0059	0.0880
Li	0.2872	0.0343	0.0149	0.1182	−0.1631
Lu	0.3014	0.0344	0.1560	−0.1236	−0.0398
Hmax	0.2170	−0.1034	−0.0543	0.3656	−0.1169
Hmin	0.0664	−0.0941	−0.1183	0.4424	−0.3506
Cg	0.1308	−0.2425	0.0677	0.0028	0.2681
Re	−0.0340	−0.1992	0.4039	0.0072	0.0603
Rf	0.2180	0.1309	0.2002	−0.1898	−0.3142
Ia	0.0510	−0.3141	−0.1533	0.1096	0.3442
J	0.1269	0.1467	−0.0457	0.4034	−0.1466
tgα	0.2759	−0.0721	0.1407	−0.2389	−0.0850
Dd	0.0720	0.3174	0.1766	−0.0248	0.3345
j	−0.0484	−0.2118	0.4035	0.0593	0.0789
a	0.1271	0.1028	−0.4225	−0.0973	0.1850
TcK	0.2483	0.0091	−0.2925	−0.1728	−0.0094
Scp	0.0694	0.2821	−0.2272	−0.1822	0.0258
Qp	0.3031	−0.1104	0.0908	−0.0989	−0.1442
T	0.2642	0.1539	0.1430	0.0413	−0.0230
Ru	0.0641	−0.3193	−0.1582	0.1261	0.3286
Fu	0.1737	0.2478	0.1291	0.1961	0.2488
Rn	0.2305	0.1724	0.0819	0.2424	0.1779
Rh	−0.2207	0.2512	0.0133	0.1126	−0.0453
Rr	−0.2415	0.1661	−0.0986	0.2273	0.0458
Lo	0.0720	0.3174	0.1766	−0.0248	0.3345
Di	0.1911	0.0795	0.0241	0.2992	0.1179

GP = Geomorphometric parameter, PC1 = Principal component 1, PC2 = Principal component 2, PC3 = Principal component 3, PC4 = Principal component 4, PC5 = Principal component 5, Lc = Length of main channel, Lb² = Length of watershed, Li = Length of contour lines, Lu = Length of channels, Hmin = Minimum height, Hmax = Maximum height, Cc = Gravelius compactness coefficient, Re = Elongation ratio, Rf = Form factor, Ia = Elongation index, J = Mean slope of watershed, tgα = Mass coefficient, Dd = Drainage density, j = mean slope of the main channel, a = Medium distance, TcK = Kirpich concentration time, Scp = Sinuosity of the main channel, Qp = Average peak flow, T = Texture ratio, Ru = Unit shape factor, Fu = River frequency, Rn = Resistance number, Rh = Relief ratio, Rr = Relative relief, Lo = General flow length, Di = Drainage intensity. Bold letters indicate the dominant coefficient.

Appendix B

Table A3

Abbreviations.

Abbreviation	Meaning
ANOVA	Analysis of variance
Cp	Compound parameter
CRB	Conchos River Basin
DTM	Digital terrain model
GIS	Geographic information systems
GA	Group analysis
INEGI *	National Institute of Statistics, Geography and Informatics
MANOVA	Multivariate analysis of variance
CONABIO *	National Commission for the Knowledge and Use of Biodiversity
PCA	Principal component analysis
PCs	Principal components

* Acronyms in Spanish.

Figures and Tables

Generate image description

Figure 1. (a) Land-use/land-cover types of the Conchos River Basin (CRB), (b) Delimitation of the 31 watersheds of the CRB, (c) Location of the CRB in Mexico.

Generate image description

Figure 2. Classification of 31 watersheds based on the values of five principal components resulting from a principal component analysis. Conchos River Basin, Chihuahua, Mexico. (a) PC1, (b) PC2, (c) PC3, (d) PC4, (e) PC5. PC1 = Principal component 1, PC2 = Principal component 2, PC3 = Principal component 3, PC4 = Principal component 4, PC5 = Principal component 5.

Generate image description

Figure 3. Dendrogram classifying 31 watersheds by group analysis. Conchos River Basin, Chihuahua, Mexico. The red line was drawn to define the groups.

Generate image description

Figure 4. Geospatial distribution of watershed groups classified by group analysis. Very High (●), High (●), Moderate (●), Low (●), Very Low (●). Conchos River Basin, Chihuahua, Mexico.

Generate image description

Figure 5. (a) Watersheds and their compound parameter (Cp). (b) Geospatial distribution of watershed groups by Cp reclassification. Very High 1 (●), High 2 (●), Moderate 3 (●), Low 4 (●), Very Low 5 (●). Conchos River Basin, Chihuahua, Mexico.

Table 1

Geomorphometric parameters.

Geomorphometric Parameter	Equation	References
Basic Parameters 1
Area (A)	A = Watershed surface area [km2]	Horton [15]
Perimeter (P)	P = Watershed perimeter [km]	Horton [15]
Length (Lb2)	Lb2 = Watershed length [km]
Stream order (u)	u = Stream order [unitless]	Strahler [16]
Main Channel Length (Lc)	Lc = Main flow channel length [km]
All Channel Lengths (Lu)	Lu = Length of all the flow channels in the watershed [km]	Horton [15]
Contour Length (Li)	Li = Contour lines’ length [km]
Number of Flow Channels (Nu)	Nu = Number of flow channels [unitless]
Number of First-Order Flow Channels (No1)	No1 = Number of total first-order flow channels in the watershed [unitless]
Maximum Height (Hmax)	Hmax = Watershed maximum height [m]
Minimum Height (Hmin)	Hmin = Watershed minimum height [m]
Medium Height (Hmed)	Hmed = Watershed medium height [m]
Shape Parameters
Gravelius Compactness Coefficient (Cc)	Cc = P/2√πA	Zavoianu [24]
Elongation Ratio (Re)	Re = 1.1284 (√A/Lc)	Schumm [25]
Shape Factor (Rf)	Rf = A/Lb2	Horton [15]
Elongation Index (Ia)	Ia = Lb2/Wwhere: W = watershed width (Km)	Horton [15]
Unit Shape Factor (RU)	RU = Lb2/A0.5	Horton [15]
Circularity Ratio (Rc)	Rc = 4πA/P2	Miller [56]
Relief Parameters
Mean Watershed Slope (J)	J = (ΣLi E/A) × 100where: E = equidistance among contour lines (km)	Horton [15]
Massivity Coefficient (tgα)	tgα = Hmed/A
Relief Relationship(Rh)	Rh = Hmax/Lb	Schumm [25]
Relative Relief (Rr)	Rr = Hmax/P	Schumm [25]
Orographic Coefficient (Co)	Co = Hmed × tgα
Linear parameter
Drainage Density (Dd)	Dd = ΣLu/A	Horton [15]
Mean Slope of the Main Channel (j)	j = (Hmax − Hmin)/Lc × 100	Horton [15]
Mean Distance (Am)	Am = Lc/(√A)	Horton [15]
Sinuosity of the Main Flow Channel (Scp)	Scp = Lc/Lb2	Mueller [57]
Kirpich Concentration Time (TcK)	TcK = 0.066(Lb2/j)0.77	Kirpich [58]
Average Peak Flow (Qp)	Qp = 43A0.522	Sen [59]
Texture Ratio (T)	T = Nu/P	Horton [15]
Rivers Frequency (Fu)	Fu = Nu/A	Horton [15]
Resistance Number (Rn)	Rn = Hmax × Dd	Schumm [25]
General Flow Length (Lo)	Lo = 1/2 × Dd	Schumm [25]
Drainage Intensity (Di)	Di = Fu/Dd	Faniran [60]

¹ The basic parameters were calculated with the GIS software.

Table 2

Principal components and geomorphometric parameters explaining the greatest portion of variance in the watersheds of the Conchos River Basin, Chihuahua, Mexico.

PC	Eigenvalue	Variance	Accumulated Variance	Geomorphometric Parameter (First)	Geomorphometric Parameter (Second)
1	9.0010	0.3462	0.3462	Average peak flow (Qp)	Channel lengths (Lu)
2	6.0848	0.2340	0.5802	Unit shape factor (Ru)	Drainage density (Dd)
3	3.4222	0.1316	0.7119	Mean distance (Am)	Elongation ratio (Re)
4	2.8620	0.1101	0.8219	Minimum height (Hmin)	Mean slope (J)
5	1.6246	0.0625	0.8844	Elongation index (Ia)	General flow length (Lo)

PC = Principal component, First = Primary importance, Second = Secondary importance.

Table 3

Average values of the geomorphometric parameters by group.

Gid	Lb 2	Lc	Li	Lu	Hmin	Hmax	Cc	Re	Rf	Ia	J	tgα	Dd
1	46.12	31.31	371.43	1038.50	1495	995	3.22	0.74	20.03	1.69	6.43	0.49	1.97
2	80.16	65.56	1820.63	2660.91	2489	1415	3.70	0.91	23.28	3.11	13.62	0.79	1.97
3	102.47	75.16	2765.79	6206.70	2294	1060	3.95	1.99	34.73	2.41	12.07	1.48	2.59
4	146.56	88.94	5556.08	6201.08	2970	1545	3.86	0.40	30.30	3.10	21.45	1.17	2.36
5	192.52	99.49	8176.13	12,155.08	2680	1230	4.07	0.44	51.66	2.20	17.31	2.45	2.58
Gid	j	a	TcK	Scp	Qp	T	Ru	Fu	Rn	Rh	Rr	Lo	Di
1	1.42	1.92	11.74	1.51	194.51	8.83	1.24	2.56	2673.67	67.01	13.48	0.99	1.30
2	2.44	2.17	18.48	1.37	313.28	16.48	1.69	3.03	4835.28	47.81	11.26	0.98	1.51
3	3.70	2.18	22.19	1.69	417.90	27.78	1.50	3.86	5915.17	34.49	6.97	1.30	1.51
4	0.98	2.87	31.18	1.68	433.61	37.17	1.73	4.87	7040.26	34.58	8.63	1.18	2.04
5	0.82	2.77	42.44	1.94	589.36	49.69	1.44	4.99	6806.71	28.86	5.50	1.29	1.95

G_id = Group identification, 1 = Group 1 (Very low erosion susceptibility), 2 = Group 2 (Low erosion susceptibility), 3 = Group 3 (Moderate erosion susceptibility), 4 = Group 4 (High erosion susceptibility), 5 = Group 5 (Very high erosion susceptibility), Lc = Length of main channel, Lb² = Length of watershed, Li = Length of contour lines, Lu = Length of channels, Hmin = Minimum height, Hmax = Maximum height, Cc = Gravelius compactness coefficient, Re = Elongation elation, Rf = Form factor, Ia = Elongation index, J = Mean slope of watershed, tgα = Mass coefficient, Dd = Drainage density, j = Mean slope of main channel, a = Medium distance, TcK = Kirpich concentration time, Scp = Sinuosity of main channel, Qp = Average peak flow, T = Texture ratio, Ru = Unitary shape factor, Fu = River frequency, Rn = Resistance number, Rh = Relief ratio, Rr = Relative relief, Lo = Length of general flow, Di = Drainage intensity.

Table 4

p-values resulting from the ANOVA analyses performed to detect differences among the groups defined by each classification method.

	Geomorphometric Parameters
	Lc	Lb 2	Li	Lu	Hmax	Hmin	Cc	Re	Rf	Ia	J	tgα	Dd
GA	0.0091	0.0865	<0.0001	<0.0001	0.1954	0.3315	0.4011	0.2675	0.0005	0.4273	0.0750	<0.0001	0.0224
Cp	0.3711	0.3576	0.3207	0.1462	0.7602	0.1029	0.4368	0.4983	0.3214	0.8336	0.5236	0.1225	0.1563
	j	a	TcK	Scp	Qp	T	Ru	Fu	Rn	Rh	Rr	Lo	Di
GA	0.3316	0.7712	0.0356	0.5116	<.0001	<.0001	0.5327	0.0069	0.0004	0.1205	0.0086	0.0224	0.0265
Cp	0.5902	0.6454	0.4377	0.2158	0.2287	0.1442	0.8277	0.0174	0.1626	0.4913	0.1979	0.1563	0.0302

GA = Group analysis, Cp = Compound parameter, Lc = Length of main channel, Lb² = Length of watershed, Li = Length of contour lines, Lu = Length of channels, Hmin = Minimum height, Hmax = Maximum height, Cc = Gravelius compactness coefficient, Re = Elongation elation, Rf = Form factor, Ia = Elongation index, J = Mean slope of watershed, tgα = Mass coefficient, Dd = Drainage density, j = mean slope of main channel, a = medium distance, TcK = Kirpich concentration time, Scp = Sinuosity of main channel, Qp = Average peak flow, T = Texture ratio, Ru = Unitary shape factor, Fu = River frequency, Rn = Resistance number, Rh = Relief ratio, Rr = Relative relief, Lo = Length of general flow, Di = Drainage intensity.