**Assessing Land-Cover E**ff**ects on Stream Water Quality in Metropolitan Areas Using the Water Quality Index**

### **TaeHo Kim , YoungWoo Kim, Jihoon Shin , ByeongGeon Go and YoonKyung Cha \***

School of Environment Engineering, University of Seoul, 163, Seoulsiripdae-ro, Dongdaemun-gu, Seoul 02504, Korea; willy1995@uos.ac.kr (T.K.); youngwoo0508@uos.ac.kr (Y.K.); sjh3473@uos.ac.kr (J.S.); rhqudrjs7@uos.ac.kr (B.G.)

**\*** Correspondence: ykcha@uos.ac.kr

Received: 7 October 2020; Accepted: 21 November 2020; Published: 23 November 2020

**Abstract:** This study evaluated the influence of different land-cover types on the overall water quality of streams in urban areas. To ensure national applicability of the results, this study encompassed ten major metropolitan areas in South Korea. Using cluster analysis, watersheds were classified into three land-cover types: Urban-dominated (URB), agriculture-dominated (AGR), and forest-dominated (FOR). For each land-cover type, factor analysis (FA) was used to ensure simple and feasible parameter selection for developing the minimum water quality index (WQImin). The chemical oxygen demand, fecal coliform (total coliform for FOR), and total nitrogen (nitrate-nitrogen for URB) were selected as key parameters for all land-cover types. Our results suggest that WQImin can minimize bias in water quality assessment by reducing redundancy among correlated parameters, resulting in better differentiation of pollution levels. Furthermore, the dominant land-cover type of watersheds, not only affects the level and causes of pollution, but also influences temporal patterns, including the long-term trends and seasonality, of stream water quality in urban areas in South Korea.

**Keywords:** urban stream; factor analysis; land-cover type; metropolitan area; minimum water quality index; pollution

### **1. Introduction**

Global urbanization is an ongoing trend, with 55% to 68% of the world's population projected to reside in urban areas by 2050 [1]. Urbanization induces multiple stressors, especially land-use/land-cover changes such as deforestation and the growth of industrial and residential areas, resulting in increased impervious surfaces [2–5]. Consequently, urbanization leads to a deterioration of water quality in streams through an increase in pollution sources and various hydromorphological changes [6–8]. Despite their at-risk status, streams in urban areas are crucial water resources with a number of designated uses, such as drinking water supply, recreation, and wildlife conservation [9–12].

Therefore, it is vital to establish management strategies for preventing or alleviating water quality problems; this requires efforts to monitor and assess stream water quality in urban areas. The water quality index (WQI), an approach that quantitatively integrates a number of chemical, physical, and biological water quality parameters, has been widely used to assess the water quality status of both surface and groundwater systems [13–17]. In recent years the advent of big data and the accumulation of monitored multivariate data has prompted a substantial increase in the application of WQI to environmental and ecological studies [18–20]. In many of these studies, the developed WQI has been used to capture long-term trends [21,22], seasonal fluctuations [23,24], or spatial variations [25,26] in

the overall stream water quality in urban areas. As well as determining the spatiotemporal patterns of stream water quality in urban areas, previous WQI-based research has also determined pollution sources and anthropogenic effects [27–29] and selected the key parameters that represent variations in water quality [30–33].

Recent assessments of urban stream water quality have increasingly employed parameter selection using a number of statistical methods, highlighting the advantages of this process for cost and time saving during assessment. For example, Wu et al. [33] used stepwise multiple regression, which assumes linearity between the WQI and each parameter, to select five parameters representing the water quality of streams in the highly developed area of Lake Taihu Basin, China. Tripathi and Singal [31] used both principal component analysis (PCA) and correlation analysis to select nine parameters to develop a WQI for the Ganga River, which flows through some highly polluted cities of India. Moreover, linear discriminant analysis was applied by Han et al. [34] to select parameters that most effectively differentiate temporal groups (wet versus dry period) and spatial groups (east vs. west parts of the lake) in the Fu River and Baiyangdian Lake, both of which are located in a highly populated region of northern China.

However, the spatial scales of previous parameter selection studies have been limited to single water bodies or single basins; thus, the parameters selected in these studies have limited applicability to other urban stream ecosystems. Furthermore, the effects of different types of anthropogenic activities (e.g., industry, cultivation, or forestation), on stream water quality in urban areas has rarely been considered [26,35]. To overcome these limitations, this study presents the first attempt, to our knowledge, to explicitly account for the effects of different land-cover types on the water quality response and key water quality parameters of urban streams. This study was conducted on a national scale, encompassing a wide range of hydromorphological and geographical characteristics and socioeconomic backgrounds, which are also key factors influencing water quality [36–40]. Therefore, this study aimed to provide parameter selection results that are both informative and applicable to other unexplored streams in urban areas of South Korea.

Streams across ten major metropolitan areas of South Korea were investigated. Cluster analysis was performed to classify stream watersheds based on their land-cover characteristics. Then, the objective WQI (WQIobj) was calculated for each land-cover type using all available water quality parameters. The long-term trends of WQIobj were evaluated using the seasonal Mann-Kendall (SMK) test, and only periods exhibiting temporal stability were used in further analyses. For each land-cover type, key parameters were selected using factor analysis (FA) to develop the minimum WQI (WQImin). The objectives of this study were: (1) To assess the long-term trends and seasonality of the overall stream water quality in metropolitan areas in South Korea; (2) to analyze how different land-cover types affect stream water quality in urban areas and key water quality parameters; and (3) to evaluate the correlation between WQIobj and WQImin and relationships between WQImin and land-covers.

### **2. Materials and Methods**

### *2.1. Study Area and Data Description*

