Integrated Framework for Detecting the Areas Prone to Flooding Generated by Flash-Floods in Small River Catchments

: In the present study, the susceptibility to flash-floods and flooding was studied across the Izvorul Dorului River basin in Romania. In the first phase, three ensemble models were used to determine the susceptibility to flash-floods. These models were generated by a combination of three statistical bivariate methods, namely frequency ratio (FR), weights of evidence (WOE), and statistical index (SI), with fuzzy analytical hierarchy process (FAHP). The result obtained from the application of the FAHP-WOE model had the best performance highlighted by an Area Under Curve— Receiver Operating Characteristics Curve (AUC-ROC) value of 0.837 for the training sample and another of 0.79 for the validation sample. Furthermore, the results offered by FAHP-WOE were weighted on the river network level using the flow accumulation method, through which the valleys with a medium, high, and very high torrential susceptibility were identified. Based on these valleys’ locations, the susceptibility to floods was estimated. Thus, in the first stage, a buffer zone of 200 m was delimited around the identified valleys along which the floods could occur. Once the buffer zone was established, ten flood conditioning factors were used to determine the flood susceptibility through the analytical hierarchy process model. Approximately 25% of the total delimited area had a high and very high flood susceptibility.

vulnerability and risk to flash floods. It should be noted that most of the procedures regarding the evaluation of flash-flood and flood risk assessment, which were adopted by the European countries includes the use of several traditional methods such as hydraulic and hydrological modeling. These techniques are time consuming and very expensive [5,6]. In this context, the need to find faster, more accurate, and cheaper techniques for determining the flood hazard has significantly increased.
In recent years, the scientific field of flash-food and flood susceptibility assessment has had a high dynamic due to the fast development of the techniques and software used to perform these analyses [7]. Thus, to assess the flood susceptibility, Geographic Information System (GIS) techniques, complex models of bivariate statistics, and machine learning are used [8]. The most used bivariate statistical techniques for assessing susceptibility to natural hazards are weights of evidence [9], frequency ratio [10], evidential belief function [11], certainty factor [12], statistical index [13], and index of entropy [14]. The most well-known machine learning models used in the study of susceptibility to floods are decision trees [15], multilayer perceptron [16], logistic regression [17], support vector machine [18], bagging [19], k-nearest neighbor [20], naïve Bayes [21], Decorate [22], Dagging [15], and adaptive neuro-fuzzy inference system [23]. Many researchers have assessed the risk of flash-floods and floods by using ensemble models resulting from the combination of several methods [15,21,24].
Nevertheless, in all the research papers where machine learning and bivariate statistics were used, the susceptibility was estimated separately for flash-floods and flooding. Up to now, there is no approach in which the susceptibility to these two phenomena, which are strongly related, can be estimated together. A first attempt to identify the torrential valleys, based on the flash-flood susceptibility, was done by Costache et al. [25]. In that study, the authors managed to detect the river valleys prone to flash-flood propagation using four hybrid models and the flow accumulation method. Nevertheless, the flooding susceptibility was not estimated along the torrential river valleys, this fact being a shortcoming that should be addressed.
In this context, we aimed to propose an integrated approach for estimating the surface runoff susceptibility and the susceptibility to floods. Thus, in the first stage, we follow the identification of areas susceptible to flash-floods by applying three overall models generated by combining frequency ratio, statistical index, and weights of evidence bivariate statistics models, on the one hand, and fuzzy analytical hierarchy process on the other hand. The models' performances were evaluated utilizing the ROC curve. The second stage of the study aims to identify the torrential valleys susceptible to the propagation of the upstream flash-floods by applying the flow accumulation method. Once the valleys with a medium, high, and very high potential for flash-flood propagation are identified, the flood susceptibility is calculated to determine the areas exposed to floods. Flood susceptibility is determined through the analytical hierarchy process stand-alone model.
It should be mentioned that this is the first time in the literature when the susceptibility of these two phenomena, flash-floods and flooding generated by them, were analyzed in an integrated way and in a spatial causal relationship. The previous studies carried out in Romania as well as in any part of the globe were focused on the estimation of flooding or flash-flood susceptibility without taking into account their strong spatial relationship.

