A Multivariate Geomorphometric Approach to Prioritize Erosion-Prone Watersheds

: 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 runo ﬀ source. However, the CRB is subjected to signiﬁcant 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 classiﬁed 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 deﬁned linear-type parameters as the main source of variability among the watersheds. The GA and the Cp were e ﬀ ective for grouping the watersheds in ﬁve groups, and provided the information for the spatial analysis. The GA methodology best classiﬁed 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. 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 e ﬃ cacy compared.


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 basins according to their erosion susceptibility. For the grouping, multivariate techniques and the compound parameter (Cp) were used and their efficacy compared.

Study Area
The CRB is located in the state of Chihuahua and Durango, Mexico, and is part of the 24 th Rio Bravo-Conchos Hydrological Region [47] (Figure 1). The basin has an area of 67,800 km 2 [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].

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;

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 2 ), 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 Sustainability 2019, 11, 5140 5 of 21 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].

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): 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 2 ) [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]: 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].

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).

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 2 and 50.59 km, respectively. Meanwhile, the watershed with the largest area and perimeter is RH24Lb, with 5428.29 km 2 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.

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).  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 2 = 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.

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).

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).

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.

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.