Ten major metropolitan areas across South Korea, with populations of greater than one million, were included in this study (Seoul, Busan, Incheon, Daegu, Daejeon, Gwangju, Suwon, Ulsan, Changwon, and Goyang (Figure 1)) [41]. Within the study area, 81 water quality monitoring sites were selected at tributaries that directly or indirectly flow into either the Han, Geum, Nakdong, or Yeongsan Rivers, the four major rivers of South Korea. The selected monitoring sites covered 35 standard watersheds with the range of watershed area from 39 to 294.9 km<sup>2</sup> , and a mean area of 103.29 km<sup>2</sup> , the smallest unit of the drainage area division system in South Korea (http://wamis.go.kr). Water quality data were provided by the National Institute of Environmental Research of the Ministry of Environment (http://water.nier.go.kr). The data spanned the time period from 2007 to 2018, and the monitoring frequency varied by site from weekly to monthly. Among the 54 water quality parameters

initially included in the data, heavy metals and other toxic chemicals, such as mercury, cadmium, arsenic, and cyanide, were not included because at least 99.5% of the values for these parameters were either missing or below the detection limit. Furthermore, parameters without available reference values were not included in the analyses. The reference values (i.e., normalization factors and weights) required to develop the Bascarón WQI were provided by previous studies [27,42–44].

**Figure 1.** Location of monitoring sites in ten major metropolitan areas of South Korea.

Fourteen water quality parameters were included in the analyses: Water surface temperature (Temp), electrical conductivity (EC), pH, dissolved oxygen (DO), five-day biochemical oxygen demand (BOD5), chemical oxygen demand (COD), suspended solids (SS), total nitrogen (TN), ammonium nitrogen (NH<sup>4</sup> <sup>+</sup>-N), nitrate nitrogen (NO<sup>3</sup> −-N), total phosphorus (TP), orthophosphate phosphorus (PO<sup>4</sup> <sup>3</sup>−-P), total coliform (TC), and fecal coliform (FC) (Table 1). Among the 81 monitoring sites initially selected for our study, 58 were included for the water quality assessment as they had measurements for all 14 water quality parameters. Land-cover data were provided by the Environmental Geographic Information System; the year of data collection varied from 2010 to 2018 depending on the region (https://egis.me.go.kr). The land-cover data involved seven categories: urban (or built-up) land, agricultural land, forested land, grassland, wetland, barren land, and water. For each of the 35 watersheds, the relative proportions of the seven land-cover categories were calculated using QGIS 2.18.16 [45] and ArcGIS 10.3 software [46].




### *2.2. Statistical Analyses*

### 2.2.1. Cluster Analysis (CA)

CA is an unsupervised pattern recognition technique, whereby individual objects are grouped into a number of clusters whose objects are more similar than those in other clusters. Among the available CA methods, hierarchical agglomerative CA (HACA) was used in this study. HACA is a successive process, in which two objects in the closest proximity form a cluster at the lowest hierarchy. In the next step, the newly generated two clusters in the closest proximity form a combined cluster. Merging continues until all objects are linked to form a single cluster at the highest hierarchy. The squared Euclidean distance was used as a measure for calculating the proximity between objects/clusters.

Furthermore, we employed the Ward's minimum variance linkage function, which uses distance information to merge objects into a hierarchical cluster tree and is visually represented by a dendrogram [47]. As HACA results in a single cluster, the dendrogram needs to be divided at a specific height to generate multiple clusters. The height in the dendrogram can be defined as (Dlink/Dmax)·100, where Dlink is the linkage distance for a pair of objects/clusters and Dmax is the maximum linkage distance. According to previous studies, the height for dendrogram partitioning was set to 60; that is, (Dlink/Dmax)·100 > 60 [48,49]. The CA was performed using 'dendrogram' function from the 'SciPy' library [50] in Python 3.6 [51]. To generate clusters based on land-cover type, HACA was performed using the relative proportions of the six land-cover types for each standard watershed (excluding water) as variables.

The differences in water quality parameters and WQI among different clusters were assessed using summary statistics and non-parametric tests, i.e., Kruskal Wallis *H* and Mann-Whitney *U* tests. Non-parametric tests were selected due to the non-normality of water quality parameters. The Kruskal Wallis *H* test examined the differences in distributions for the three clusters. When the significant differences occurred, as a post-hoc analysis, the Mann-Whitney *U* test was used to identify which cluster(s) revealed the significant difference in distribution from the other cluster(s). The Kruskal Wallis *H* and Mann-Whitney *U* tests were performed using 'kruskal' and 'mannwhitneyu' from the 'SciPy' library [50] Python 3.6 [51]. Statistical significance was indicated by *p*-value < 0.05.

### 2.2.2. Water Quality Index (WQI) Development

The method for WQI development used in this study builds on the WQIobj [43], a modification of the *Bascarón* WQI, also known as subjective WQI [13], which excludes the constant term multiplied to WQIobj, which reflects the subjective judgment of overall water quality. The WQIobj is calculated as follows,

$$\text{WQI}\_{\text{obj}} = \frac{\sum\_{i=1}^{n} \text{C}\_{\text{i}} \text{P}\_{\text{i}}}{\sum\_{i=1}^{n} \text{P}\_{\text{i}}} \tag{1}$$

where n is the number of available water quality parameters, C<sup>i</sup> is a normalization factor that converts the value of a parameter into a common scale ranging from 0 to 100 with an interval of 10 (Table 1), Pi is the weight indicating the relative importance of parameters, which ranges from 1 to 4 (Table 1), and WQImin is a simplification of WQIobj indicating the minimum WQI [42,43] and is calculated as,

$$\text{WQI}\_{\text{min}} = \frac{\sum\_{i=1}^{n\_{\text{min}}} \mathbf{C}\_{i}}{\mathbf{n}\_{\text{min}}} \tag{2}$$

Note that Equation (2) for WQImin does not include the weight term, indicating that the parameters included in WQImin assessment are considered equally important. Here, nmin is the number of key parameters, which is a subset of all n available parameters. The WQIobj and WQImin scores were graded into five classes to indicate the overall water quality status: excellent (91–100), good (71–90), medium (51–70), bad (26–50), and very bad (0–25) [42,43,52]. Also, when comparing WQImin with

WQIobj for evaluating whether they are well-correlated, linear regression (WQIobj = a·WQImin + b) was performed using 'linregress' function from the 'SciPy' library [50] in Python 3.6 [51].

### 2.2.3. Seasonal Mann-Kendall (SMK) Test

The Mann-Kendall (MK) test is a non-parametric test that assesses if the temporal trend of a variable exhibits a monotonic increase or decrease [53,54]. The SMK test is an extension of the MK test that accounts for the effect of seasonality by performing the test separately for each pre-defined season [55]. In this study, the SMK was employed to identify the point in time at which WQI no longer shows a significant increasing or decreasing trend; this stabilized time period was divided into training and test sets, and further analyses were performed. CA results were used to assess the trend of monthly averaged WQIobj for each land-cover cluster, initially using the data for the entire time period (2007–2018). When the trend showed a significant increase or decrease (*p*-values < 0.05), the SMK was performed excluding a 1-year period of data from the starting year. The test was iteratively performed until the trend appeared to be insignificant. The SMK test was performed using 'seasonal\_test' function from the 'pyMannKendall' library [56] in Python 3.6 [51].

### 2.2.4. Factor Analysis (FA)

FA attempts to account for the structure (i.e., correlation and variation) of data, consisting of measured variables with a reduced number of factors, which are also termed latent variables. Exploratory FA (EFA), which does not assume a priori relationships among the factors and measured variables, was used to reveal the underlying factors behind the correlations among measured water quality parameters. Contrary to confirmatory FA, EFA does not posit any relationship between specific factors and measured variables. Therefore, it suits the purpose of this analysis.

Prior to analysis, the values of all water quality parameters except Temp and pH were log-transformed, and the values of all parameters were standardized to have a distribution with a mean of zero and standard deviation of one. To examine whether water quality data were suitable for FA, the Kaiser-Mayer-Olkin (KMO) test [57] and Bartlett's test [58] were performed. The FA was assumed to be valid when the KMO value exceeded 0.5 and the Bartlett's test result was significant (*p*-value < 0.05).

To determine the number of factors retained in the FA, Horn's parallel analysis (PA) was used [59]. PA compares the eigenvalues (which indicate the relative importance of a factor in explaining the variance of measured variables) from measured data with the eigenvalues from random data, which have the same sample size and number of variables as the measured data and are obtained using a Monte-Carlo simulation. The differences between the eigenvalues from measured data and the mean eigenvalues from random data were calculated. Factors with differences greater than zero were retained in the FA.

As a method for factor extraction, principal component analysis (PCA) was used [60,61]. The maximum likelihood method, another common method for factor extraction, was not selected because of its multivariate normality requirement, which is often not met for water quality parameters even after the transformation (e.g., log-transformation) of values. Squared factor loading, which reflects the proportion of variance in a measured variable explained by each factor, was calculated as a result of PCA implementation. The communality was calculated by summing the squared factor loadings of a given variable across all factors to indicate the proportion of variance in a measured variable explained by all factors. Moreover, the uniqueness was calculated by subtracting the communality from the total variance of a variable.

Factor rotation (i.e., the change in the axes of factors) was implemented to yield interpretable factors by attaining a simple structure for factor loadings. Without rotation, most variables load heavily onto the first and early factors, whereas rotation yields a simple structure in which each variable loads heavily onto only one factor, while loading lightly onto the other factors. Varimax rotation was used as a rotation method, which is a common type of orthogonal rotation. Orthogonal rotation assumes

that factors remain uncorrelated with one another. FA was performed using 'principal' function from 'psych' [62] packages in R 3.5.3 [63].

For each land-cover type, the water quality parameter that showed the highest loading factor associated with each retained factor was interpreted as a key water quality parameter. Accordingly, the number of retained factors corresponded to the number of key parameters representing the overall stream water quality in urban areas. The FA procedure was performed for each land-cover type determined by CA, and the key parameters for different land-cover types were used for the WQImin calculation.

### **3. Results**

### *3.1. Land-Cover Characteristics of Metropolitan Areas in South Korea*

Using the HACA, three clusters were generated based on the land-cover characteristics of 35 watersheds in ten major metropolitan areas (Figure 2a). Notably, each of the watersheds included in each of the three clusters had a single dominant land-cover: Urban, agriculture, and forest, respectively (Figure 2b, Table S1). The mean proportion of urban land for the 15 watersheds with urban-dominated land-cover (URB) was 0.50 (± one standard deviation of 0.12), which was higher than that of agricultural (0.06 ± 0.05) and forested (0.30 ± 0.10) land. In contrast, the five watersheds with agriculture-dominated land-cover (AGR) had a mean relative area of 0.44 (± 0.08) for agricultural land-cover, which was more dominant than urban (0.16 ± 0.07) and forested (0.24 ± 0.07) land-cover. The 15 watersheds with forest-dominated land-cover (FOR) were mainly composed of forested land, with a mean proportion of 0.60 (± 0.08), whereas the proportion of urban (0.12 ± 0.06) and agricultural (0.16 ± 0.04) land was relatively minor.

**Figure 2.** Clustering results of 35 watersheds, named metropolitan area with numbering, based on six land-cover categories. (**a**) Dendrogram exhibiting three clusters generated from hierarchical agglomerative cluster analysis. The horizontal dashed gray line represents the height for dendrogram partitioning, (Dlink/Dmax)·100 > 60. (**b**) Percentage (%) of the dominant land-cover type for each of the three clusters. The red circle, yellow triangle, and green square denote watersheds that are urban-dominated, agriculture-dominated, and forest-dominated, respectively.

The three land-cover types (URB, AGR, and FOR) were unevenly distributed across the metropolitan areas. Among the URB, 73.3% were concentrated in Seoul (nine watersheds) and its adjacent cities, Suwon (one watershed) and Incheon (one watershed). Three of the five AGR were located in Gwangju, whereas the other two were located in Busan and Changwon. The spatial distribution of FOR was also concentrated, with 33.3% in Daejeon and 26.7% in Daegu.

### *3.2. Land-Cover E*ff*ects on Stream Water Quality in Urban Areas*

The long-term trends of overall water quality calculated using all available parameters (WQIobj), based on the results of SMK tests, differed by land-cover type (Figure 3). For URB, WQIobj values gradually improved until becoming stable in 2015 (Figure 3a). In comparison, WQIobj values for AGR showed a greater improvement in early years before becoming stable in 2012 (Figure 3b). For FOR, WQIobj values did not show any significant trend during the entire period from 2007 to 2018 (Figure 3c). In more recent years (2015–2018), during which all land-cover types exhibited a stable trend, the overall water quality was worst for URB (*p*-values < 0.05 as a result of Kruskal Wallis *H* and Mann-Whitney *U* tests), as indicated by lower WQIobj values (75.04 ± 9.90) than those for AGR (78.91 ± 8.31) and FOR (82.82 ± 7.97). Regardless of the land-cover type and time period, WQIobj values tended to be lower during the wet season (July to September) than during the dry season (Figure 3).

**Figure 3.** Long-term (2007–2018) trends of objective water quality index (WQIobj) for watersheds with (**a**) urban-dominated, (**b**) agriculture-dominated, and (**c**) forest-dominated land-cover. Blue and red circles denote the mean monthly WQIobj for dry and wet seasons, respectively, and vertical lines denote one standard deviation of monthly WQIobj. The gray area represents a period exhibiting no significant increase or decrease in WQIobj based on the results of seasonal Mann-Kendall tests.

The land-cover types of the watersheds influenced most water quality parameters in urban streams except for pH, EC, DO, and PO<sup>4</sup> <sup>3</sup>−-P, which were similar regardless of the dominant land-cover (Table 2). Compared with URB and AGR, FOR exhibited the lowest level of contamination for the majority of water quality parameters. The level of contamination between URB and AGR differed depending on the water quality parameter. In terms of nitrogen (TN and NO<sup>3</sup> −-N) and microbiological indicators (TC and FC), the streams in URB exhibited significantly worse conditions than those in AGR (Table 2). On the other hand, indicators for organic matter (BOD<sup>5</sup> and COD) and turbidity (SS) indicated significantly higher levels of water contamination in AGR than URB (Table 2).




### *3.3. Key Water Quality Parameters for Di*ff*erent Land-Cover Types*

The water quality data were suitable for the application of FA, as indicated by the results of the KMO test (0.82 for URB, 0.67 for AGR, and 0.73 for FOR) and Barlett's test (*p*-value < 0.05 for all land-cover types). To perform the FA, the data measured during the more recent years (2015–2018), when the WQIobj values stabilized for all land-cover types, were divided into training (2015–2016) and testing (2017–2018) data sets. The results of FA using the training data indicated that three factors apiece should be retained for URB, AGR, and FOR (Table S2). For each land-cover type, the water quality parameters with the highest factor loading, associated with each of the retained factors, were selected as the key parameters for the WQImin calculation (Table 3). Frequently, for a given factor, more than one water quality parameter had a factor loading greater than 0.75 [64], which is indicative of a strong correlation between the factor and the parameter (Table 3). In such cases, the parameters were generally highly correlated to each other, with a Pearson's correlation coefficient ranging from 0.49 to 0.88 (Figure S1). Consequently, the three key parameters selected for URB were COD, FC, and NO<sup>3</sup> −-N, in order of corresponding factors (Table 3). Three parameters were selected for AGR were FC, COD, and TN (Table 3). The three parameters selected for FOR were COD, TN, and TC (Table 3).

### *3.4. Comparison between WQIobj and WQImin*

Using the test data, the relationships between monthly WQImin and WQIobj values were assessed; WQImin and WQIobj generally exhibited moderate to strong, linear relationships with R<sup>2</sup> values of 0.66 for URB, 0.78 for AGR, and 0.73 for FOR (Figure 4). For both WQIobj and WQImin, URB was generally associated with the poorest overall water quality, with mean WQI values of 75.79 and 67.20, respectively. Further, based on both WQIobj and WQImin, the overall water quality for AGR (mean WQI values of 78.86 and 73.39) was generally poorer than that for FOR (mean WQI values of 82.41 and 77.41). The location of intersection, where the regression line and one-to-one line cross, differed by land-cover type: 87.48 for URB, 81.98 for AGR, and 88.14 for FOR (Figure 4). Below the intersection, WQIobj values tended to be higher than WQImin scores, whereas the opposite was true above the intersection (Figure 4). As the proportion of values below the intersection was greatest for URB, the positive difference between the mean WQIobj and WQImin values for URB (8.59) was greater than that for AGR (5.47) and FOR (5.00). Within each land-cover type, the variation of WQImin values, with one standard deviation of 13.61 for URB, 13.05 for AGR, and 12.09 for FOR, was greater than the variation of WQIobj values, with one standard deviation of 9.27 for URB, 8.62 for AGR, and 7.97 for

FOR. Note that the degree of variation in WQI values, in descending order, was URB, AGR, and FOR for both WQIobj and WQImin.

**Table 3.** Factor loadings for 14 water quality parameters for watersheds with urban-dominated (URB), agricultural-dominated (AGR), and forest-dominated (FOR) land-cover. Asterisks (\*) indicate a factor loading greater than 0.75 or the highest factor loading in the factor. Var (%) represents the explained variance of total variance for each factor.


**Figure 4.** Relationships between objective water quality index (WQIobj) and minimum WQI (WQImin) for watersheds with, (**a**) urban-dominated, (**b**) agriculture-dominated, and (**c**) forest-dominated land use. Black circles denote WQI values calculated using the testing data set (2017–2018). Black dotted and blue dashed lines represent one-to-one, and regression lines, respectively. Red square represents the point of intersection between the one-to-one line and regression line.

### *3.5. Spatial Distribution of Overall Stream Water Quality in Urban Areas*

WQI values by site, calculated for 2015–2018, indicated that WQIobj and WQImin values were highly linearly correlated, with an R<sup>2</sup> value of 0.84 (Figure 5b). However, there was a clear tendency for WQIobj values to be higher than WQImin values (Figure 5a,b). The difference in values between WQIobj and WQImin led to differences in WQI classification in 25.9% of the 58 monitoring sites (Figure 5a). In Seoul, the change in calculation method from WQIobj to WQImin yielded a change of classification from good to medium in 33.3% of 18 monitoring sites. In the other five metropolitan areas (i.e., Daejeon, Gwangju, Daegu, Busan, and Ulsan), a change in classification occurred in one or two sites, accounting for 7.7–40.0% of the sites in each area (Figure 5a). In the remaining five metropolitan areas (i.e., Goyang, Suwon, Incheon, and Changwon), no change in WQI classification occurred in response to application of the WQImin (Figure 5a).

### *3.6. Seasonality of Overall Stream Water Quality in Urban Areas*

From 2015 to 2018, the monthly patterns of overall water quality calculated using WQImin differed by land-cover type (Figure 6). For URB, which exhibited the worst overall water quality, the proportion of WQImin values corresponding to equal to or worse than medium status increased during the wet season (July to September), whereas the proportion of good to excellent status sites increased during the dry season (all other months) (Figure 6a). For FOR, the WQImin status was consistently better than or equal to medium, and the proportion of medium status sites increased during the wet season (Figure 6c). For AGR, the WQImin status tended to worsen during the wet season, with an increase in the proportion of medium status sites; however, this seasonality was less consistent compared with other land-cover types (Figure 6b).

**Figure 5.** Spatial distribution of the water quality index (WQI) in ten major metropolitan areas of South Korea. (**a**) Mean objective WQI (WQIobj) and minimum WQI (WQImin) values and grades from 2015 to 2018 for each of the 58 monitoring sites. (**b**) Relationship between mean WQIobj and WQImin values.

**Figure 6.** Monthly distribution (%) of minimum water quality index (WQImin) grades from 2015 to 2018 for watersheds with, (**a**) urban-dominated, (**b**) agriculture-dominated, and (**c**) forest-dominated land-cover. Month names for dry and wet seasons are colored blue and red, respectively.

### **4. Discussion**

### *4.1. Suitability of FA as a Parameter Selection Method*

In this study, FA, which involves factor extraction and rotation processes, was used to reduce multiple intercorrelated physical, chemical, and biological water quality parameters into a smaller number of latent factors, and to select key water quality parameters that had the strongest correlation with a given latent factor. In previous studies, along with subjective judgments [27,33,65–69], multivariate statistical techniques were employed to select parameters on an objective basis. For example, stepwise multiple regression has been used [33,69,70] to determine the set of parameters that could best explain the variance of WQIobj. Compared with unsupervised learning (e.g., FA), regression is a supervised method that requires reference values; in this case, WQIobj values for training data. However, because of the multi-collinearity and the resulting bias, WQIobj is not often a suitable reference.

Furthermore, previous studies have used PCA at the first step followed by Pearson's correlation analysis to extract water quality parameters that showed high contributions to selected components and low correlations with other parameters [31,66]. Post-hoc correlation analysis was required, since few first factors derived from PCA are strongly associated with most of the correlated parameters. Therefore, the application of PCA alone is not sufficient to attain key parameters that represent extracted factors. To address this limitation, in this study PCA was conducted in conjunction with factor rotation, which yields a simple structure for the factor loading matrix, in which only a small number of variables have high loadings onto a given factor and do not overlap among the factors. As a result, parameters with high loadings on a given factor appear to be more distinct and homogeneous. Therefore, a set of parameters with high loadings across all factors are expected to represent multifaceted aspects of water quality. Furthermore, the use of varimax rotation as a factor rotation method ensures the extracted factors are uncorrelated with one another, facilitating the selection of key parameters, the relationships among which can be assumed to be independent. Therefore, factor rotation used in conjunction with PCA does not require subsequent correlation analysis, which simplifies parameter selection to a single-step process.

### *4.2. Key water Quality Parameters*

Selected key water quality parameters were similar among different land-cover types (COD, FC, and NO<sup>3</sup> −-N for URB; COD, FC, and TN for AGR; COD, TC, and TN for FOR), indicating that the relationships among parameters were consistent regardless of land-cover type. For example, across all land-cover types, COD, BOD5, and SS were closely correlated (Figure S1) and had high loadings with the same factor (Table 3). The high correlations were shown, since the three parameters commonly account for biodegradable organic matter. In addition, for one being the subset of the other, FC and TC, and NO<sup>3</sup> −-N and TN, were closely related to each other (Figure S1) and had the highest loadings onto the same factor for all land-cover types (Table 3). Note that phosphorus parameters showed moderate to strong associations with TC and FC within the same factor for all land-cover types (Table 3). Therefore, rather than phosphorus parameters, either TC or FC, which showed higher loadings with the factor than the phosphorus parameters, was selected as the key parameter. A possible speculation over this co-occurrence tendency is that phosphorus and fecal indicator bacteria may originate from the same pollution source (e.g., domestic sewage and agricultural runoff) or the same mechanism (e.g., sediment release), but future research will be necessary for interpreting the causal relationships.

The presence of multiple parameters with almost equally high loadings onto a given factor necessitated comparisons between WQImin and modified WQImin, in which a key parameter (e.g., COD) is replaced by its surrogate parameter (e.g., BOD5) that was strongly related to the key parameter within the same factor. The results illustrated that modified WQImin was generally in close agreement with WQImin (Figure S2), suggesting that a set of parameters that shows high loadings within the same factor can be used interchangeably. Note that, compared with other sets of parameters, linear relationships

between WQImin and modified WQImin for fecal indicator bacteria were weaker because of the large variability inherent in FC and TC concentrations. Nonetheless, given the marginal differences in factor loading between TC and FC regardless of land-cover type (Table 3), between the two parameters, the key parameter should be selected depending on management focus or data availability.

The results of FA need to be interpreted and applied with care. The factor extraction process of FA determines the factors worth retaining, and the subsequent factor rotation, whereby the factors become least correlated with each other, yields the proportion of variance explained by a given factor to be distributed more evenly among the factors. Therefore, it is not particularly valid to prioritize the factors and the consequent key parameters. Instead, the selected key parameters should be considered independently of each other and as equally important. In this regard, assigning different weights to key water quality parameters with equal importance should not be included as a step for WQImin development. Previous studies reported that using weights improved the linearity between WQImin and WQIobj [33,70]. In contrast to these findings, we found that the use of weights, which were estimated based on two methods, the relative weight [33,70] and the percentage of variance explained by the given factor (Table 3), yielded only slight differences in the WQImin-WQIobj relationships (Figure S3).

It should be acknowledged that the water quality data, used in this study, did not include several widely measured parameters, such as parameters for minerals, salts, metals and flow rate. If such parameters were added to the data, FA may include additional factors and key parameters. Moreover, the results of parameter selection did not contain the basic water quality parameters of Temp and pH in the key parameter list for any land-cover type. In addition, despite being frequently included as a key parameter [42,43,68–70] in previous studies, DO was not selected for any land-cover type in this study (Table 3). Variations in Temp, pH, and DO may be influenced by anthropogenic activities but are also attributable to natural variability. That is, they exhibit diurnal fluctuations and are strongly influenced by meteorological conditions [33,65]. Our results suggest that Temp, pH, and DO, whose patterns are substantially influenced by natural variations, may not successfully capture the total variance of stream water quality in urban areas, and may not be suitable for being included as key parameters.

### *4.3. Comparison between WQImin and WQIobj*

Our results of test data showed that WQImin and WQIobj have close linear relationships across all land-cover types (Figure 4), suggesting that WQImin can be used to predict WQIobj using the established regression model. However, WQImin values tended to be higher than WQIobj above a certain threshold and lower than WQIobj below this threshold. This tendency indicates that the use of WQImin eliminates the "eclipse effect" [71], which arises from the redundancy inherent in WQIobj; accordingly, WQIobj is subject to overestimating bad water quality status and underestimating good water quality status. The removal of redundancy was also evidenced by the larger variance of WQImin compared with that of WQIobj for all land-cover types (Figure 4). Therefore, the development and use of WQImin is expected to improve the identification of the overall water quality status and the level of water pollution in streams across urban areas. Our results demonstrate that the method selection for WQI assessment has important resource and management implications. Changing the method from WQIobj to WQImin altered the spatial distribution of the overall water quality status; this status change occurred in a minor to substantial portion of monitoring sites, depending on the metropolitan area (Figure 5). This change suggests that the use of WQImin instead of WQIobj, which may involve a status change from "good" to "medium" or vice versa, may affect priority setting and resource allocation among individual watersheds or groups of watersheds.

### *4.4. Land-Cover E*ff*ects on Stream Water Quality in Urban Areas*

Our results indicate that the dominant land-cover affected the overall stream water quality in urban areas, with mean values of both WQIobj and WQImin decreasing in the order: FOR > AGR > URB (Figure 4). The dominant land-cover type also contributed to the deterioration of differing water quality parameters (i.e., nitrogen and microbiological indicators for URB, but organic matter and turbidity for AGR) (Table 2). The long-term trends of overall water quality differed by land-cover type (Figure 3). Over the last decade, WQIobj trends for URB and AGR exhibited early improvement before becoming stable, whereas the trend for FOR did not change significantly (Figure 3). These patterns support that, across the country, management programs implemented to control point or non-point sources for URB and AGR were effective in improving overall stream water quality [72–75]. Moreover, the implementation of conservation measures against continuing development pressures in metropolitan areas played a role maintaining the water quality in FOR. Furthermore, the land-cover type exerted an influence on the seasonality of overall water quality (Figure 6). In recent years (2015–2018), the seasonal patterns of WQImin have differed for URB and FOR, whereas AGR exhibited less obvious seasonality. The less consistent seasonality for AGR may be partly attributable to the small sample size (*n* = 287, compared with n for URB = 1881 and *n* for FOR = 1162) corresponding to AGR. During the wet season, both URB and FOR exhibited a negative change in overall water quality with an increase in the proportion of "medium" and "good" status sites relative to "excellent" status sites (Figure 6). For URB with typically high proportions of impervious surfaces, stormwater runoff may play a significant role in decreasing overall water quality during the wet season [76–78]. Moreover, an increase in sediment discharge as well as sediment perturbation with rainfall events may facilitate the release of pollutants into surface water [79–82], resulting in a decrease in overall water quality during the wet season in both URB and FOR. In contrast, subsequent to the wet season, when dilution effects can occur [83–85], URB alone exhibited an increase in the proportion of "bad" status sites relative to "medium" and "good" status sites (Figure 6). This indicates that, not only non-point sources, but also point sources, such as wastewater treatment plant effluent, are significant forms of pollution for URB.

### **5. Conclusions**

This study provided a statistical framework for implementing parameter selection in order to develop an objective WQImin in a single-step process. Comparisons between WQIobj and WQImin suggested that WQImin calculated with the key parameters yielded comparable results to WQIobj. Furthermore, WQImin reduced the eclipse effects arising from the use of correlated parameters for water quality assessment to result in a better differentiation between good and bad water quality statuses. These results have implications for management authorities, especially those motivated to launch their own monitoring network system but who have limited available resources. In this context, our results can be used to reduce monitoring demands by prioritizing the monitoring importance of a minimal number of water quality parameters. The results of WQImin confirmed that the dominant land-cover type of watersheds influence multidimensional aspects of urban stream water quality; namely, the overall degree and level of pollution as well as long-term and seasonal patterns. To confirm our results, future studies should expand the number of water quality parameters exhibiting various characteristics.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/11/3294/s1. Figure S1: Matrices of the Pearson's correlation coefficient for the period 2015–2016 among 14 water quality parameters for (a) urban-dominated (URB), (b) agricultural-dominated (AGR), and (c) forest-dominated (FOR) land-cover. Water quality parameters with high factor loadings (>0.75) on the same factor are outlined in the same color, Figure S2: Relationships between the minimum water quality index (WQImin) and modified WQImin from 2015 to 2018. To develop the modified WQImin, key parameter values were predicted using the established linear relationship between a key parameter and a surrogate parameter. Then, predicted values were converted into normalization factors for WQImin calculation. In the x-axis label, WQImin (COD → BOD<sup>5</sup> ) indicates that biochemical oxygen demand (BOD<sup>5</sup> ) was used as the surrogate for the key parameter of chemical oxygen demand (COD). Black dotted lines indicate 1:1 lines. Figure S3: Relationships between objective and minimum water quality indices (WQIobj and WQImin) from 2017 to 2018. Weights were determined using two methods; for a-c, a relative weight was assigned to each key parameter and for d-f, the percent variance explained by a given extracted factor was assigned to each key parameter. Black dotted lines and blue dashed lines indicate 1:1 lines and regression lines, respectively. Table S1: Proportions of three land-cover categories (urban, agricultural, and forested land) for urban-dominated watersheds (URB), agricultural-dominated watersheds (AGR), and forest-dominated watersheds (FOR). Table S2: Parallel analysis results comparing eigenvalues and simulated mean eigenvalues for urban-dominated (URB), agriculture-dominated (AGR), and forest-dominated (FOR) land-cover. The simulated mean eigenvalue indicates the mean eigenvalue calculated from randomly generated simulation data. Asterisks (\*) indicate that the eigenvalue is higher than the corresponding simulated mean eigenvalue.

**Author Contributions:** Conceptualization, T.K. and Y.C.; data curation, T.K., Y.K., J.S. and B.G.; formal analysis, T.K.; funding acquisition, T.K., Y.K., J.S., B.G. and Y.C.; Investigation, T.K., Y.K., J.S. and B.G.; methodology, T.K. and Y.C.; project administration, T.K. and Y.C.; resources, T.K., Y.K., J.S. and B.G.; software, T.K. and Y.C.; supervision, T.K. and Y.C.; validation, T.K. and Y.C.; visualization, T.K. and Y.C.; writing—original draft, T.K. and Y.C.; writing—review and editing, T.K. and Y.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2020R1A2C1009961).

**Acknowledgments:** The authors would like to thank the editors and reviewers for useful comments which are helpful in improving the manuscript quality.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

**Romulus Costache 1,2,3,† , Alina Barbulescu 3,† and Quoc Bao Pham 4,\***


**Abstract:** 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.

**Keywords:** flooding; flash-floods; bivariate statistics; fuzzy multicriteria decision-making; small catchments; Romania

### **1. Introduction**

According to Hu et al. [1], a total number of 2.5 billion peoples were affected by flashfloods and floods between 1994 and 2013. In the same period, 0.16 billion fatalities occurred due to the same natural risk phenomena. Therefore, since flash-floods are extremely severe phenomena, they are also very dangerous for human life [2,3]. These phenomena appear most frequently in small river basins characterized by a high slope. Additionally, areas with smaller slopes favor the accumulation of the transported water and materials [4]. In this context, the identification of sections that favor the surface runoff occurrence, torrential valleys on which the flash-floods are propagated, and the flood susceptibility assessment in those regions, is one of the most important measures to combat the negative effects of these phenomena on water quality and human society. Additionally, the results provided by this type of analysis are very important in assessing a region's vulnerability and risk to flash floods. It should be noted that most of the procedures regarding the evaluation of flashflood and flood risk assessment, which were adopted by the European countries includes

### **Citation:** Costache, R.;

Barbulescu, A.; Pham, Q.B. Integrated Framework for Detecting the Areas Prone to Flooding Generated by Flash-Floods in Small River Catchments. *Water* **2021**, *13*, 758. https://doi.org/10.3390/w13060758

Academic Editors: Marco Franchini and Salvador García-Ayllón Veintimilla

Received: 24 November 2020 Accepted: 9 March 2021 Published: 11 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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

### **2. 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<sup>2</sup> , 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. 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.

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

*Water* **2021**, *13*, x FOR PEER REVIEW 3 of 25

**Figure 1.** The study area. (**a**) Study area location in Romania; (**b**) torrential area location. **Figure 1.** The study area. (**a**) Study area location in Romania; (**b**) torrential area location.

#### **3. Methods 3. Methods**

*3.1. Background of the Models 3.1. Background of the Models*

#### 3.1.1. Statistical Index 3.1.1. 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: 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:

$$\mathcal{W}\_{lj} = \ln\left(\frac{f\_{ij}}{f}\right) = \ln\left(\frac{\frac{Npix(S\_i)}{Npix(N\_l)}}{\frac{\sum Npix(S\_i)}{\sum Npix(N\_l)}}\right) \tag{1}$$

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.

### 3.1.2. 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]:

$$FR = \frac{\frac{Np\dot{x}(1)}{Np\dot{x}(2)}}{\frac{\sum Np\dot{x}(3)}{\sum Np\dot{x}(4)}}\tag{2}$$

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; ∑ *Npix*(3) is the total number of torrential pixels within the study area; and ∑ *Npix*(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.

### 3.1.3. Weights of Evidence

Weights of evidence (WOE) is a widely used statistical model for landslide, flood, and fire forest susceptibility assessment [29–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<sup>+</sup>* ) 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–31]:

$$\mathcal{W}^+ = \ln \frac{P\{B|S\}}{P\{B|\overline{S}\}} \tag{3}$$

$$\mathcal{W}^- = \ln \frac{P\{\overline{B} | \mathcal{S}\}}{P\{\overline{B} | \overline{\mathcal{S}}\}} \tag{4}$$

where *B* and *B* are the presence and absence of flood conditioning parameters, respectively; *P* is the probability; and *S,* and *S* 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]:

$$\mathcal{W}^{+} = \ln \frac{\frac{Npi\dot{x}\_1}{Npi\dot{x}\_1 + Npi\dot{x}\_2}}{\frac{Npi\dot{x}\_3}{Npi\dot{x}\_3 + Npi\dot{x}\_4}} \tag{5}$$

$$\mathcal{W}^{-} = \ln \frac{\frac{Np \text{ix}\_2}{Np \text{ix}\_1 + Np \text{ix}\_2}}{\frac{Np \text{ix}\_4}{Np \text{ix}\_3 + Np \text{ix}\_4}} \tag{6}$$

where *W<sup>+</sup>* and *W*<sup>−</sup> are the positive and negative weights, respectively; *Npix*<sup>1</sup> and *Npix*<sup>2</sup> are the number of pixels with flood points inside *and* outside of the class, respectively; and *Npix*<sup>3</sup> and *Npix*<sup>4</sup> 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]:

$$\mathcal{W}f = \mathcal{W}plus + \mathcal{W}mintotal - \mathcal{W}min \tag{7}$$

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.

### 3.1.4. Fuzzy Analytical Hierarchy Process

The analytical hierarchy process (AHP) is an algorithm used for flood, landslide, and fire forest susceptibility mapping [32–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]:

$$A' = \begin{bmatrix} 1' & a'\_{12} & \cdots & a'\_{1n} \\ a'\_{21} & 1' & \cdots & a'\_{1n} \\ \vdots & \vdots & \ddots & \vdots \\ a'\_{n1} & a'\_{n2} & \cdots & 1' \end{bmatrix} = \begin{bmatrix} 1' & a'\_{12} & \cdots & a'\_{1n} \\ 1/a'\_{21} & 1' & \cdots & a'\_{1n} \\ \vdots & \vdots & \ddots & \vdots \\ 1/a'\_{n1} & 1/a'\_{n2} & \cdots & 1' \end{bmatrix} \tag{8}$$

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:

$$\sigma'\_i = \left(a'\_{i1} \otimes a'\_{i2} \otimes \dots \otimes a'\_{in}\right)^{\frac{1}{n}},\tag{9}$$

$$w\_i' = r\_i' \otimes \left(r\_1' \otimes \dots \otimes r\_n'\right)^{-1},\tag{10}$$

where *a* 0 *in* is the fuzzy comparison value between the pair criterion *i* and criterion *n*; and *r* 0 1 is the geometric mean of the fuzzy comparison values for criterion *i* compared to each of the other criteria; *w* 0 *i* is the fuzzy weighting of the *i*th criterion; and *w* 0 *i* = (*lw<sup>i</sup>* , *mw<sup>i</sup>* , *uw<sup>i</sup>* ), where *lw<sup>i</sup>* , *mw<sup>i</sup>* and *uw<sup>i</sup>* are the values of the lower, middle, and upper, fuzzy weighting of the *i*th 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]:

$$A' = \begin{pmatrix} a'\_{lj} \\ a'\_{lj} \end{pmatrix}\_{n \times n} = \begin{bmatrix} (1,1,1) & (l\_{12}, m\_{12}, u\_{12}) & \cdots & (l\_{1n}, m\_{1n}, u\_{1n}) \\ (l\_{21}, m\_{21}, u\_{21}) & (1,1,1) & \cdots & (l\_{2n}, m\_{2n}, u\_{2n}) \\ \vdots & \vdots & \ddots & \vdots \\ (l\_{n1'}, m\_{n1}, u\_{n1}) & (l\_{n2'}, m\_{n2'}, u\_{n2}) & \cdots & (1,1,1) \end{bmatrix} \tag{11}$$