Study Area
The present study focused on the Izvorul Dorului River basin located in the mountainous area of the central part of Romania. The surface of the study area is 33 km 2 , which falls into the category of small-area basins that are frequently affected by rapid floods. The altitude inside the study zone varies from 763 m to 2202 m (Figure 1a). This high difference in altitude on a small area creates favorable premises for flash-flood genesis and their propagation from the upper to the lower part of the river basin. The river basin is characterized by an average high slope of 15.6°, which is another indicator of the high potential for flash-flood propagation along the valleys in the study area. According to the existing information and, as can be seen in Figure 1b, the afforestation degree of the river basin is around 50%. Additionally, from Figure 1b, one can remark that in the perimeter of the deforested surfaces located at the highest altitude also exists a very high potential for a rapid surface runoff on the slopes. This is another element indicating that the genesis of the flash-floods is related to the high-altitude region of the river basin from where they are propagated along the steep river valleys toward the lowland area. The lower part of the study area corresponds to the built space of Sinaia city, the most famous mountain tourist resort in Romania. This locality has been affected by floods multiple times, caused by Izvorul Dorului River and its tributaries. One of the most violent flash-floods took place in August 2010 when several dozens of buildings were affected as well as National Road 1, National Road 71, and the railroad between Bucharest and Brasov cities. Additionally, as a result of different strong floods, several landslides were activated and affected the houses from Sinaia.

Statistical Index
Proposed by van Westen [26], statistical index (SI) is a bivariate method widely used in natural risk susceptibility evaluation studies [13,27]. According to this model, the score of a predictor class can be computed by applying the natural algorithm to the ratio between the density of pixels associated with the phenomenon presence in the predictor class and the density of the same pixels across the study area [28]. Thus, a well-known formula to estimate the SI weight is the following: where Wij is the weight of the class/category i of predictor j; fij is the density of the phenomenon in class i of predictor j; f is the density of phenomenon in the study; Npix(Si) is the number of pixels associated with the phenomenon in class i; and Npix(Ni) is the sum of pixels of the same parameter class.

Frequency Ratio
Frequency ratio (FR) is a bivariate statistical model, widely applied to evaluate flood and landslide susceptibility mapping worldwide [9,10,13]. The relationship between food occurrences and conditioning parameters is used to analyze and calculate the frequency ratio. The mathematical expression of frequency ratio (i.e., the frequency ratio of class i of factor j) is given in Equation (2) [10]: where Npix (1) is the total number of torrential points contained by a class/category of factor; Npix (2) is the total number of pixels contained by each class/category; ∑ (3) is the total number of torrential pixels within the study area; and ∑ (4) is the total number of pixels within the study area.
After calculating the frequency ratio, each controlling factor summed up all the values to generate a map of flood vulnerability. If the frequency ratio is greater than 1, the conditioning factors strongly influence flooding, otherwise, there is a negative relationship between conditioning factors and flood occurrence.

Weights of Evidence
Weights of evidence (WOE) is a widely used statistical model for landslide, flood, and fire forest susceptibility assessment [29][30][31]. This method was first introduced for geological studies in 1992, then adopted for the analysis of different hazards (e.g., fire forest, flood, landslides) [27]. This method estimates the weights of evidence coefficients based on the relationship between each class of factors and the flood absence/presence. The positive weight (W + ) and the negative weight (W − ) are necessary for the computation. These weights reflect the presence and absence of areas affected by the flood, respectively, and can be computed using the following [29][30][31]: where B and are the presence and absence of flood conditioning parameters, respectively; P is the probability; and S, and ̅ are the presence and absence of flooding, respectively.
The output of the performed processes is used to implement Equations (3) and (4) in ArcGIS. Subsequently, the mathematical representation of these two equations are [29]: where W + and W − are the positive and negative weights, respectively; Npix 1 and Npix 2 are the number of pixels with flood points inside and outside of the class, respectively; and Npix 3 and Npix 4 are the number of pixels without flooding inside and outside of the class, respectively.
The final weights of evidence coefficients (Wf) assigned to each factor class can be obtained as follows [29]: where (Wf) is the final weight of evidence coefficients; Wmintotal is the total of all negative weights in a multiclass map; and Wplus and Wmin are the positive negative weights of a class factor, respectively.

