Flash-Flood Susceptibility Assessment Using Multi-Criteria Decision Making and Machine Learning Supported by Remote Sensing and GIS Techniques

: Concerning the signiﬁcant increase in the negative e ﬀ ects of ﬂash-ﬂoods worldwide, the main goal of this research is to evaluate the power of the Analytical Hierarchy Process (AHP), ﬁ (kNN), K-Star (KS) algorithms and their ensembles in ﬂash-ﬂood susceptibility mapping. To train the two stand-alone models and their ensembles, for the ﬁrst stage, the areas a ﬀ ected in the past by torrential phenomena are identiﬁed using remote sensing techniques. Approximately 70% of these areas are used as a training data set along with 10 ﬂash-ﬂood predictors. It should be remarked that the remote sensing techniques play a crucial role in obtaining eight out of 10 ﬂash-ﬂood conditioning factors. The predictive capability of predictors is evaluated through the Information Gain Ratio (IGR) method. As expected, the slope angle results in the factor with the highest predictive capability. The application of the AHP model implies the construction of ten pair-wise comparison matrices for calculating the normalized weights of each ﬂash-ﬂood predictor. The computed weights are used as input data in kNN–AHP and KS–AHP ensemble models for calculating the Flash-Flood Potential Index (FFPI). The FFPI also is determined through kNN and KS stand-alone models. The performance of the models is evaluated using statistical metrics (i.e., sensitivity, speciﬁcity and accuracy) while the validation of the results is done by constructing the Receiver Operating Characteristics (ROC) Curve and Area Under Curve (AUC) values and by calculating the density of torrential pixels within FFPI classes. Overall, the best performance is obtained by the kNN–AHP ensemble model.


Study Area
The location of the upper and middle parts of the Prahova river basin is in the south-central part of Romania. The study area is defined by the parallels of 45 •  The study area is characterized by the elevations between 128 m and 2505 m on the highest mountain peaks. The study area is spread on both the Carpathian region and the Subcarpathian region and it is characterized by a mean slope angle of 12.43 • [37]. The lithology of the study area contains hard rocks, flysch and sedimentary rocks such as sand, gravel and loess. The amount of precipitation in the Subcarpathian region is approximately 600 mm/year, while in the Carpathian region it increased to 1200 mm/year. During the summer, heavy rainfalls frequently generate severe flash-flood events. Land-use also substantially influences the surface runoff genesis. Thus, the built-up surfaces that create favorable conditions for flash-flood occurrence are spread approximately on 10% of the selected basin area. The afforestation represents 48.2%. Regarding the soil characteristics, the hydrological soil group A covers 801 km 2 , hydrological soil group B has an area of 727 km 2 and the hydrological soil group D covers 682 km 2 . It should be noted that the hydrological soil group D, which has the highest potential in the genesis of surface runoff, covers around 25% of the selected basin area.

Data
Given the fact that this study aims at the mapping of areas that are highly susceptible to flashflooding, only the geospatial data were used for the analysis. Particularly, it should be mentioned that the RS and GIS were applied to process the input dataset used in the present research.

Inventory of Torrential Areas
The inventory of torrential areas is mandatory to correctly assess the relationship between the selected influencing factors and flash-flood susceptibility [37]. Based on the work of Costache and Zaharia [38], the surfaces, such as ravines and gullies, were used for analysis, representing the actual torrential microforms. The remote sensing techniques were involved in this stage of the study and the aerial imagery from the Google Earth application was used to delineate the areas affected by torrential processes within the study area ( Figure 1). The validation of these areas was performed by field surveys using a professional GNSS device. Finally, a total surface of 260 km 2 , representing torrential areas and containing about 289,000 pixels, was delineated ( Figure 1). This surface represents approximately 10% of the entire territory covered by the study area (2600 km 2 ).

Flash-Flood Conditioning Factors
To compute the FFPI for the study area, ten flash-flood variables were selected. Seven morphometric factors based on the digital elevation model (DEM) and the other three factors were obtained from vector databases. The DEM for the study area was generated from the Shuttle Radar Topography Mission (SRTM), which has a spatial resolution of 30 m, using RS techniques [39]. The DEM included in the SRTM database was obtained based on measurements carried out through the Synthetic Aperture Radar (SAR) interferometry [40]. Based on the DEM, the following variables were derived: slope angle, aspect, plan curvature, profile curvature, convergence index, topographic position index (TPI) and topographic wetness index (TWI). Thus, the application of RS techniques

Data
Given the fact that this study aims at the mapping of areas that are highly susceptible to flash-flooding, only the geospatial data were used for the analysis. Particularly, it should be mentioned that the RS and GIS were applied to process the input dataset used in the present research.