where *a* 0 *ij* = (*lij*, *mij*, *uij*) and *a* <sup>0</sup>−1 *ij* = (1/*lij*, 1/*mij*, 1/*uij*) for *i, j* = 1, . . . , n and *i* 6= *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 *A* 0 in a first stage, as follows:

$$RS\_{\bar{i}} = \sum\_{j=1}^{n} a'\_{\bar{i}j} = \left(\sum\_{j=1}^{n} l\_{\bar{i}j}, \sum\_{j=1}^{n} m\_{\bar{i}j}, \sum\_{j=1}^{n} u\_{\bar{i}j}\right), \ i = 1, \ldots, n \tag{12}$$

Then, the value of the fuzzy synthetic extent in terms of the *i*th object is obtained through the normalization of the above relation, as follows [32]:

$$S'\_i = \sum\_{j}^{n} a'\_{ij} \otimes \left[ \sum\_{k=1}^{n} \sum\_{j=1}^{n} a'\_{kj} \right]^{-1} = \left( \frac{\sum\_{j=1}^{n} l\_{ij}}{\sum\_{k=1}^{n} \sum\_{j=1}^{n} u\_{kj}}, \frac{\sum\_{j=1}^{n} m\_{ij}}{\sum\_{k=1}^{n} \sum\_{j=1}^{n} m\_{kj}}, \frac{\sum\_{j=1}^{n} u\_{ij}}{\sum\_{k=1}^{n} \sum\_{j=1}^{n} l\_{kj}} \right), i = 1, \dots, n. \tag{13}$$

The computation of the degree of possibility of *S* 0 *<sup>i</sup>* ≥ *S* 0 *j* represents the third step and is achieved through Equation (14):

$$V(S\_{l}' \ge S\_{j}') = \begin{cases} \begin{array}{c} 1, \text{if } m\_{l} \ge m\_{j}, \\ \frac{(u\_{l} - m\_{i}) + (m\_{j} - l\_{j})}{(u\_{l} - m\_{i}) + (m\_{j} - l\_{j})}, l\_{j} \le u\_{i}, \; i, j = 1, \ldots, n; j \ne i \\\ 0, \text{ otherwise} \end{array} \tag{14}$$

where *S* 0 *<sup>i</sup>* = (*l<sup>i</sup>* , *m<sup>i</sup>* , *ui*) and *S* 0 *<sup>j</sup>* = *lj* , *m<sup>j</sup>* , *u<sup>j</sup>* . Considering that:

$$w'(a\_i) = \min\{V(S\_i' \ge S\_k')\}, k = 1, 2, \dots, n; k \ne i \tag{15}$$

the weight vector values can be calculated by:

$$w'(a\_i) = \begin{bmatrix} w'(a\_1), \ w'(a\_2), \dots, \ w'(a\_n) \end{bmatrix}^T. \tag{16}$$

The weight vectors were obtained using the following equation after a normalization process:

$$w(a\_1) = \begin{bmatrix} w(a\_1), \ w(a\_2), \dots, \ w(a\_n) \end{bmatrix}^T \tag{17}$$

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.

### *3.2. Data Used*

### 3.2.1. 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.

### 3.2.2. 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–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).

Topographic position index (TPI) is a predictor calculated from the DEM, which shows the relative position of a point in the research area in relation to the immediately neighboring regions [50]. The next TPI classes were established using the natural breaks method: −7.8–−1.8; −1.7–−0.5; −0.4–0.5; 0.6–1.9; 2–8.6 (Figure 3c). The following five classes of Topographic Wetness Index (TWI) were delimited using the natural break method: −4.4–4.7; 4.8–8.4; 8.5–11.8; 11.9–15; 15.1–23.1 (Figure 3d).

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

*Water* **2021**, *13*, x FOR PEER REVIEW 8 of 25

**Figure 2.** Flash-flood and flood predictors. (**a**) Slope; (**b**) Land use; (**c**) Hydrological Soil Group; (**d**) **Figure 2.** Flash-flood and flood predictors. (**a**) Slope; (**b**) Land use; (**c**) Hydrological Soil Group; (**d**) Convergence index.

Convergence index.

method: −7.8–−1.8; −1.7–−0.5; −0.4–0.5; 0.6–1.9; 2–8.6 (Figure 3c). The following five classes of Topographic Wetness Index (TWI) were delimited using the natural break method:

−4.4–4.7; 4.8–8.4; 8.5–11.8; 11.9–15; 15.1–23.1 (Figure 3d).

Topographic position index (TPI) is a predictor calculated from the DEM, which

**Figure 3.** Flash-flood and flood predictors. (**a**) Profile curvature; (**b**) Aspect; (**c**) Topographic position index (TPI); (**d**) Topographic Wetness Index TWI. **Figure 3.** Flash-flood and flood predictors. (**a**) Profile curvature; (**b**) Aspect; (**c**) Topographic position index (TPI); (**d**) Topographic Wetness Index TWI.

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

succeeded at a difference of 200 m (Figure 4a).

study area, the range from 763.1 m to 2202 m was split into seven classes that generally

**Figure 4.** Flash-flood and flood predictors. (**a**). Elevation; (**b**) Plan curvature; (**c**) Lithology; (**d**) Distance from rive. **Figure 4.** Flash-flood and flood predictors. (**a**). Elevation; (**b**) Plan curvature; (**c**) Lithology; (**d**) Distance from rive.

Plan curvature shows the difference between the surfaces on which the convergent and divergent runoff is manifested. Three classes were delimited for the plan curvature values (Figure 4b): −3–−0.1; 0–0.1; 0.2–1.9.

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.

### *3.3. Methodological Steps Implemented in the Present Study*

### 3.3.1. 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.

### 3.3.2. 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.

### 3.3.3. 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 flashflood 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.

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

3.3.4. Step 4: Computation the Flood Potential Index (FPI) Based on the Most Performant

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

*Water* **2021**, *13*, x FOR PEER REVIEW 12 of 25

FFPI Result

**Figure 5.** Flowchart of the present research. **Figure 5.** Flowchart of the present research.

#### **4. Results 4. Results**

#### *4.1. Bivariate Statistics Coefficients 4.1. 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. The values of bivariate statistics coefficients highlight the spatial relationships between the location of torrential areas and the classes/categories of flash-flood predictors. According to Table 1, the lowest weights of evidence coefficients were assigned to hydrological soil group D (−15.04), lithological category of sandy flysch, marls shale (−10.3), lithological category of clays, limestone (−9.27), and hydrological soil group C (−8.91). The highest weights of evidence values were attributed to slope angles higher than 45◦ (3.34), grassland land use (1.52), lithological category of conglomerates, breccias (1.48),

and slope angles between 15.1◦ and 25◦ (0.95). In terms of frequency ratio coefficients, the lowest values (0) were associated with agricultural zones, built-up areas, bare rocks, lithological categories of sandstone, gravels, clays, and limestone, zones with slope angles lower than 3◦ , and hydrological soil groups C and D. The highest frequency ratio values were assigned to zones with slope angles higher than 45◦ (13.5), grassland land use (2.07), lithological category of conglomerates, breccias (2.04), areas with slope angles between 15.1◦ and 25◦ (1.82), and convergence index class between −86 and −3 (1.78). In the case of SI coefficients, the lowest values were calculated for hydrological soil group D (−7.83), zones with slopes lower than 3◦ (−5.99), lithological category of sandstone and gravels (−5.25), lithological category of sandy flysch and marls shale (−5.12), and built-up areas (−4.94). The highest SI coefficients were obtained by zones with slope angles higher than 45◦ (2.6), grassland land use (0.73), lithological category of conglomerates, breccias (0.71), zones with slope angles between 15.1◦ and 25◦ (0.6), and convergence index class between −86 and −3 (0.58).


**Table 1.** Bivariate statistics of flash-floods conditioning factors classes.

### *4.2. 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:

$$\left[\sum\_{k=1}^{n}\sum\_{j=1}^{n}a\_{kj}^{\prime}\right]^{-1} = \left(88.48\ 130.16\ 182.17\right)^{-1} = \left(0.005\ 0.008\ 0.011\right)\tag{18}$$


**Table 2.** Fuzzy analytical hierarchy process evaluation matrix.


**Table 2.** *Cont.*

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:

$$w'(a\_i) = \{1.0.71\,0.32\,0.68\,0.68\,0.32\,0\,0\,0.32\,0.71\}^T \tag{19}$$

$$w(a\_i) = \{0.211\,0.15\,0.066\,0.143\,0.143\,0.066\,0\,0\,0.066\,0.15\}^T \tag{20}$$

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.

Flash-flood potential index values were mapped using the map algebra capability from ArcGIS software. All three flash-flood potential indices were standardized between 0 and 1 and then reclassified into five classes using the natural break method. In the case of FFPIFAHP-SI, very low values, situated from 0 to 0.25, were found in about 2.82% of the study area (Figure 6a). The values, ranging from 0.26 to 0.62, were mainly associated with the southern half and represent 28.9% of the entire river basin. The medium FFPIFAHP-SI class corresponded to approximately 13.86% of the Izvorul Dorului catchment. The high and very high potential were spread over a total of 54.43% of the entire catchment surface. The analysis of FFPIFAHP-FR revealed that the very low potential spanned 18.48% of the total study area and was present mainly in the southern half. The low flash-flood potential accounted for approximately 16.29% of the catchment surface, while the medium FFPIFAHP-WOE, with values from 0.32 to 0.48, covered 26.66% of the study zone. A zone including 38.57% of the research area was characterized by a high and very high flash-flood potential (Figure 6b). Following the application of the FAHP-WOE ensemble, only 0.62% of the Izvorul Dorului catchment had a very low flash-flood potential (Figure 6c). The low flash-flood potential, with values between 0.26 and 0.44, covered around 10.62% of

the entire territory, while the medium values quantified approximately 30.81% of the river basin. A percentage of 57.95% of the study area had high and very high FFPIFAHP-WOE values ranging from 0.64 to 1.

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


### *4.3. 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.

### *4.4. 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). *Water* **2021**, *13*, x FOR PEER REVIEW 17 of 25

**Figure 6.** Flash-flood potential index (FFPI) values. (**a**) Fuzzy Analytical Hierarchy Process—Statistical Index (FAHP-SI); (**b**) Fuzzy Analytical Hierarchy Process—Frequency Ratio (FAHP-FR); (**c**) Fuzzy Analytical Hierarchy Process—Weights of Evidence (FAHP-WOE). **Figure 6.** Flash-flood potential index (FFPI) values. (**a**) Fuzzy Analytical Hierarchy Process—Statistical Index (FAHP-SI); (**b**) Fuzzy Analytical Hierarchy Process—Frequency Ratio (FAHP-FR); (**c**) Fuzzy Analytical Hierarchy Process—Weights of Evidence (FAHP-WOE).

*4.3. 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

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

*Water* **2021**

*4.4. Flood Potential Index Computation* 

**5. Discussion** 

ent sources [57].

**Figure 7.** Receiver operating characteristic (ROC) Curves. (**a**) Success rate; (**b**) Prediction rate. **Figure 7.** Receiver operating characteristic (ROC) Curves. (**a**) Success rate; (**b**) Prediction rate. , *13*, x FOR PEER REVIEW 21 of 25

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

**Figure 8.** (**a**) River valleys torrential degree of the river; (**b**) Flood potential index classes. ogy proposed by Costache et al. [25], the valleys in the study area were identified and **Figure 8.** (**a**) River valleys torrential degree of the river; (**b**) Flood potential index classes.

In the last ten years, the study of the susceptibility to hydrological risk phenomena

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 sta-

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 methodol-

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 differ-

tistical bivariate methods and the fuzzy AHP model.



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.


**Table 5.** Properties of comparison matrices in the previous table.

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

### **5. 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 Trotus, 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.

### **6. 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 flash-flood 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.

**Author Contributions:** Conceptualization, R.C., Q.B.P. and A.B.; Data curation, R.C. and Q.B.P. Methodology, R.C. and A.B.; Writing—original draft R.C., Q.B.P. and A.B.; Writing—review and editing, R.C., Q.B.P. and A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a grant of the Romanian Ministry of Education and Research, CNCS-UEFISCDI, project number PN-III-P1-1.1-PD-2019-0424, within PNCDI III.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Alina Bărbulescu <sup>1</sup> , Cristina S, erban 2,\* and Marina-Larisa Indrecan <sup>2</sup>**


**Abstract:** This article proposes a new approach for determining the optimal parameter (β) in the Inverse Distance Weighted Method (IDW) for spatial interpolation of hydrological data series. This is based on a genetic algorithm (GA) and finds a unique β for the entire study region, while the classical one determines different βs for different interpolated series. The algorithm is proposed in four scenarios crossover/mutation: single-point/uniform, single-point/swap, two-point/uniform, and two-point swap. Its performances are evaluated on data series collected for 41 years at ten observation sites, in terms of mean absolute error (MAE) and mean standard error (MSE). The smallest errors are obtained in the two-point swap scenario. Comparisons of the results with those of the ordinary kriging (KG), classical IDW (with β = 2 and the optimum beta found by our algorithm), and the Optimized IDW with Particle Swarm Optimization (OIDW) for each study data series show that the present approach better performs in 70% (80%) cases.

**Keywords:** genetic algorithm (GA); IDW; spatial interpolation

### **1. Introduction**

Evaluating and predicting the effects of atmospheric factors dynamics, like precipitation and temperature, are of major importance for human activity, especially for zones with arid or rainy climates. Since water scarcity impacts billions of people worldwide, it is important to assess the water resources availability at ungauged locations [1]. Spatial interpolation methods are utilized for estimating the values of environmental variables using data recorded at neighbor locations. The most utilized approaches are classified as deterministic, geostatistical, and combined (or hybrid) [2,3]. The Inverse Distance Weighting (IDW) is a deterministic (mechanical) technique. The attribute values of any pair of points are related to each other, their similarity being inversely proportional to the distance between the two locations [4,5].

Since IDW does not involve advanced computational knowledge, researchers widely utilized it for spatial interpolation problems. Different authors presented comparable IDW performances with other spatial interpolation methods [6–11]. In [6,7], it is shown that IDW provided better or comparative results as ordinary kriging (OK) in the spatial interpolation of precipitation in Taiwan and Norfolk Island. Ly et al. [8] reported that OK and IDW provided the smallest root mean squared error in a study concerning the daily rainfall at the catchment scale in Belgium. Dong et al. [9] found that Ordinary CoKriging (OCK) performed better than OK and IDW when interpolating daily rainfall in a river basin from China. IDW, Thiessen Polygons Method (TPM), and kriging have been evaluated against the Most Probable Precipitation Method (MPPM) on annual, monthly, seasonal, and annual monthly maximum precipitation series from ten stations of 41 data [10]. IDW over performed TPM and OK, but underperformed MPPM. Chen et al. [11] proposed an improved regression-based scheme (PCRR) that was superior to IDW and multiple linear

**Citation:** B˘arbulescu, A.; S, erban, C.; Indrecan, M.-L. Computing the Beta Parameter in IDW Interpolation by Using a Genetic Algorithm. *Water* **2021**, *13*, 863. https://doi.org/ 10.3390/w13060863

Academic Editors: Giuseppe Pezzinga and Elias Dimitriou

Received: 14 January 2021 Accepted: 20 March 2021 Published: 22 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

regression (MLR) interpolation methods on data from the mesoscale catchment of the Fuhe River.

Even if the classical IDW (with the value of the parameter β = 2) was successfully employed for a long period for spatial interpolation problems, being easy to use, improving its performances was targeted by scientists. For example, Lu and Wong [4] proposed the weights' modification depending on the neighboring locations' distribution density around the unsampled place. Golkhatmi et al. [12] introduced altitude as a new variable in the IDW interpolation (keeping β = 2) and reported good results in the case study.

Another direction is finding the best β. This is an optimization problem by itself, targeted by many scientists [13–19]. For example, Noori et al. [13] employed IDW for estimating the distribution of precipitation in Iran, the value of the parameter (β) being recursively searched in the interval (1, 5], increasing its value each time. However, this grid-search procedure is time-consuming for small step sizes [5]. To avoid this drawback, Mei et al. [14] designed and implementing parallel adaptive inverse distance weighting (AIDW) interpolation algorithms by using the graphics processing unit (GPU) for accelerating the parameter finding. Gholipour et al. [15] propose a hybridization of IDW with a harmony search, which improves the convergence rate and reduces the search time.

In the same idea, hybrid methods have been proposed. Zhang et al. [16] combined Support Vector Machines (SVM) with IDW obtaining the SVM residual IDW, obtaining superior results by comparison to IDW and OK for the spatial interpolation of the multiyear average annual precipitation in the Three Gorges Region basin. Nourani et al. [17] used a two-stage framework for spatial interpolation of precipitation, employing, in the first stage, three artificial intelligence models that generate the input for the second stage, where they utilize IDW for spatial interpolation. Bărbulescu et al. [18] proposed a Particle Swarm Optimization approach (called OIDW) for finding a single β in IDW interpolation of maximum annual precipitation from the Dobrogea region (Romania). Chang et al. [19] applied a genetic algorithm (GA) to find the optimal distances between the gauged stations to minimize the estimation errors in IDW. Still, based on our knowledge, no attempt to optimize the choice of β parameter of IDW using a GA has been made so far.

On the other hand, GAs are widely used for solving real-life problems. For example, Ratnam et al. [20] improved seasonal air temperature forecasts using a genetic algorithm. Nasseri et al. [21] presented an optimized scenario for rainfall forecasting using a genetic algorithm coupled with an artificial neural network using rainfall hyetograph of recording rain gauges in the Upper Parramatta catchment (Sydney, Australia). Using the ability of GAs to search complex decision spaces, Sen and Ôztopal [22] utilized such an algorithm for optimizing the classification of rainy and non-rainy day occurrences using atmospheric data (temperature, humidity, dew point, vertical velocity). Heat conduction and control problems have also been solved by utilizing GAs [23,24].

In this context, this article proposes a new approach that optimizes the finding of the beta parameter of IDW. This is based on a genetic algorithm and finds a unique β for the entire study region, while the classical one determines different βs for different interpolated series. The algorithm is proposed in four scenarios crossover/mutation: singlepoint/uniform, single-point/swap, two-point/uniform, and two-point swap. Comparisons of its performances with those of the classical IDW (with β = 2 and the optimal beta found in our algorithm), ordinary kriging, and two versions of the optimized IDW by using Particle Swarm Optimization (OIDW) are also provided.

### **2. Methodology and Data Series**

### *2.1. IDW Interpolation*

The study problem is estimating a variable's values at ungagged locations employing the same variable's known values, registered at the neighboring observation sites [18]. In terms of mathematics, one can formulate the problem as follows. Given a set of spatial data of a variable z at *n* observation sites, *s*1, . . . , *s<sup>n</sup>* determine the same variable's values at the study site, *s*0.

The IDW interpolation formula is:

$$\widehat{z(s\_0)} = \sum\_{i=1}^{n} \frac{1/d(s\_0, s\_i)^{\beta}}{\sum\_{i=1}^{n} \left(1/d(s\_0, s\_i)^{\beta}\right)} \, z(s\_i), \,\,\beta > 1\tag{1}$$

where *z*[(*s*0) is the value computed for the site *s*0, *z*(*si*) is the value recorded at the site *s<sup>i</sup>* , *d*(*s*0, *si*) is the distance between *s*<sup>0</sup> and *s<sup>i</sup>* , and β is a parameter whose value is either given or determined by different optimization methods. In the original algorithm, β = 2 [25,26].

The interpolation quality depends on β which is generally determined after running a grid search. The time spent for finding the parameters is inverse proportional with the step of the grid.

### *2.2. Genetic Algorithms*

A genetic algorithm (GA) is a metaheuristic method inspired by natural selection laws that try to find optimal solutions to complex problems to which deterministic approaches usually cannot find a good result. The genetic operators, selection, crossover, and mutation establish a balance between the exploration and exploitation of the search space [27,28]. Exploration means that the algorithm searches for new solutions in new regions, while exploitation refers to making refinement to existing solutions to improve their quality. A function called fitness measures the quality of the solutions, which are represented by chromosomes.

A GA starts with a population of some random chromosomes and (by applying the principle of 'survival of the fittest') produces multiple generations by selecting in each one the fittest individuals for breeding. The mutation is then applied to increase the population diversity. Along with the generations, better individuals, i.e., better approximations to the solution, are obtained. The process continues until the fittest individual (the optimal solution) is found or the maximum number of generations is reached.

Using a genetic algorithm to solve a problem means finding the representation of the problem's solutions (encoding of the chromosomes), the fitness function, and the genetic operators. A chromosome is a feasible solution to the problem. In our case, a chromosome represents a real value of the parameter *a* ≤ β ≤ *b*. Thus, we apply a value encoding and get a binary string with the length *l*, calculated using the following formula (the default encoding of real values to binary strings):

$$2^{l-1} < (b-a) \* 10^z < 2^l \tag{2}$$

where *z* represents the given number of β's decimals. In this study, *l* = 9 bits.

The decimal value, *val*, of the binary chromosome representation, is computed by (3). We get the real value of a chromosome (β) by applying (4).

$$val = (\pounds - a) \times (2^l - 1) / (b - a) \tag{3}$$

or

$$\beta = a + val \times (b - a) / (2^l - 1) \tag{4}$$

The fitness function controls the possibility of individuals' reproduction. The better chromosome is (i.e., the better fitness is), the more likely it is to be selected for breeding the next generation. Since our goal is to minimize the error between the results obtained by the spatial interpolation and those recorded at the meteorological stations, the fitness function will record the mean standard error (MSE) between the known data and those computed by IDW. A GA performs best when a feasible solution maximizes the fitness function. Hence, we apply one of the most commonly adopted fitness mapping (inversion scaling), which does not alter the minimum location, but converts a minimization problem to an equivalent maximization one.

We use the mean standard error and mean absolute deviation (MAE) to evaluate the GA's performance. The lowest the MAE or MSE is, the better the algorithm performs.

Genetic operators are used for producing new generations of individuals with more diverse properties. There are three operators, selection, crossover, and mutation, which can set and, most of the time, find a good ratio between exploration and exploitation of the search space.

A selection operator determines the best individuals' regions that will exchange information to create a new generation. In this paper, the roulette wheel selection method [27] is used.

A crossover operator combines two or more parents to generate one or two offspring. It implements the idea that a swap of information between good individuals will generate an even better one. In this paper, the single-point crossover and the two-point crossover are used to create new offspring [29].

A mutation operator randomly modifies chromosomes with a given probability, *pm*, called mutation rate, leading to an increased population's structural diversity. Thus, a mutation operator facilitates the recovery of genetic material lost during the selection step and exploring new solutions. Here we used the uniform mutation [27] (one gene is randomly chosen and its value is modified) and swap mutation [27] (two positions on the chromosome are randomly selected, and their values are interchanged).

One may configure several control parameters in a genetic algorithm to achieve a balance between exploration and exploitation. If the population size is large, the search space is more explored than when the population size is small [28]. However, the runtime of the algorithm would increase. If crossover and mutation rates are high, the search will explore much of the solutions space, but there is a high chance of missing good solutions, the GA acting more like a random search. If crossover and mutation rates are low, the search space remains unexplored, and in this case, the GA resembles the hill-climbing algorithms. Therefore, we investigated the influence of the population size, crossover rate, mutation rate, and stop condition on the GA results. We performed each test ten times and averaged the results to increase their precision (as suggested in the literature). We implemented two crossover operators and two mutation operators to find the ones which are best suited for our problem. We also ran several tests for each pair of operators to see the relationship between the control parameters and the fitness value. Details are presented in the following sections.

### *2.3. New Approach for Estimating Beta*

The genetic algorithm we implemented is presented in the following.

**Input**: The distances between stations and the precipitation series recorded at these stations. **Output**: The optimal parameter value of β Begin

	- a. **Select** some chromosomes for crossover operation (the number of selected individuals is defined by the crossover rate)
	- b. **Apply** one of the crossover operators described in 2.2 to generate two new offspring
	- c. **Copy** the remaining chromosomes (that were not recombined) to the next generation

**End**

For a better understanding, a flowchart of the procedure of determination of the beta parameter is presented in Figure 1.

In order to find the best parameters settings for our problem, we fine tune our algorithm, which means creating several GA variants to test and find the best one, by slightly changing of GA parameters (population size, number of generations, crossover and mutation rates). We change only one parameter at a time, and try out several evolutionary literature-based test values. For example, for the crossover rate, the most used values in applications are in interval [0.6,1), whereas the mutation rate should be less than 10%. We run these GA variants on our problem, accept that parameter value at which the GA performs best, and continue to the next GA parameter, and so forth, to the last one. More precisely, to select the best population size, stop condition, and crossover rate the following steps are done.


**Figure 1.** The procedure flow chart. *nGen* represents the maximum number of generations.

Although the complexity of GAs has a probabilistic convergence time [30], the settings of our genetic algorithm are not complex, and, based on the experimental results that show that it converges in a short time, we may state that the best convergence time is logarithmic, *O*(*log*(*n*)), whereas the worst is linear, *O*(*n*). In the Results and Discussion, we present the recorded execution time (in seconds), which shows that the algorithms stop in short time for each test we did.

### *2.4. Data Series*

Dobrogea is a region covering a surface between the Romanian Littoral of the Black Sea, the lower Danube River, and the Danube Delta, situated in the southeast part of Romania and characterized by long droughts periods. Records show the absence of precipitation for 4–6 months per year after 1961, which affects agricultural activities. Researchers analyzed precipitation and temperature evolution in this zone, especially after 2010, to mitigate the drought effects [31–33].

The data we are working with is formed by the maximum annual precipitation series recorded during a period of 41 years at 10 main meteorological stations from the Dobrogea region (Figure 2).

**Figure 2.** Maximum annual precipitation series.

### **3. Results and Discussion**

Firstly, we ran several tests to find the settings of control parameters that are most likely to produce the best results. We started with predefined values from the literature [27,28] for the stop condition (10 generations), crossover rate (0.75), and mutation rate (0.015). Then, we varied the population size (from 10 to 80, with a step of 5) and computed the fitness function's corresponding values, run time, and β. Table 1 shows the relationship between the fitness value and population size.

For each pair of operators (single-point/uniform, single-point/swap, two-point/uniform, two-point/swap), we chose the optimal population size to be the lowest value from which population growth does not significantly influence the modification of the fitness value. This is 45, 40, 30, and 35 individuals, respectively, and the fitness function value is 0.0317. The corresponding β values obtained in the four scenarios are 1.256, 1.372, 1.308, and 1.336, respectively. These results are highlighted in Table 1.


**Table 1.** The impact of population size on the GA accuracy. The best results are highlighted.

Since we used a predefined number of generations as the stop condition, in the second stage, we had to determine its optimal value. To find it, for each pair of operators, we ran tests with several values of the number of generations, the population size previously estimated (45, 40, 30, 35, respectively—Table 1), keeping the mutation rate set to 0.015, and the crossover rate set to 0.75. Table 2 and Figure 3 show that the fitness value does not improve after a certain number of generations (which is the optimal number of generations).

**Table 2.** The impact of the number of generations on the GA performance. The best results are highlighted.


**Figure 3.** The impact of the number of generations on fitness value in (**a**) single-point/uniform, (**b**) single-point/swap, (**c**) two point/uniform, (**d**) two point /swap scenarios.

Based on the highest fitness value, the optimum number of generations determined was nine for the single-point/uniform mutation scenario, whereas, for the other scenarios, the number of generations was five. The corresponding β-values obtained are 1.228, 1.242, 1.122, and 1.362, respectively (highlighted in Table 2, together with the fitness and beta values). Since in our algorithm, β is represented by a chromosome, and several genetic operators have been used, different chromosomes (β) could produce the same fitness value.

The next step is the determination of the optimal crossover rate. For this aim, we ran the algorithm in each scenario with different values of the crossover rate (from 0.6 – value suggested in the literature to 0.95, with a step of 0.05), the population size and number of generations previously determined (45, 40, 30, and 35 individuals, respectively; 9, 5, 5, 5 generations, respectively), and the mutation rate kept to 0.015. We chose the optimal crossover rate for each pair of operators to be the value that gives the best (highest) fitness.

From Table 3 it results that the best crossover rates are 0.8 when using a singlepoint/uniform scenario, 0.6 for single-point/swap, 0.75 for two-point/uniform, and 0.7 for two-point/swap. These values correspond to the highlighted sequences of values in Table 3.

The last step was the determination of the best mutation rate. Therefore, we analyzed the impact the mutation rate has on the GA's results. We considered the population size, the number of generations, and the crossover rates we established in previous stages, and we performed new tests aiming at detecting the value of the mutation rate. For example, for single-point uniform mutation, we took the population size = 45, the number of generations = 9, the crossover rate = 0.80, and ran the tests for a mutation rate from 0.02 to 0.1, with a step size of 0.01. For each pair of operators, we search the optimal mutation rate for which the fitness value evolves to a maximum along with the generations. Table 4 contains the values of the fitness function obtained after running the algorithm in the four scenarios, with different mutation rates. For example, in Table 4a we present the values of the fitness function obtained for each generation (from 1 to 9) and the mutation rates from 0.02 to 0.1, in the single-point crossover/uniform mutation scenario. The highest fitness value is obtained after nine generations in the single-point crossover/uniform mutation, with a mutation rate of 0.06 (the sixth column–the highlighted values).


**Table 3.** The impact of crossover rate on the GA accuracy. The best results are highlighted.

**Table 4.** The impact of mutation rate on the GA accuracy. The best results are highlighted.



**Table 4.** *Cont.*

In the single-point crossover/swap mutation (Table 4b), the highest fitness value is 0.0315, obtained after 5 generations, with a mutation rate of 0.08 (the eighth column of Table 4b). For the two-point crossover/uniform mutation and two-point crossover/swap mutation (Table 4c,d), the best mutation rates are 0.04 and 0.05, respectively, and the corresponding value of the fitness function is 0.0313 (contained in the highlighted columns the fourth and the fifth, respectively).

After setting the optimal parameters, determined in the previous stages, we finally ran the algorithm to determine the optimum beta parameters. Table 5 summarizes the parameters used to implement the proposed genetic algorithm (columns 2–5), the fitness function obtained after running the algorithm with these parameters (column 6), the execution time (column 7), and the value obtained for the IDW's parameter (last column).

Remark that the values of β are different when using different scenarios, even if the fitness value is the same. This is due to the specifics of the individuals' selection and operations in GAs.

The lowest execution time (0.6188) is obtained when using the single-point/swap scenario and the highest one (10.5875 s) when using a two-point/swap procedure. Even if in the two-point/uniform case, the population size and the number of generations are the smallest, the execution time is high (the second-highest).



Table 6 contains the MSE and MAE for each station and the average (the last row of the table) computed after running the algorithm in each scenario. Comparing the MSEs in the two-point/swap and single-point/uniform (single-point/swap, and two-point/uniform) scenario, they are smaller in 70% (70%, 70%) cases, so our algorithm, in two-point/swap scenario, performs better in 70% cases compared to the other three scenarios. The MSEs' averages (31.5874, 31.5306, 31.5188, 31.5153) are comparable, the smallest being obtained in the two-point/swap scenario, followed by the third one.

Comparing the MAEs in the two-point/swap and single-point/uniform (singlepoint/swap, and two-point/uniform) scenario, they are smaller in 80% (80%, 80%) cases, so our algorithm, in two-point/swap scenario, performs better in 80% cases compared to the other three scenarios. The MAEs' averages (23.6352, 23.5542, 23.5308, 23.5228) are comparable, the smallest being obtained in the two-point/swap scenario, followed by the third one.

The corresponding values computed for beta in the best two cases are 1.042 (in twopoint/swap) and 1.064 (in two-point/uniform).


**Table 6.** MSEs and MAEs in GA. The best results are highlighted.

For comparison reasons, we performed the classical IDW, with β = 2 (the value used in most applications) and β = 1.042. The MAE and MSE values computed for each station are presented in Table 7.

Comparing the MSEs in the two-point/swap algorithm (Table 6, column 8) with those from the IDW with β = 2 (Table 7 column 2), they are smaller in 70% cases (the first four stations, the sixth, eighth, and ninth), so our algorithm performs better in 70% cases. Comparing the MSEs in the two-point/swap algorithm with those from the IDW with β = 1.042 (Table 7 column 4), they are smaller in 60% cases (the second, third, fourth, sixth, eighth, and ninth stations), so our algorithm better performs in 60% cases.

In terms of the average MSEs, that in the two-point/swap approach is smaller than those of the IDW (β = 2), IDW (β = 1.042), and slightly higher than in KG (Table 7, the last column). Still, our approach is preferable against KG since it is difficult to determine the kriging parameters, requiring special knowledge of geostatistics.


**Table 7.** MSE and MAE in the classical IDW for β = 2 and β = 1.042. MSE in ordinary kriging (KG).

\* Results from [18], Table 2; \*\* Results from [18], Table 1.

The MAEs in the two-point/swap algorithm are smaller than those from the IDW with β = 2, in 80% cases (all, but the first and sixth station), and comparable for the first station. The MAEs in the two-point/swap algorithm are smaller than those from the IDW (β = 1.042), they in 60% cases (all, but the first and sixth station), and comparable for the first station.

The average MAE in the two-point/swap approach (23.5228) is significantly smaller than those in the IDW with β = 2 (26.90), and IDW with β = 1.042 (26.08).

From the computational viewpoint, the highest computational time in our experiment was 10.5875 (in the two-point scenario), while in the grid search to estimate beta with 3 decimals takes 60 seconds for each series, so, a total of 60\*10 stations = 600 seconds, which is 56.67 times higher than in our approach.

The last two columns of in Table 7 contains the MSE in the optimized IDW, denoted by OIDW [18], in two scenarios, as described in [18]—with different beta, found using a Particle Swarm Optimization (PSO) approach (column 7), or with a single best beta found by the same approach.

In term of MSE, our GA algorithm (Table 7, column 8) performs better that OIDW (Table 7, columns 7 and 8) in 80% of cases. Therefore, we can say that a significant improvement of the interpolation performances are obtained, that may reflect in the water management policy.

### **4. Conclusions**

In this article, we presented a new approach to finding the beta parameter in IDW, using a GA implemented in four scenarios. The settings of this GA were optimized for finding the best fitness function and, by consequence, the best parameter beta, for all the study sites, not only for some of them.

It is shown that the algorithm proposed here performs better (in all scenarios) than the classical one (with β = 2 and β = 1.042) in terms of average MSE and MAE. When compared the MSEs and MAEs for the individual stations, the following results are obtained:


The algorithm performs faster than the classical IDW, for which the running time on the same problem is 60s for each interpolated data series (so 600s for all ten series). It is easy to be implemented and used and can be applied to similar problems only by changing the input data.

Compared with other artificial intelligence methods used for finding beta (OIDW) our approach shows superior performances in 80% of cases.

Another advantage is that our algorithm provides a single beta for all the stations, optimizing the interpolation.

The results obtained in all four GA's scenarios are comparable. Since the execution time is the highest in the best scenario (Table 5), the other alternatives can be successfully used for the spatial interpolation when the number of series or the number of records per station is very high.

**Author Contributions:** Conceptualization, A.B. and C.S, .; methodology, A.B. and C.S, .; software, C. S, . and M.-L.I.; validation, A.B. and C.S, .; writing—review and editing, A.B., C.S, . and M.-L.I., visualization, C.S, . and M.-L.I.; supervision, A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-6358-9