Fuzzy Analytical Hierarchy Process
The analytical hierarchy process (AHP) is an algorithm used for flood, landslide, and fire forest susceptibility mapping [32][33][34][35]. Through a pairwise comparison matrix constructed based on expert knowledge, AHP was used to calculate the weights of relevant criterion map layers. Since AHP has several advantages such as its fuzzy extension, the fuzzy analytical hierarchy process (FAHP) was proposed and applied to solve the hierarchical fuzzy problems. It can be employed to increase the analysis quality, reducing the subjectivity in the estimation of weights criteria by a combination of the fuzzy set theory and the analytical hierarchy process [36]. The following steps show how to determine the weights of criteria in the FAHP.
The pairwise comparison matrices are constructed from flood conditioning factors (elevation, slope angle, stream density, curve number, rainfall, lithology, land use, soil texture, etc.). Linguistic terms are assigned to the pairwise comparison (Equation (8)) to establish the most important criteria [37]: where a'ij indicates a pair of criteria i and j. The Buckley method [38] is utilized to calculate the fuzzy geometric mean and fuzzy weight of each criterion by: where is the fuzzy comparison value between the pair criterion i and criterion n; and is the geometric mean of the fuzzy comparison values for criterion i compared to each of the other criteria; is the fuzzy weighting of the ith criterion; and = ( , , ), where , and are the values of the lower, middle, and upper, fuzzy weighting of the ith criterion, respectively [37,39].
The extent analysis algorithm was applied to determine the final values of the flood conditioning factor weights. The construction of a fuzzy triangular comparison matrix is the first step. This matrix is done by [40]: where = ( , , ) and = (1/ , 1/ , 1/ ) for i, j = 1, …, n and i ≠ j. Next, we computed the priority vector of the triangular matrix. Then, the fuzzy arithmetic function was employed to sum up each row of the matrix in a first stage, as follows: Then, the value of the fuzzy synthetic extent in terms of the ith object is obtained through the normalization of the above relation, as follows [32]: The computation of the degree of possibility of ≥ represents the third step and is achieved through Equation (14): where = ( , , ) and = ( , , ). Considering that: the weight vector values can be calculated by: The weight vectors were obtained using the following equation after a normalization process: where w is a non-fuzzy number. The present study was carried out by completing several methodological steps, as presented in Figure 5 and also briefly described below.

Torrential Areas Inventory
Identifying the areas previously affected by a natural risk phenomenon is vital for detecting other zones where that phenomenon has a high probability of occurrence [41]. The appearance of any phenomenon will be favored in areas with characteristics similar to those where the phenomenon has already occurred [42]. For this reason, to estimate the susceptibility to the occurrence of rapid floods, torrential areas were inventoried and mapped. These areas were generated by the rapid surface runoff on the slopes. The modality of identification of such zones is presented in the studies [43]. Torrential areas are zones characterized by the unified presence of a torrential microform of relief such as ravines and gullies generated by surface runoff. Thus, through the satellite images made available through the Google Earth application (Figure 1), an area affected by intense torrential processes of about 170 hectares was vectorized, which is located in the upper part of the river basin where the absence of vegetation and the high slopes favor the apparition of such phenomena.