Inventory of Torrential Areas
The inventory of torrential areas is mandatory to correctly assess the relationship between the selected influencing factors and flash-flood susceptibility [37]. Based on the work of Costache and Zaharia [38], the surfaces, such as ravines and gullies, were used for analysis, representing the actual torrential microforms. The remote sensing techniques were involved in this stage of the study and the aerial imagery from the Google Earth application was used to delineate the areas affected by torrential processes within the study area ( Figure 1). The validation of these areas was performed by field surveys using a professional GNSS device. Finally, a total surface of 260 km 2 , representing torrential areas and containing about 289,000 pixels, was delineated ( Figure 1). This surface represents approximately 10% of the entire territory covered by the study area (2600 km 2 ).

Flash-Flood Conditioning Factors
To compute the FFPI for the study area, ten flash-flood variables were selected. Seven morphometric factors based on the digital elevation model (DEM) and the other three factors were obtained from vector databases. The DEM for the study area was generated from the Shuttle Radar Topography Mission (SRTM), which has a spatial resolution of 30 m, using RS techniques [39]. The DEM included in the SRTM database was obtained based on measurements carried out through the Synthetic Aperture Radar (SAR) interferometry [40]. Based on the DEM, the following variables were derived: slope angle, aspect, plan curvature, profile curvature, convergence index, topographic position index (TPI) and topographic wetness index (TWI). Thus, the application of RS techniques was crucial for obtaining the accurate DEM and, consequently, seven flash-flood conditioning factors. Moreover, the RS techniques, Remote Sens. 2020, 12, 106 5 of 26 namely the supervised classification method, had a vital contribution in the delimitation of land-use categories. These land-use categories were further used to obtain the Curve Number, which is one of the most critical flash-flood predictors. The flash-flood predictors, along with their influence on runoff generation, are described in the following text.
Lithology has an important influence on surface runoff genesis because this factor controls the water infiltration rate. The lithological classes for the research area were extracted from the digital version of the Geological Map of Romania (1:200,000). Around 40% of the study area is built with marls, flysch and conglomerates. The lithological map was subsequently converted to 30 m raster (Figure 2a). Regarding the convergence index (Ci), its importance for the computation of flash-flood potential was highlighted in many of the previous studies [38,[41][42][43]. The river valleys are highlighted by the negative values of this morphometric parameter, while the positive ones indicate the interfluvial zones. The Ci values were classified into the following five classes: (−96)-(−3), (−3)-(−2), (−2)-(−1), (−1)-0 and 0-95 ( Figure 2b). Profile curvature differentiates the areas prone to runoff (negative values) from those (positive values) which are less susceptible to surface runoff [43]. More than 50% of the torrential pixels belong to the class between −8.2 and 0 ( Figure 2c). Curve Number (CN) (Figure 2d) is widely used to estimate the direct runoff based on a given rainfall event. Concerning a specific territory, the CN values, which may range between 0 to 100, are derived based on the land-use categories and hydrological soil groups (HSG) [37,[44][45][46]. The HSG control the infiltration rate, the soil moisture and surface runoff, while the land-use has a high influence on surface runoff due to the different Manning coefficients [47,48]. The HSG were derived from the vector format of the Soil Map of Romania (1:200,000). Moreover, the land-use/land-cover was derived from the CORINE Land Cover (2018) database. It should be mentioned that the CORINE Land Cover dataset was obtained by the supervised classification of Sentinel-2 and Landsat 8 images. Using the natural breaks (Jenks) classification, the resulting CN values were classified as follows: 0-43, 43-52, 52-69, 69-80, 80-89, 88-100 ( Figure 3d). The values between 43 and 80 cover 70% of the total basin area and approximately correspond to 45% of the torrential areas. Plan curvature is a flash-flood predictor which highlights the surfaces which are characterized by convergent or divergent runoff. Thus, this morphometric factor is very useful in analyzing the areas susceptible to surface runoff. The highest share (58%) of the study area corresponds to the plan curvature values from 0.1 to 0.5 ( Figure 2e). Modified Fournier Index (MFI), as another flash-flood conditioning factor, describes the spatial variability of rainfall intensity within a specific region [37,49]. The following Equation (1) can be used to calculate the MFI values [50]: where P i is the monthly average amount of precipitation for month i in mm, and P is the average annual rainfall. The average amount of precipitation from 23 meteorological stations located inside and around the study area was used ( Figure 1). Based on the work of Kourgialas and Karatzas [50], the spatial variability of the precipitation was modeled with the use of Spline interpolation. Finally, the MFI values were derived by using the precipitation amount in raster format and Equation (1). These values were reclassified into five classes using natural breaks (Jenks) algorithm: 53. 56-57.43, 57.44-62.83, 62.84-68.69, 68.7-75.14, 75.15-83.46 (Figure 2f). The class between 53.76 and 57.43 quantifies around 28% of the study area. The highest percentage of torrential pixels (36.6%) is included in the MFI class between 68.7 and 75.14. According to the vast majority of studies [48,[51][52][53][54][55][56], slope degree has the highest impact on the generation of surface runoff. This morphometric flash-flood predictor was derived from the DEM (Figure 3a). The following classes were established for slope angle: 0-3 • , 3-7 • , 7-15 • , 15-25 • and >25 • [49]. Although the slope class between 15 • and 25 • contains almost 40% of the torrential areas, the highest density of torrential pixels was recorded in the class higher than 25 • . Thus, the slope class higher than 25 • , which quantifies 7.08% of the total study area, includes 20.7% of the total surface of torrential areas ( Figure 4). This situation can be explained by the high velocity of water runoff that is recorded on slopes higher than 25 • . The torrential relief microforms also are present in a high percentage (32.8%) on slopes between 7 • and 15 • because the gravitational force is high enough to generate rapid surface runoff.
Remote Sens. 2020, 12, 106 6 of 28 of the torrential areas, the highest density of torrential pixels was recorded in the class higher than 25°. Thus, the slope class higher than 25°, which quantifies 7.08% of the total study area, includes 20.7% of the total surface of torrential areas ( Figure 4). This situation can be explained by the high velocity of water runoff that is recorded on slopes higher than 25°. The torrential relief microforms also are present in a high percentage (32.8%) on slopes between 7° and 15° because the gravitational force is high enough to generate rapid surface runoff. The values of TPI denote the elevation difference between a grid cell and the neighboring cells [57]. The resulting TPI values were classified as follows: (-34.82)-(-2.08); (-2.08)-(-0.84); (-0.83)-0.14; 0.15-1.12; 1.13-27.96. The highest share (45%) is represented by the TPI class from (-0.83) to 0.14 ( Figure 3b). TWI is an important runoff predictor which was computed based on Equation (2) [58]:  The values of TPI denote the elevation difference between a grid cell and the neighboring cells [57]. The resulting TPI values were classified as follows: (−34.82)-(−2.08); (−2.08)-(−0.84); (−0.83)-0.14; 0.15-1.12; 1.13-27.96. The highest share (45%) is represented by the TPI class from (−0.83) to 0.14 ( Figure 3b). TWI is an important runoff predictor which was computed based on Equation (2) [58]: where α is the cumulative upslope area draining through a point (per unit contour length), and tan β is the slope angle at the point. The classification of the TWI values is as follows: −9.7-4.5; 4.6-8.4; 8.5-12; 13-15; 16-25 ( Figure 3c). The largest areas correspond to the medium TWI class (8.5-12) which has a total share of 29% ( Figure 4).
Aspect is a flash-flood predictor, derived from the DEM, which values were divided into ten classes. Slope aspect has an indirect role on the surface runoff due to its influence on the other factors such as rainfall regime, soil humidity, and solar radiation. Having a percentage of 15%, eastern and south-western slopes occupy the most extensive surfaces (Figure 3d). Eastern surfaces contain approximately 20% of the total torrential pixels ( Figure 4).
Aspect is a flash-flood predictor, derived from the DEM, which values were divided into ten classes. Slope aspect has an indirect role on the surface runoff due to its influence on the other factors such as rainfall regime, soil humidity, and solar radiation. Having a percentage of 15%, eastern and south-western slopes occupy the most extensive surfaces (Figure 3d). Eastern surfaces contain approximately 20% of the total torrential pixels ( Figure 4).

Analytical Hierarchy Process (AHP)
The AHP belongs to multi-criteria decision-making methods that solve complex problems in a simpler way [59]. The AHP model is based on the active participation of decision-makers within the entire methodological workflow [60]. Many studies applied this method for the evaluation of susceptibility to other natural hazards [26,36,52,61,62]. By using this method, an unstructured problem is broken into many components. The workflow applied in the AHP model can be described within the following five main steps [63]: (i) Defining objectives and dividing the unstructured problem into its components.
(ii) Determining the detailed criteria and alternatives.
(iii) Creation of the pair-wise comparison matrix, which is constructed based upon the expert opinion regarding the influence of each factor or factor class/category to flash-flood occurrence. The relative value of a factor can range from 1 to 9 when that factor is more important than another and conversely from 1/2 to 1/9 when that factor is less important than another [64].
(iv) Computing the relative weights of each criterion (flash-flood predictors) using the eigenvalue technique.
where λ is the largest eigenvalue of the matrix which can be determined from the comparison matrix, and n is the number of flash-flood predictors or the number of factor classes/categories.
where RI is the random consistency index which relates to the number of factors included in the comparison matrix [65]. When the CR value is less than 0.1, then it can be stated that the comparison is consistent [65].

k-Nearest Neighbor (kNN)
According to Marjanovic et al. [66], this simple machine learning algorithm classifies a point in n-dimensional input space, taking into account the value of class containing the k-closest neighboring elements using the training dataset. Therefore, this method is used especially in predictive analysis [67]. Basically, in the kNN algorithm, the elements of the same geographical site will have the same characteristics if they are located near each other [68]. This algorithm uses a simple voting system to attribute, to a spatial object, a new class value that is most common in the neighboring instances [66]. The new class value is assigned according to a distance metric that exists between the new spatial object and the k-nearest neighbor. Usually, to define this metric, the Euclidean distance is calculated as follows (Equation (5)) [69]: where x i is the coordinate for x point and z i is the coordinate for z point.
Other distance measures such as Manhattan, Chebyshev, and Hamming can be more suitable for other model settings.
Regarding the present study, the kNN algorithm estimated the conditional probability for each class of flash-flood predictor (X) to belong to the torrential class (y = 1) or non-torrential class (y = 0). The conditional probability can be determined using Equation (6) [70]: where i is the torrential or non-torrential pixel within training dataset, T, k is the number of torrential or non-torrential pixels, which are located in the proximity of each class of flash-flood predictors.
The k value represents a positive integer which is equal to the number of neighbors used to calculate the new class value for a specific spatial object ( Figure 5). Therefore, the determination process of an optimal k value is critical due to influences on a considerable measure of the results of the kNN model and the accuracy of the flash-flood susceptibility map. A low value of k may lead to a large prediction variance while a large value of k may cause a large model bias.
Remote Sens. 2020, 12, 106 10 of 28 where is the torrential or non-torrential pixel within training dataset, , is the number of torrential or non-torrential pixels, which are located in the proximity of each class of flash-flood predictors.
The value represents a positive integer which is equal to the number of neighbors used to calculate the new class value for a specific spatial object ( Figure 5). Therefore, the determination process of an optimal value is critical due to influences on a considerable measure of the results of the kNN model and the accuracy of the flash-flood susceptibility map. A low value of may lead to a large prediction variance while a large value of may cause a large model bias.

Lazy K-Star (KS) Algorithm
To the best of our knowledge, the presented approach represents the first attempt to determine the susceptibility to a specific natural hazard through the use of a lazy K-Star (KS ) algorithm. The KS algorithm is an instance-based learning algorithm which creates explicitly the shortest string that connects two instances using the Kolmogorov distance [71] . The total probability of all paths from instance a to an instance b represents the probability function used in a KS algorithm and is defined as follows (Equation (7)): where P * is the probability function, t represents the value of T, which is the set of transformations predefined. Thus, the KS function is written as follows (Equation (8)): Apart from the fact that it is not considered exactly a distance function, the KS algorithm also has the following properties highlighted in Equation (9) and Equation (10) [72]: The implementation of the KS model requires, at the beginning of the training process, a correct estimation of optimal parameters for X0 (real numbers) and s (symbolic attributes) [73]. Concerning each type of parameter, the probability distribution includes a number of instances varying from 1 to N. The number of instances is equal to N in that case, if all the instances are equally weighted [73].

Lazy K-Star (KS) Algorithm
To the best of our knowledge, the presented approach represents the first attempt to determine the susceptibility to a specific natural hazard through the use of a lazy K-Star (KS) algorithm. The KS algorithm is an instance-based learning algorithm which creates explicitly the shortest string that connects two instances using the Kolmogorov distance [71]. The total probability of all paths from instance a to an instance b represents the probability function used in a KS algorithm and is defined as follows (Equation (7)): where P * is the probability function, t represents the value of T, which is the set of transformations predefined. Thus, the KS function is written as follows (Equation (8)): Apart from the fact that it is not considered exactly a distance function, the KS algorithm also has the following properties highlighted in Equations (9) and (10) [72]: The implementation of the KS model requires, at the beginning of the training process, a correct estimation of optimal parameters for X 0 (real numbers) and s (symbolic attributes) [73]. Concerning each type of parameter, the probability distribution includes a number of instances varying from 1 to N. The number of instances is equal to N in that case, if all the instances are equally weighted [73]. Different functions of conditional probability (P * ) can be used to estimate this number based on Equation (11) [71]: where N is the total number of training instances and n 0 is the smallest distances from a.
The main advantage of a KS algorithm is represented by its high effectiveness when it works with large data sets and, also, by the fact that this model is robust to noisy training data [74].
To compute the FFPI values using the KS algorithm, Weka 3.9 software [75] was used.

Proposed Methodology for Predicting Flash-Flood Potential
The entire methodological workflow is graphically illustrated in Figure 6. As it can be observed, the flash-flood susceptibility was computed through two stand-alone and two ensemble methods.
Remote Sens. 2020, 12, 106 11 of 28 Different functions of conditional probability (P * ) can be used to estimate this number based on Equation (11) [71]: where N is the total number of training instances and n is the smallest distances from a.
The main advantage of a KS algorithm is represented by its high effectiveness when it works with large data sets and, also, by the fact that this model is robust to noisy training data [74].
To compute the FFPI values using the KS algorithm, Weka 3.9 software [75] was used.

Proposed Methodology for Predicting Flash-Flood Potential
The entire methodological workflow is graphically illustrated in Figure 6. As it can be observed, the flash-flood susceptibility was computed through two stand-alone and two ensemble methods.

Establishment of Flash-Flood Database
During the first step of the methodological workflow, the flash-flood database for the study area was established applying the ArcGIS 10.5 software. The database that was used in this study includes the spatial extension of the 10 flash-flood predictors described in Sub-Section 3.2 and a total surface of 260 km 2 (around 289,000 pixels) containing areas with torrential phenomena. According to literature [16,32,35], to obtain a higher performance of the machine learning models, another sample with non-torrential areas was generated for the training process, which has the same number of pixels

Establishment of Flash-Flood Database
During the first step of the methodological workflow, the flash-flood database for the study area was established applying the ArcGIS 10.5 software. The database that was used in this study includes the spatial extension of the 10 flash-flood predictors described in Section 3.2 and a total surface of 260 km 2 (around 289,000 pixels) containing areas with torrential phenomena. According to literature [16,32,35], to obtain a higher performance of the machine learning models, another sample with non-torrential areas was generated for the training process, which has the same number of pixels (289,000) as in the case of the torrential regions. These surfaces were randomly selected from territories with a slope angle below 3 • . The slopes lower than 3 • were chosen due to their very low potential for surface runoff; these surfaces were generally not influenced by torrential phenomena. Therefore, it was likely that flash-flooding would occur on these surfaces.
It should be mentioned that all three aforementioned elements of the flash-flood database were processed into 30 m resolution rasters.

Selection of Flash-Flood Predictors Applying Information Gain Ratio (IGR) Method
The uncorrelated flash-flood predictors, i.e., predictors having very low predictive capability, with flash-flood phenomena may generate noisy input data, which may result in decreased predictive capability of the applied models [76]. Concerning this, the IGR method was employed to evaluate the predictive capability of the selected flash-flood predictors [77]. Particularly, a high value of IGR demonstrates that the flash-flood predictor has a high predictive capability.
IGR for each flash-flood predictor was computed with the use of Equations (12) where D is the training dataset composed of n input samples, n (Y i , D) is the number of samples in the training data D belonging to the class label Y i (torrential, non-torrential). Weka 3.9 software was employed to estimate the predictive ability of each flash-flood predictor.

Computing the AHP Weights for Factor Classes/Categories
The computation of the AHP weights for each factor class/category was carried out after the construction of 10 pair-wise comparison matrices, each one corresponding to each flash-flood predictor. Therefore, by assigning a relative dominant value, each factor and class/category in the same flash-flood predictor was rated against every other. The entire workflow and all the matrices therein used Microsoft Excel 2016 according to the steps described in Section 4.1. The quality of the comparisons was tested by using the consistency ratio calculated according to Equation (4). The AHP weights were used as input in the kNN and KS models.

Training and Validation Dataset Preparation
Concerning the majority of the studies approaching the assessment of susceptibility to natural hazards using machine learning algorithms, the initial dataset consists of training and validating pixels. According to literature [34,[78][79][80], a percentage of 70% and 30% of the dataset was used to train and validate the models, respectively. Regarding the present case, both the torrential and non-torrential samples were divided according to the aforementioned percentages. Thus, the training dataset contains 202,300 torrential pixels (≈182 km 2 ) or 70% and 202,300 non-torrential pixels (≈182 km 2 ), also 70%.
The validating dataset includes 86,700 torrential pixels (≈78 km 2 ) or 30%, and 86,700 non-torrential pixels (≈78 km 2 ), also 30%. The Subset Features tool of ArcGIS 10.5 software was employed to establish the training and validation datasets.
Furthermore, each pixel from the training and validation samples was attributed the normalized values (between 0 and 1) of the initial numerical factors and the categories of the categorical factors. According to literature [81], the machine learning models require that the input data are normalized using the same range, since the bias may occur in the results due to the bigger magnitude of the initial untransformed data. The normalization was performed in terms of the work by Amiri et al. [82]. This dataset was applied for training the kNN and KS stand-alone models. Regarding terms of kNN-AHP and KS-AHP ensembles, the pixels from the training and validation samples were assigned the AHP weights calculated according to the methodology described in Sections 4.1 and 5.3.

Configuration and Training of Flash-Flood Potential Models
The entire training process of the kNN stand-alone and kNN-AHP ensemble was performed with the help of XLSTAT Microsoft Excel 2016 software. Regarding this, the training and validation data in GIS format was converted into tabular format to be imported into Microsoft Excel 2016. Further, an essential step in the computation of FFPI through the kNN stand-alone model and kNN-AHP ensemble was the determination of the optimal number of Nearest Neighbors or k-number. Thus, a 10-fold cross-validation procedure was employed, based on the work of Nguyen et al. [83], to find the best k-value. Therefore, the best k-number was estimated by modifying the k-values until the highest classification accuracy, calculated for both kNN and kNN-AHP, was reached. The best k-values were 19 for kNN and 18 for kNN-AHP.
Concerning the training of the KS stand-alone and KS-AHP ensemble, Weka 3.9 software was used. Regarding this, in the first step, the same data used for kNN models were converted into Comma Separated Values (CSV) format to be imported into the Weka 3.9 software. Using the Weka software, the performance of the KS models is dependent on the quality of training data, but it is highly influenced also by the selection of a global blending parameter. Similar to the procedure applied by Naji et al. [84], the best value of a global blending parameter for both the KS and KS-AHP was established after a trial and error procedure, taking into account the classification accuracy. Concerning terms of the KS stand-alone model, the best global blending parameter was established at 16 while for the KS-AHP ensemble model, the best global blending value was 20.

Evaluation of the Model Performance
After configuring and running the models, the quality of the results was evaluated through several statistical indices like Sensitivity, Specificity, and Accuracy. The description and detailed definitions of these indices were included in several studies regarding the flash-flood or landslide susceptibility [4,85]. The statistical measures used in this study were calculated through the following Equations (16)- (18): Speci f icity = TN FP + TN (17) where FP (false positive) and FN (false negative) are the number of pixels erroneously classified.

Flash-Flood Potential Mapping and Results Validation
Once the best configuration for the four models was established, the algorithms were run, and the flash-flood predictor importance was determined. These importance values were finally used to derive the FFPI. The importance flash-flood predictors determined through the kNN and KS stand-alone models were multiplied with the rasters without AHP weights, while the importance derived through the application of the kNN-AHP and KS-AHP ensembles was multiplied with the rasters having assigned the AHP weights. Regarding these operations, the ArcGIS 10.5 software was applied.
During the first stage, the validation of flash-flood potential mapping was made by calculating the share of training and validating samples in the computed FFPI classes. This is a widely used method for results validation [38], by which the total percentage of training and validation pixels that correspond to the high and very high FFPI class can be quantified.
During the second stage, the results validation was done through the Receiver Operating Characteristic (ROC) Curve, which belongs to widely used methods for the assessment of the model accuracy [86,87]. The acquired results and existing torrential areas were compared when creating the ROC Curve. Therefore, both the success and prediction rates were created to assess the reliability and accuracy of the FFPI maps. Regarding this, the training 70% of the torrential areas were used to create the success rate, while the validating 30% of torrential areas were used to construct the prediction rate. The Area Under Curve (AUC) in the case of the success rate follows how the model classified the results. Considering the case of the prediction rate, it highlights the accuracy of the results. The model is defined as efficient when the value of AUC is close to 1. When the AUC value is close to 0, it highlights a non-informative model [88]. The AUC value can be estimated through the following Equation (19) [3]: where TP (true positive) and TN (true negative) are the numbers of pixels that are correctly classified, P is the total number of pixels with torrential phenomena and N is the total number of pixels without torrential phenomena.

Predictive Ability of Flash-Flood Conditioning Factors
The results of the IGR method, which was applied to determine the predictive ability of the selected flash-flood conditioning factors, are shown in Figure 7.
Remote Sens. 2020, 12, 106 14 of 28 During the second stage, the results validation was done through the Receiver Operating Characteristic (ROC) Curve, which belongs to widely used methods for the assessment of the model accuracy [86,87]. The acquired results and existing torrential areas were compared when creating the ROC Curve. Therefore, both the success and prediction rates were created to assess the reliability and accuracy of the FFPI maps. Regarding this, the training 70% of the torrential areas were used to create the success rate, while the validating 30% of torrential areas were used to construct the prediction rate. The Area Under Curve (AUC) in the case of the success rate follows how the model classified the results. Considering the case of the prediction rate, it highlights the accuracy of the results. The model is defined as efficient when the value of AUC is close to 1. When the AUC value is close to 0, it highlights a non-informative model [88]. The AUC value can be estimated through the following Equation (19) [3]: where (true positive) and (true negative) are the numbers of pixels that are correctly classified, P is the total number of pixels with torrential phenomena and N is the total number of pixels without torrential phenomena.

Predictive Ability of Flash-Flood Conditioning Factors
The results of the IGR method, which was applied to determine the predictive ability of the selected flash-flood conditioning factors, are shown in Figure 7.   Based on the results of the analysis, slope angle was the most significant flash-flood conditioning factor with an average predictive ability of 0.9, followed by curve number (0.87) profile curvature (0.79), Modified Fournier Index (0.69), lithology (0.63), plan curvature (0.58), TWI (0.51), convergence index (0.46), TPI (0.39) and aspect (0.29). It can be observed that the lowest value of average merit was 0.29, which signifies that all the factors are important in a flash-flood occurrence process and, therefore, were taken into account for the presented analysis.

AHP Weights Results
Table S1 from the Supplementary Materials showed all possible pair-wise comparisons between the classes of the selected flash-flood predictors. The weight values that indicate the importance of each class/category were also included in Table 1. Among all class/categories, the highest importance, equal to 0.633, was achieved by the profile curvature class between 0.9 and 9.7, followed by the plan curvature class between 0.1 and 0.5 (0.512), slope angles higher than 25 • (0.507) and convergence index class between −96 and −3 (0.485).
Mentioned in Section 4.2, the evaluation of the consistency of expert judgments was assessed through the Consistency Ratio (CR) values, which were computed for each comparison matrix. According to Table 1, all CR values are less than 0.1, which demonstrates that all the comparisons between the class/categories within the same predictor are consistent. Table 2 also highlights the values of several parameters used to determine the weights of each factor or the CR value.

Application of kNN and kNN-AHP Ensemble Model
Based on the cross-validation procedure, the optimal k-value was estimated to be used for training purposes of the kNN individual model and kNN-AHP ensemble. Regarding the kNN stand-alone model, it can be seen that the highest accuracy (82.  (Table 2).
Stated in Section 5.7, the importance of each factor, in terms of Flash-Flood Potential Index-k-Nearest Neighbor(FFPI kNN ) and Flash-Food Potential Index-k-Nearest Neighbot-Analytical Hierarchy Process (FFPI kNN-AHP ), were derived. Regarding terms of FFPI kNN , the slope angle result was the factor with the highest importance, having a weight of 20.9%. The rest of the factors resulted in the following weights: curve number (19.8%), lithology (12.2%), MFI (10.9), convergence index (9.9%), TPI (8.2%), TWI (7.4%), plan curvature (6.1%), profile curvature (2.9%) and aspect (1.7%). Furthermore, the importance was included in map algebra to calculate the FFPI kNN . The FFPI kNN values (0.04-0.96) (Figure 9a), were reclassified into five classes based on the natural breaks (Jenks) method. This grading method is considered the most appropriate for the classification of values into the classes [89]. Very low flash-flood potential values ranged between 0.04 and 0.26 and their share of the study area is approximately 11.3%. According to Figure 9a, these surfaces can be found mostly along the main river valleys. The second class has a share of approximately 33.9% and relates to the surfaces with low flash-flood potential. The medium FFPI kNN values are spread homogenously across the study area and are included in the range 0.43-0.56. Concurrently, the middle class of flash-flood potential covers 31% and spreads across the whole study area. Altogether, the fourth and fifth FFPI kNN classes have a share of 23.8% (Figure 10 (Table 2). Stated in Sub-Section 5.7, the importance of each factor, in terms of Flash-Flood Potential Indexk-Nearest Neighbor (FFPIkNN) and Flash-Food Potential Index -k-Nearest Neighbot-Analytical Hierarchy Process (FFPIkNN-AHP), were derived. Regarding terms of FFPIkNN, the slope angle result was the factor with the highest importance, having a weight of 20.9%. The rest of the factors resulted in the following weights: curve number (19.8%), lithology (12.2%), MFI (10.9), convergence index (9.9%), TPI (8.2%),  Concerning terms of the kNN-AHP ensemble, the highest weight was reached again by the slope angle (21.4%). The weights of the rest of the factors are as follows: curve number (19.3%), lithology (13.2%), MFI (11.2%), convergence index (8.1%), TPI (7.4%), TWI (6.6%), plan curvature (5.7%), profile curvature (4.3%) and aspect (2.8%). Similar to FFPI AHP and FFPI kNN , the values of flash-flood potential, which this time vary between 0.03 and 0.35, were classified into five classes based on the natural breaks (Jenks) method (Figure 9b). The values in the class between 0.03 and 0.14 relate to the areas with a very low potential for flash-flooding. A percentage of 7.42% of the study area is covered by surfaces where the flash-floods are the least probable. Low FFPI kNN-AHP values recorded a share of around 26.74%, while the medium values had the share of approximately 34.98%. FFPI kNN-AHP values higher than 0.21 belong to areas which have high to very high flash-flood potential. These areas recorded a share of 30.8% and can be found, especially, in the northern half of the study area.
Concerning terms of the kNN-AHP ensemble, the highest weight was reached again by the slope angle (21.4%). The weights of the rest of the factors are as follows: curve number (19.3%), lithology (13.2%), MFI (11.2%), convergence index (8.1%), TPI (7.4%), TWI (6.6%), plan curvature (5.7%), profile curvature (4.3%) and aspect (2.8%). Similar to FFPIAHP and FFPIkNN, the values of flashflood potential, which this time vary between 0.03 and 0.35, were classified into five classes based on the natural breaks (Jenks) method (Figure 9b). The values in the class between 0.03 and 0.14 relate to the areas with a very low potential for flash-flooding. A percentage of 7.42% of the study area is covered by surfaces where the flash-floods are the least probable. Low FFPIkNN-AHP values recorded a share of around 26.74%, while the medium values had the share of approximately 34.98%. FFPIkNN-AHP values higher than 0.21 belong to areas which have high to very high flash-flood potential. These areas recorded a share of 30.8% and can be found, especially, in the northern half of the study area.

Application of KS and KS-AHP Ensemble Model
According to Table 2

Application of KS and KS-AHP Ensemble Model
According to Table 2 FFPIKS-AHP was spatially modelled by assigning the following weights for each flash-flood conditioning factor: slope angle was 20.8%, curve number was 19.2%, lithology was 16.5%, MFI was 13.3%, convergence index was 11.4%, TPI was 5.9%, TWI was 5.1%, plan curvature was 3.5%, profile curvature was 2.7% and aspect was 1.6%. Therefore, the values in the range between 0.05 and 0.53 were reclassified into five classes based on the natural breaks (Jenks) method. The interval between 0.05 and 0.14 covered an area equal to 6.58% of the studied river basin and corresponded to the very FFPI KS-AHP was spatially modelled by assigning the following weights for each flash-flood conditioning factor: slope angle was 20.8%, curve number was 19.2%, lithology was 16.5%, MFI was 13.3%, convergence index was 11.4%, TPI was 5.9%, TWI was 5.1%, plan curvature was 3.5%, profile curvature was 2.7% and aspect was 1.6%. Therefore, the values in the range between 0.05 and 0.53 were reclassified into five classes based on the natural breaks (Jenks) method. The interval between 0.05 and 0.14 covered an area equal to 6.58% of the studied river basin and corresponded to the very   Table 3). Table 3. Share of torrential pixels in Flash-Flood Potential Index (FFPI) classes.   Table 3).

Conclusions
The accurate identification of highly susceptible areas is the basis for adopting the main non-structural measures intended to prevent the negative impacts of flash-floods. Taking this context, the presented article proposed a complex methodology for identifying the areas susceptible to flash-flooding in the upper and middle parts of the Prahova river basin. The methodology adopted was based on the computation of FFPI using two machine learning models (k-Nearest Neighbor (kNN) and KS (K-star)) and their novel ensemble with an Analytical Hierarchy Process (AHP). The prediction power of the ensemble models, namely kNN-AHP and KS-AHP, also was tested. The weights determined by the pair-wise comparison matrices were used as input data for kNN-AHP and KS-AHP ensembles. Furthermore, ten flash-flood predictors were selected through the IGR method, along with 70% of torrential areas. These predictors, as well as 70% of torrential area, were used to train the models. The validation of the results was performed with the use of the other 30% of torrential areas. The assessment of the model performances was done by the computation of several statistical measures. Regarding the validation of the results, it was carried out in two steps: i) by the share of torrential areas in the FFPI classes; ii) through the Receiver Operating Characteristics ROC Curve. Concerning the ROC Curve, the best model was the kNN-AHP ensemble for both success rate (AUC = 0.901) and prediction rate (AUC = 0.896). It should be remarked that when the kNN was used in combination with other models for determining the susceptibility to different natural hazards, it also obtained very good performances. An example in this regard was provided by Bui et al. [90] who calculated the landslide susceptibility in Vietnam through a combination of a fuzzy kNN algorithm with differential evolution optimization. The AUC of the success rate in that case was 0.944, while for the prediction rate it was 0.841.
Areas which resulted in high to very high flash-flood potential had a share between 23.8% (FFPI kNN ) and 36.12% (FFPI KS-AHP ). These surfaces were found, especially, in the north-western part characterized by the presence of tourist localities, which creates the premise of a high vulnerability degree to flash-flood occurrence. These localities were frequently affected by floods and flash-floods in the following years: 2005, 2006, 2010, 2014, 2017 and 2018 [91].
The main novelty of the presented research is considered to be the use of the KS model, which, to the best of our knowledge, is the first attempt to use it for the assessment of susceptibility to natural hazards, as well as the application of the two ensembles for estimating the FFPI. Even if the kNN-AHP method achieved the highest accuracy, the KS was used for FFPI computation due to its advantages presented in Section 4. These advantages are highlighted by the higher performances obtained by the KS stand-alone model in comparison with the kNN stand-alone model. Nevertheless, it was observed that, when the input data was represented by reclassified predictors, having assigned the AHP coefficients, the performance of the kNN-AHP exceeded the performance of the KS-AHP.
This study provides results which are applicable for effective spatial planning as well as for improving the flash-flood forecast and warning activities conducted by the National Institute of Hydrology and Water Management of Romania.