Flash-Flood and Flood Predictors
Whereas torrential zones represent an indicator of the rapid surface runoff on the slopes, certain geographical factors are the predictors of this phenomenon, or in other words, are the variables that generate and favor the surface runoff. Moreover, the genesis of floods generated by flash-floods also depends on the characteristics of geographical factors. Therefore, to identify as accurately as possible the surfaces favorable to the genesis of flash-floods and those susceptible to floods, twelve conditioning factors were taken into account. Eight morphometrical predictors were obtained by processing the digital elevation model, while the other four flood and flash-flood predictors were collected from the following vector databases: hydrological soil groups from the Digital Soil Map of Romania, 1:200,000; land use/cover from Corine Land Cover, 2018; lithology from the Digital Geological Map of Romania, 1:200,000; and distance from rivers was estimated with the help of the river network in an Environmental Systems Research Institute (ESRI) shapefile format. Below, the main characteristics of flood and flash-flood predictors are briefly presented.
The slope is the geographic factor that has the biggest influence on both the potential for rapid surface runoff and the flood potential [24]. Surfaces with steep slopes cause rapid water drainage, while the flat surfaces lead to the water accumulation process [44]. In our case study, the sloping relief had values between 0.1° and 54.1° ( Figure 2a). This interval was divided into six classes according to the literature [43].
Land use/cover is another predictor that influences both flash-floods and floods [45]. Lands covered with pastures or without vegetation will favor the appearance of rapid runoff on the slopes, while areas covered with forests are characterized by a lower potential for runoff and flooding [21]. In the study area, the grassland and the forest shared equally almost all of the territory ( Figure 2b). Additionally, the presence of the built space in the lower part of the Izvorul Dorului River basin was observed.
Hydrological soil group has a high influence on the flood. Thus, the flooding phenomenon will likely be over the areas with soils with high clay content such as those in hydrological group D, while water infiltration will be more pronounced on soils with a sandy texture [46][47][48]. Within the study area, the largest surface was occupied by hydrological soil group A (Figure 2c).
Convergence index (CI) is a predictor obtained from the DEM whose values show the concentration degree of the drainage network. CI values close to −100 indicate a high density of the river network whereas positive CI values are associated with the interfluvial surfaces. In the study area, the CI values are situated in the range from −86 to 84 ( Figure  2d). These were divided into five classes according to the literature [43].
Profile curvature is a predictor whose negative values show the surfaces that favor the accelerated surface runoff, while the decelerated runoff manifests itself on the surfaces with positive values. The information from the literature was used to classify profile curvature values into the next classes: −2.3-−0.1; 0-0.1; 0.2-2.6 ( Figure 3a).
The aspect factor obtained from the DEM is an indicator of the humidity potential that exists at the slope level [49]. In the case of the Izvorul Dorului basin, the southeast surfaces were the most extensive, these being followed by the southwest slopes ( Figure  3b).  The elevation is a useful indicator for detecting the surfaces exposed to flooding processes that may occur as a result of flash-flood propagation from the upper part of river basins [7]. The lower relief zones have a higher sensitivity to flooding occurrence. For the study area, the range from 763.1 m to 2202 m was split into seven classes that generally succeeded at a difference of 200 m (Figure 4a). Lithology is a predictor that influences the infiltration capacity at the ground surface, so it should be considered in the studies concerning the flood and flash-flood potential. The conglomerates, breccias, sandy flysch, and marls shale are predominant in the study area (Figure 4c).
Distance from the river was generated using the Euclidean distance tool from ArcGIS 10.3 software. This is an important parameter that indicates the distance from different surfaces to the nearest watercourse. The surfaces in the vicinity of watercourses will be more prone to flash-floods and the floods generated by them. In the present study, the distance from the river predictor was classified into eight classes.

Step 1: Flash-Flood Database Preparation
The flash-flood database used in the present research consisted of 1965 torrential points collected from the delineated torrential surfaces and ten flash-flood conditioning factors. Building and processing the flash-flood database were done through ArcGIS 10.3 software. It should be noted that the torrential points were obtained by converting the torrential areas from a raster format, with a cell size of 30 m, to a point. Therefore, each point corresponds to a raster cell. According to the literature [51], the entire sample was divided into a training dataset (70%) and a validation dataset (30%). The training dataset was used to calculate the frequency ratio, weights of evidence, and statistical index coefficients, while the validation dataset was used to evaluate the accuracy of the results achieved.

Step 2: Computation of Flash-Flood Potential Index (FFPI)
The flash-flood potential index represents a qualitative indicator of the potential for torrential surface runoff, which exists at the slope level [52]. In the first stage, the frequency ratio, weights of evidence, and statistical index coefficients were determined by analyzing the spatial correlations between the torrential points included in the training sample and the ten flash-flood predictors. In this regard, the equations from Sections 3.1.1-3.1.3 were implemented in Excel and ArcGIS. The number of pixels used in the computation of the types of bivariate statistics coefficients was 1376. Furthermore, the second stage consisted of the computation of flash-flood predictors weights by the fuzzy analytical hierarchy process method. Finally, three variants of the flash-flood potential index were computed by the weighted sum between fuzzy analytical hierarchy process weights and the values of frequency ratio, weights of evidence, and statistical index coefficients.

Step 3: Evaluation of Results Accuracy Using Receiver Operating Characteristic (ROC) Curve
The results of FFPI were assessed using the receiver operating characteristic (ROC) curve. The ROC curve represents a graphical plot that highlights the ability of a binary model to classify a given dataset used in the modeling process into the presence or the absence of a specific phenomenon [53]. This is the most frequently used algorithm to validate the outcomes provided by a model for natural hazards susceptibility [42,49,54,55]. The ROC curves were constructed by comparing the existing torrential points with the flash-flood potential index results. Both the success rate, constructed with the training sample, and prediction rate constructed with the validation sample, will be used. The area under curve (AUC) will highlight the performance of each flash-flood potential index model.

Step 4: Computation the Flood Potential Index (FPI) Based on the Most Performant FFPI Result
To identify the valleys with a high torrential degree, the best performing flash-flood potential index that resulted was used in a flow accumulation procedure ( Figure 5). Through the flow accumulation method, the flash-flood potential index values are weighted at the level of the river network within the study area. The weighted flash-flood potential index values are further classified into five categories: very low, low, medium, high, and very high. In the next stage, to select the river valleys along which the flood potential index will be calculated, the hydrographic network having assigned a medium, high, and very high flash-flood propagation susceptibility is selected. The flood potential index represents a qualitative indicator that highlights the degree to which a specific region can be affected by the flooding phenomenon [56]. The area on which the flood potential index will be computed was limited to a buffer zone of 200 m along with the selected river network. Eventually, the flood potential index values are obtained by involving the next ten flood conditioning factors in the analytical hierarchy process method: slope, land use, hydrological soil groups, convergence index, topographic position index, topographic wetness index, elevation, distance from the river, plan curvature, and lithology. The values of the FPI are then classified into five categories through which the areas prone to flooding generated by flash-floods will be detected.

Bivariate Statistics Coefficients
The values of bivariate statistics coefficients highlight the spatial relationships between the location of torrential areas and the classes/categories of flash-flood predictors.

Flash-Flood Potential Index Computation Using Fuzzy Analytical Hierarchy Process Based Ensembles
Following the methodological steps described in Section 3.1.4 the fuzzy analytical hierarchy process algorithm was applied to determine the weights of the flash-flood predictors. In the first step, the fuzzy analytical hierarchy process evaluation matrix was created based on expert judgment (Table 2). Furthermore, using the values included in the evaluation matrix, the synthesis values were calculated by using the formula:  The synthesis values calculated above were used in the following step to calculate the fuzzy numbers for each flash-flood predictor. The fuzzy numbers were then compared using the degree of possibility procedure, which is exemplified in Table 3. Utilizing the results provided by the degree of possibility method, the weight vector values were calculated using the following relations: In the next step, employing the defuzzification procedure, the Triangular Fuzzy Numbers (TFNs) were transformed into the crisp weights that will be attributed to each flash-flood predictor and multiplied with statistical index, frequency ratio, and weights of evidence values to obtain the flash-flood potential index. Table 3. The ordinate of the highest intersection point, the degree possibility for Triangular Fuzzy Numbers (TFNs), and the weights of the flash-flood predictors.

Flash-Flood Potential Index Results Validation
Results validation is a mandatory step to establish the best ensemble model whose results will be used to identify the areas prone to flood generated by flash-floods. In this regard, the success rate and prediction rate were used. The success rate revealed that the highest performance was obtained by the results provided by FAHP-WOE (AUC = 0.837), followed by FAHP-SI (AUC = 0.833) and FAHP-FR (AUC = 0.723) (Figure 7a). The same hierarchy was also revealed by the construction of the prediction rate. Thus, the FAHP-WOE ranked first (AUC = 0.79), followed by FAHP-SI (AUC = 0.787) and FAHP-FR (AUC = 0.717) (Figure 7b). Therefore, following the results validation procedure, the FFPIFAHP-WOE was selected to be used in the next step of the analysis.

Flood Potential Index Computation
According to the description provided in Section 3.3.4, the flow accumulation method was applied to FFPIFAHP-WOE to evaluate the torrential degree of the river valleys across the study area. The results showed that a percentage of 21.59% of the total river valleys identified were characterized by a low and very low torrential degree and are, therefore, considered to be less favorable for flash-flood propagation (Figure 8a). For a 200 m buffer zone along with the other 78.41% of the river valleys, the flood potential index (FPI) was calculated. In this regard, the stand-alone analytical hierarchy process (AHP) multicriteria decision-making was used. It should be mentioned that through AHP, in the first stage, the weights of flash-flood predictors and classes/categories of flash-flood predictors were calculated. In terms of flash-flood predictors, the highest weight was detected for slope (0.224), followed by land use (0.137), elevation (0.137), distance from river (0.137), lithology (0.085), plan curvature (0.081), hydrological soil groups (0.064), convergence index (0.055), TPI (0.046), and TWI (0.031) ( Table 4). The analysis of the weights attributed to the classes/categories of flash-flood predictors revealed that the highest value was obtained for hydrological soil group D (0.66), followed by the plan curvature class between −3 and −0.1 (0.539), the TPI class between −7.8 and −1.8 (0.439), the TWI class between −4.4 and 4.7 (0.433), the conglomerates and breccias lithological categories (0.423), the convergence index class between −86 and −3 (0.42), and the slope angle class lower than 3° (0.379). The consistency of judgments was evaluated using the consistency ratio (CR) values. The results from Table 5 show that the CR values were less than 0.1, indicating that all the comparisons within the matrices were consistent. Table 5 also contains the values of some parameters involved in the CR computation. To derive the flood potential index across the study area, the AHP weights, together with the raster dataset associated with the flood predictors, were used in map algebra of ArcGIS software. The normalized values of FPI were classified into five classes using the natural break method. The very low class, between 0 and 0.12, covered about 18.9% of the total area and was mainly spread along the valleys located in the southern part of the catchment. Another 23% of the delimited zone was characterized by a low flood potential. The medium FPI values (between 0.32 and 0.53) were associated with about 33.3% of the delimited perimeter. The high and very high potential was spread around 24.8% and was associated with FPI values higher than 0.53 (Figure 8b).

Discussion
In the last ten years, the study of the susceptibility to hydrological risk phenomena has developed significantly as a result of the combined application of geospatial analysis techniques with statistical models or algorithms from artificial intelligence [49]. It is well known that small river basins located in mountainous areas favor the occurrence of flashfloods and their propagation toward the areas located in the lower zones of the basins. The mountainous area of Romania is not an exception and is often affected by severe flashflood events that generate property damage and loss of life. In this context, the present study aimed to identify the areas exposed to floods generated by flash-floods within the Izvorul Dorului River basin located in the Romanian Carpathians, which could produce pollution such as the transport of polycyclic aromatic hydrocarbons resulting from different sources [57].
The present study included a first part in which the potential for rapid water runoff on the slopes was determined, the second part referred to the identification of valleys with high torrential potential, followed by the evaluation of flood susceptibility existing along these valleys. The potential for rapid surface runoff, expressed through the FFPI, was calculated by applying three ensemble models resulting from the combination of three statistical bivariate methods and the fuzzy AHP model.
The decision to apply three ensemble models was taken after a careful review of the literature according to which hybrid models have higher performance than stand-alone ones [15]. The models applied for the calculation of the FFPI showed that the hydrographic basin of the Izvorul Dorului River has a high and very high potential for a rapid surface runoff with a percentage between 38% and 58% of its surface. It was also highlighted that in particular, the upper and middle basin is characterized by these values of FFPI. Since the genesis of rapid water runoff on the slopes is finally reflected in the flooded areas along the rivers, it was decided to continue the study with the identification of valleys with a high potential for flash-flood propagation, along with the identification of the floodplains. In this regard, the three FFPI models were evaluated, and the result provided by FAHP-WOE, characterized by an AUC-ROC curve of 0.837 in the case of training data and 0.79 in the case of test data, was identified as the most accurate. Using the methodology proposed by Costache et al. [25], the valleys in the study area were identified and classified according to the degree of torrentiality. Valleys with a small and very small propagation potential were eliminated from the analysis, with only those characterized by a medium, high, and very high potential being used. The AHP model was further used to calculate the flood potential index along the torrential valleys and at the same time to determine the potential for flooding generated by the flash-flood propagation. It is worth mentioning that following the flash-flood genesis (which is facilitated by the torrential areas highlighted in Figure 1) and their propagation, the areas located along the torrential valleys are the most affected regions because the water flow from the slopes will be concentrated on the main river network. Therefore, it is very important to indicate the surfaces that are finally affected by these complex phenomena. This resulted in 24% of the delimited surface having a high and very high potential for flooding.
In a previous study, Costache et al. [58] estimated only the flooding susceptibility along the large river valleys within the Trotuș River basin from Romania, unlike the present study which analyzed the following three elements in close spatial and causal connection: (i) flash-flood potential at the slopes level; (ii) river valleys torrential degree; and (iii) flood potential along the river basins within this small catchment. In addition to this difference regarding the methodological approaches, the present study also differed from that conducted by Costache et al. [58] by the methods proposed for determining the susceptibility to the analyzed hydrological hazards. Thus, in the present study, three ensemble models of the fuzzy analytical hierarchy process with bivariate statistical methods for the estimation of flash-flood potential at the slopes level were applied, while in the previous study, three other ensemble models of the adaptive neuro-fuzzy inference system (ANIFS) were applied to determine the flooding potential at the large river valley level. Moreover, the flow accumulation procedure was applied in the present research in order to identify the torrential valleys. Another example where the fuzzy multicriteria decision making analysis was proposed to estimate the flood susceptibility was the study carried out by Azareh et al. [59]. In that research, which focused on the Haraz watershed in Iran, the authors used a combination between DEMATEL, analytical network process, and fuzzy logic in order to estimate the flood susceptibility. Like in the present case, the performance of the applied model was very good, which was revealed by an AUC-ROC curve between 0.8 and 0.9. Nevertheless, the main difference between the present study and the research work developed by Azareh et al. [59] is given by the fact that the mentioned study only included the evaluation of the terrain surface potential along the river valley, to produce the flooding phenomenon and did not also include an evaluation of the slopes for surface runoff genesis.

Conclusions
The assessment of flash-floods and flood susceptibility is an actual scientific topic due to the high potential of the studies to propose solutions for reducing the economic damage and diminishing the number of victims. The new approach developed in the present research is useful because it provides a complete overview regarding the susceptibility of the entire phenomenon composed of rapid surface runoff on the slopes, the propagation of flash-floods generated by the surface runoff, and the potential for flooding along torrential valleys. The water quality in the floodplains will be lower because the flashflood waves will be accompanied by the massive transport of materials from the slopes and inside the forest vegetation. Furthermore, the decision-makers will have a clearer image regarding the places they must implement measures to reduce the water runoff on the slopes, to arrange the torrential valleys, and to protect the areas exposed to floods.