Using GIS, Remote Sensing, and Machine Learning to Highlight the Correlation between the Land-Use/Land-Cover Changes and Flash-Flood Potential

Abstract

The aim of the present study was to explore the correlation between the land-use/land cover change and the flash-flood potential changes in Zăbala catchment (Romania) between 1989 and 2019. In this regard, the efficiency of GIS, remote sensing and machine learning techniques in detecting spatial patterns of the relationship between the two variables was tested. The paper elaborated upon an answer to the increase in flash flooding frequency across the study area and across the earth due to the occurred land-use/land-cover changes, as well as due to the present climate change, which determined the multiplication of extreme meteorological phenomena. In order to reach the above-mentioned purpose, two land-uses/land-covers (for 1989 and 2019) were obtained using Landsat image processing and were included in a relative evolution indicator (total relative difference-synthetic dynamic land-use index), aggregated at a grid-cell level of 1 km². The assessment of runoff potential was made with a multilayer perceptron (MLP) neural network, which was trained for 1989 and 2019 with the help of 10 flash-flood predictors, 127 flash-flood locations, and 127 non-flash-flood locations. For the year 1989, the high and very high surface runoff potential covered around 34% of the study area, while for 2019, the same values accounted for approximately 46%. The MLP models performed very well, the area under curve (AUC) values being higher than 0.837. Finally, the land-use/land-cover change indicator, as well as the relative evolution of the flash flood potential index, was included in a geographically weighted regression (GWR). The results of the GWR highlights that high values of the Pearson coefficient (r) occupied around 17.4% of the study area. Therefore, in these areas of the Zăbala river catchment, the land-use/land-cover changes were highly correlated with the changes that occurred in flash-flood potential.

1. Introduction

The increase in frequency and magnitude of hydrological hazards, such as flash floods, is due to the last decades of climate change at a global scale [1], as well as to the changes affecting land-use/land-cover and land management. The flash floods, generated by the surface runoff on the slopes, represent one of the most dangerous natural hazards producing the greatest damage to human communities. Consequently, it is essential to study and monitor the areas holding a high potential for surface runoff.

In the scientific literature, the topic of runoff potential and flash flooding was approached by different authors in their studies [2,3,4,5,6,7,8,9,10,11,12]. At the same time, numerous studies focused on assessing the connection between land-use changes and different features of hydrological hazards, such as the multiplication in number of floods due to the conversion of natural vegetation cover into farmland, increase in floods volume along with built-up space extension, enhancement of runoff through the loss of forested areas, or decrease in frequency of peak discharge by increasing the afforestation [13,14].

The variety of research directions is connected to numerous factors contributing to land-use change, such as urban extension, intensification of agriculture, deforestation, afforestation, land abandonment, and the like [15,16,17]. The decisive role played by these land-use and land-cover changes on the enhancement of hydrological hazards resides in the derived alteration of quantitative relationships between the water cycle elements, such as interception, infiltration, or evaporation [18,19]. The vegetation cover is highly involved in this context by influencing the evapotranspiration [20,21,22], and subsequently, the soil water balance.

The geoinformation technologies, namely geographic information systems (GIS) and remote sensing, represent useful tools for assessing the land-use changes, as well as flash-flood potential [23]. The continuous advancements in the field of remote sensing have enabled researchers to obtain satellite images with short revisit cycles and different spatial resolutions, depending on the sensor, for any part of the world. Specifically, Landsat sensors have a sufficient spatial resolution that is detailed enough to characterize the processes that influence the land-use in the study area. The role of GIS is to further process the remotely sensed data as well as other source data to create the vector or raster inputs (i.e., flash-flood inventory and flash-flood predictors) for the selected modelling approaches.

Different methods or models can be used for assessing the flash-flood potential/susceptibility.

Multi-criteria-based techniques, such as the analytical hierarchy process [24,25], analytical network process [26]; weighted linear combination [27]; and Viekriterijumsko kompromisno rangiranje (VIKOR), Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS), or Simple Additive Weighting (SAW) techniques [28,29], represent a group of methods that have been used in a number of flood susceptibility studies. However, the disadvantage of these methods is that they usually rely on subjective evaluations when weighting factors. For that reason, many researchers preferred to use more objective methods that are based on statistics, including frequency ratios, logistic regression, or other bivariate or multivariate statistical models [30,31,32,33].

In recent years, however, more sophisticated data-driven methods have been used for flash-flood susceptibility assessment because they are more robust and have a better capability to handle complicated relationships between input variables. In particular, different machine learning algorithms [34] have been used, including artificial neural networks (ANN) [35], adaptive neuro-fuzzy inference systems (ANFIS) [36], decision trees (DT) [37], or support vector machines (SVM) [38]. Furthermore, the current trend in flood susceptibility mapping is to create hybrid or ensemble machine learning models in order to achieve an even better accuracy of the resulting model [39,40,41].

The multi-layer perceptron neural network (MLP), which was used in this study, belongs to the artificial neural network techniques, having one of the highest prediction powers and efficiency in modelling [42,43]. This technique is capable of universal modelling, but mainly it is used for analyzing nonlinear, multivariate, and complex processes in the real world [44]. The MLP neural network has been used for different types of natural hazards, such as landslide assessment [45,46], land subsidence [47], or for determining the soil consolidation and compression coefficient [44,48]. Despite the advantages of this method, its application in flash-flood or flood susceptibility studies is limited to only few studies [49,50,51,52].

Although the abovementioned methods and models reach a high predictive capability, it is still difficult to accurately predict flash floods, as well as to provide an answer to the increase in flash floods frequency and its connection to the landscape properties, such as land-use.

Therefore, the aim of this study is to investigate the correlation between the land-use change and the flash-flood potential changes in the Zăbala catchment (Romania) between the years 1989 and 2019 using GIS, remote sensing, and machine learning algorithms. The procedural steps of the applied workflow are the following:

(i) A bi-temporal analysis of two satellite scenes, Landsat 5 for 1989 and Landsat 8 for 2019, using the supervised classifications in order to derive the land covers. These two years were selected in order to reflect the land-use changes before and after the change in political regime of the country.

(ii) The land-use change analysis where land-use/land-cover changes were summarized in the Markov matrix and subsequently quantized by the total relative difference – synthetic dynamic land-use index (TRDSDLUI).

(iii) Processing the flash-flood inventory, consisting of 127 flash-flood locations and 127 non-flash-flood locations, and 10 flash-flood conditioning factors.

(iv) Application of the multilayer perceptron (MLP) neural network model for deriving the flash-flood potential index (FFPI) for both years (1989 and 2019).

(v) Validation of the MLP model using 30% of flash-flood locations (testing dataset).

(vi) Calculation of the changes in flash-flood potential between 1989 and 2019 using the relative difference for flash-flood potential index (RDFFPI).

(vii) Statistical correlation between TRDSDLUI and RDFFPI using a geographically weighted regression (GWR).

From the methodological point of view, the importance of this study can be seen in an integrated approach of the connection between land-use/land-cover changes and flash-flood potential.

2. Study Area

Zăbala river catchment is located in the central-eastern part of Romania (Figure 1), measuring an area of approximately 600 km². The altitude values range from 312 m in the sub-Carpathian region, at the confluence with Putna river, to 1700 m (Figure 1) in the Curvature Carpathians zone (Vrancei Mountains), where the largest area of the basin is located. The lithology, which is represented by the internal flysch of the Carpathians, alternating between hard and impermeable rocks (such as sandstones, tuffs, or shales), favors the surface runoff [53,54].

The mean relief slope angle in the study area, calculated through the geoprocessing of the digital elevation model, is equal to 12.7°. The areas recording high slope angles (above 15°), where runoff is highly increased, cover 31% of the total catchment area, demonstrating the exposure of the study area at risk for phenomena associated with surface runoff. According to the results of the supervised classification of the Landsat 8 scene from 2019, the forest vegetation covers approximately 60% of the basin total area, forming extended and compact patches in the mountainous area. In the sub-Carpathian zone, the absence of forest vegetation and the dominance of pastures overlapping a loamy-clay soil texture increase the runoff potential.

The high slope angles, combined with the hard lithological substrate and the presence of pastures covering a soil with a fine texture determine the augmentation of the flash-flood susceptibility within the perimeters of human communities from the sub-Carpathian zone. According to the General Inspectorate for Emergency Situations of Romania (GIES), the following localities were frequently affected by flash floods: Nistorești, Nereju, Paltin, and Spulber (Figure 1). The destructive consequences of flash floods have occurred along the Zăbala river during the last few years (20.10.2009, 11.05.2010, 22.04.2016, 28.05.2017, and 05.04.2019), leading to important damages to the road network and households, as well as the loss of human lives [55].

3. Data and Methods

Due to the fact that the present study is based on the spatial analysis of land-use/land-cover changes in relation to the changes produced within the flash-flood potential across the study area between 1989 and 2019, the vast majority of the data used are geospatial. Thus, the determination of land-use/land-cover changes was based on the use of two Landsat 5 and Landsat 8 satellite scenes. The scenes from Landsat 5, acquired on 18 August 1989, and Landsat 8, acquired on 17 May 2019, were downloaded from the United States Geological Survey (USGS ) Earth Explorer platform and were used in the supervised classification procedure. It should be mentioned that the spatial resolution of Landsat scenes was 30 × 30 m. Along with the Landsat imagery, other databases were used to assess the flash-flood potential across the study area. Thus, six flash-flood predictors are morphometric factors that were derived based on the digital elevation model (DEM). The DEM for the study area was created from contour lines digitized at a 5-m equidistance based on the Topographical Map of Romania [56]. The DEM was derived using the same spatial resolution as the Landsat scenes. Another two flash-flood predictors, represented using hydrological soil group and lithology, were extracted from the Digital Soil Map of Romania, 1:200,000 [57], and the Geological Map of Romania, 1:200,000 [58]. The monthly precipitation amounts from 1961 to 1989 and 1990 to 2018, collected from 32 meteorological stations around the study area, were taken into account along with the DEM, in order to calculate the spatial variability of the flash-flood predictor represented by the modified Fournier index.

In order to analyze the relationship between land-use/land-cover changes and flash flood potential during the period 1989–2019, the collected data was inserted into the methodological workflow (Figure 2), which includes four main steps: remote sensing image processing, land-use/land-cover change analysis, flash-flood potential assessment, and geo-statistical analysis. These steps are presented in the following sub-sections.

3.1. Remote Sensing Images Processing

The two land uses/land covers, for 1989 and 2019, were extracted from remote sensing imagery. This method of data acquisition was employed in numerous studies, focusing on highlighting the effects of land-use changes on different phenomena [59,60,61,62]. In the literature, there are many different methods used to extract the land cover from remotely sensed imagery. Thus, the unsupervised classification is one of the methods used to extract the land uses/land covers. This method can achieve very good results when the values belonging to a certain land-use/land-cover class are very close in terms of the spatial distribution [63]. Regarding the supervised classification, one of the most known methods is modified k-nearest neighbor [64], which allows the users to validate the value of each pixel according to its neighbors. Before the development of the PCs processing power, the minimum-distance-to-means classification method was considered faster than other methods because it is based on a simple mathematical algorithm [65]. The advantage of the support vector machine classification method lies in the fact that it allows for the analysis of high dimensional datasets [66]. Maximum likelihood is another very popular classification method that is based on a statistical rule that analyses the probability of each pixel belonging to a particular class [67]. In a study carried out on Istanbul city, Erbek et al. [68] showed that the maximum likelihood algorithm can obtain a higher classification accuracy than other methods (learner vector quantization) or a lower classification accuracy (than multilayer perceptron) depending on different factors. Nevertheless, given the strengths of the maximum likelihood algorithm, such as the advantages of both the mean vectors and the multivariate spreads of each class and the identification of those elongated classes [69], we decided to also use the maximum likelihood in the present analysis.

The primary source for land-use/land-cover classification was the two aforementioned Landsat scenes (Figure 3a,b) (spatial resolution 30 × 30m), acquired on 18 August 1989 and 17 May 2019. The scenes, provided by USGS Earth Explorer, were pre-processed and then classified using ENVI 5.0 software, made by L3HARRIS Geospatial, Broomfield, Colorado – United States of America. Only the reflective bands (1–5 and 7) of the sensor were used.

In a first stage, in the image pre-processing procedure, the radiometric and relative atmospheric corrections were performed in order to eliminate the effects of multiple factors, such as the changes in sensor characteristics, atmospheric condition, solar angle, and sensor view angle [70,71,72,73]. In terms of radiometric corrections, it was mandatory to convert the digital number of pixels into spectral radiance values, and further, the radiance into reflectance [74,75,76]. Furthermore, the dark object subtraction method was used for the relative atmospheric correction [77,78].

After pre-processing, the supervised maximum likelihood image classification was performed [79]. The first step included the identification of many ground truth regions of interest belonging to seven land-use/land-cover categories detected in Landsat imagery. In this respect, we tried to extract, based on the satellite image, the best training and testing areas to extract the seven land-use/land-cover datasets. We correlated the information found in the image with different land-use/land-cover information available online (we used the official Corine Land Cover dataset [80]), and we assigned the different reflectance values of the satellite imagery to the land-use/land-cover class that they overlaid/corresponded with. In this way, different reflectance values that corresponded with the same land use/land cover were assigned correctly in a more accurate final classification. In accordance with Makantasis et al. [81], the regions of interest were divided into training areas (80%) and testing areas (20%). The partitioning into training and testing areas were based on a random selection of the pixels within region of interest.

The training areas were selected based on spectral responses from various combinations of spectral bands (3,2,1; 4,3,2; and 5,4,3 for Landsat 5, and 4,3,2; 5,4,3; and 6,5,4 for Landsat 8). The training sets were created by digitizing polygons on the image. Then, the image processing software system performed a signature analysis that involved a statistical characterization of the training areas. Once a statistical characterization was achieved for each information class, the image was classified by examining the reflectance for each pixel and by making a decision about which of the signatures it resembled most [82].

The post-processing operations included exporting the classified data as a vector format (polygon) and further analysis was performed in ArcGIS for Desktop10.5, developed by Environmental Systems Research Institute (ESRI), Redlands, CA, United States. In order to test the accuracy of the classifications, which is a crucial step in the aerial imagery classification process [83], two confusion matrices, one for each year using ground truth testing areas, were constructed. Based on the confusion matrices, the kappa index was also determined to estimate the accuracy. It should be noted that we decided to use a “standard” post classification comparison method (using a confusion matrix), and the changes obtained were mostly visible for areas converted from forests to pastures. The resulting land uses/land covers were introduced using two analyses: calculation of the land-use/land-cover change indicator (${\mathrm{TRD}}_{\mathrm{SDLUI}}$) and the FFPI for 1989 and 2019.

3.2. Land-Use/Land-Cover Change Analysis

The land-use changes were summarized in the Markov matrix, which was created according to the following steps. First, the raster datasets for 1989 and 2019 land uses/land covers were reclassified using codes between 10 and 80 for the year 1989 and between 1 and 8 for the year 2019; second, the transitions were obtained in ArcGIS for Desktop 10.5 software via the addtion of the two raster datasets, through the Raster Calculator tool.

At the same time, the changes were quantized and spatialized by means of ${\mathrm{TRD}}_{\mathrm{SDLUI}}$ (Equation (1)), derived from an annual ratio for land-use/land-cover change, namely the synthetic dynamic land use index ${\mathrm{TRD}}_{\mathrm{SDLUI}}$ [70]:

$${\mathrm{TRD}}_{\mathrm{SDLUI}}=\frac{{}_{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left|{\Delta \mathrm{LU}}_{\mathrm{i}-\mathrm{j}}\right|}}{2{\displaystyle \sum }{\mathrm{LU}}_{\mathrm{i}}}\times 100(\%)$$

(1)

where: ${\mathrm{TRD}}_{\mathrm{SDLUI}}$—total relative difference derived from synthetic dynamic land use index;

${\mathrm{LU}}_{\mathrm{i}}$—the area of the land use/land cover i at the starting date (1989);

${\Delta \mathrm{LU}}_{\mathrm{i}-\mathrm{j}}$—transition area from land use/land cover i to other land uses/land covers during the studied period (1989–2019).

This index measures the intensity of land-use/land-cover changes in each spatial unit by summarizing all transitions recorded among different land uses/land covers. The index was spatially computed by means of a grid overlapping the study area, having the cell size of 1 km² and data geo-processing for calculating this indicator required a workflow [84] including elementary spatial analysis in ArcGIS for Desktop 10.5 software.

3.3. Flash-Flood Potential Assessment

The flash-flood potential index (FFPI) was calculated for the first time by Smith [85], and has subsequently been used and adapted by numerous researchers [5,53,85,86]. Generally, the FFPI is defined as a dimensionless indicator whose values can be obtained by overlapping many geographical factors that influence the surface runoff process [87]. The geographical factors that were widely used in the aforementioned studies for FFPI computation are: slope angle, soil, land use/land cover, lithology, profile curvature, and convergence index. The FFPI values are easily obtained in GIS software by using Map Algebra [53,87,88]. For the present study, the calculation and spatialization FFPI for the Zăbala catchment (for 1989 and 2019) were performed by using the multilayer perceptron model combined with GIS techniques.

3.3.1. Flash-Flood Inventory

In order to generate the flash-flood potential index (FFPI) for the two reference years, an essential step was the inventory of the previous flash-flood events that took place in the study area. It should be mentioned that the location of historical flash-flood events was based on data provided by the General Inspectorate for Emergency Situations (GIES) of Romania. Thus, two data sets were established, one for 1989 and the other for 2019. For the year 1989, the phenomena produced between 1965 and 1989 were taken into account, while for the year 2019, the events produced between 1990 and 2019 were considered. It should be mentioned that the year 1965 was chosen because this year was the first in which the flash-flood events were quantified by the national authorities. The year of 1990 was chosen because it is the first year of the period 1990–2019 being also the one immediately after the first reference period 1965–1989, namely the collectivization period. Within the first period (1965–1989) the most destructive flash-flood events occurred in 1970, while within the second period (1990–2019), the most destructive flash floods were recorded in 2005. Due to this fact, the vast majority of flash-flood locations belong to these two years. For a greater objectivity of the present study, from the total number of flash floods collected, an equal number of flash-flood locations for each of the two analyzed periods were selected. Thus, 127 flash-flood locations were quantified (Figure 1), which were subsequently divided into training (70%) and validating (30%) datasets [89]. Further, in order to ensure a high accuracy of the results [90], 127 non-flash-flood locations were generated inside of areas where the flash floods were absent and in which there was a very low probability of flash flood occurrence. These were located in the forested areas with very low slopes in general, and were divided into training (70%) and validating (30%) datasets [39]. It should be mentioned that for a higher objectivity of the results, the same dataset regarding the absence of the phenomena was used for both periods.

3.3.2. Flash-Flood Conditioning Factors

Further, in order to calculate the FFPI in the GIS environment, the following 10 geographical factors that influence the surface runoff were derived: slope angle (Figure 4a), profile curvature (Figure 4b), hydrological soil group [91] (Figure 4c), lithology (Figure 4d), topographic position index (TPI) (Figure 4e), topographic wetness index (TWI) (Figure 4f), convergence index (Figure 5a), aspect (Figure 5b), modified Fournier index for 1989 and 2018 [92] (Figure 5c,d respectively), and land use/land cover for 1989 and 2019 (Figure 7a,b, respectively). The first two morphometric factors, namely the slope and the profile curvature, were derived from the digital elevation model (DEM). The relief slope is the morphometric factor that highlights the gravitation force influencing the surface runoff [86,87,88,93]. Consequently, the runoff intensifies as the geo-declivity value increases. The profile curvature is another morphometric factor that differentiates convex areas (having negative values), affected by the accelerated runoff form the concave ones (having positive values), facing a decelerated runoff [94]. The hydrological soil groups and lithology were converted from the vector format into raster datasets with a 30-m spatial resolution. As was mentioned before, the spatial distribution of hydrological soil groups was extracted from the Romanian Soils Map (2002) [95]. This property of the soil has a strong influence on the surface runoff [96] as it directs water infiltration into the soil. Therefore, a fine clay texture, which includes the hydrological group D, is supposed to favor the water runoff due to the reduction of infiltration, while a sandy texture, which includes the hydrological group A, will create an increase in the infiltration ratio and a decrease in runoff potential. The classification of soils into four groups are made based on their hydraulic conductivity, which is influenced by the texture [96]. The soil group A is characterized by a hydraulic conductivity higher than 40 μm/s, while the soil group B, with a loamy texture, has a hydraulic conductivity between 10 and 40 μm/s. The soil group C, with a hydraulic conductivity between 1 and 10 μm/s, includes soils with a loamy-clay texture. The soil group D has a hydraulic conductivity lower than 1 μm/s [94]. The lithological layer for the Zăbala catchment was extracted from the 1:200,000-scale Geological Map of Romania. TPI (Figure 4e) is a key factor that shows the difference between the altitude of a raster cell and the altitude of the neighboring cells [97]. The TPI map was constructed by reclassifying its values into five classes taking into account natural breaks intervals. TWI (Figure 4f) is a morphometric factor that highlights the areas from the ground that are most favorable to water flow accumulation [98]. Like in the case of TPI, the TWI map was created by dividing its values into five intervals. The convergence index (Figure 5a), derived from the DEM, is a widely used indicator within studies that approach the flash-flood potential [99]. The map of the convergence index was obtained after the reclassification of its values into five classes according to the literature [50]. Aspect (Figure 5b) is an important indicator mainly for soil moisture status. It was taken into account in this study because it has an indirect influence on flash-flood phenomena [97]. The amount of rainfall also plays an important role as its increase determines the rise of rapid runoff potential on the slopes. One of the most important characteristics of the rainfall that influences the flash-flood occurrence is the intensity. In order to assess the spatial variability of the rainfall intensity across the Zăbala catchment, the modified Fournier index was used. This index is widely used to estimate the rainfall intensity [97]. The modified Fournier index was calculated according to Equation (2):

$$MFI={\displaystyle \sum }_{i=1}^{12}\frac{P_i^2}{P}$$

(2)

where P_i is the mean monthly precipitation (mm), and P is the mean annual precipitation (mm).

In order to determine the MFI values for the two years considered for the analysis, the monthly precipitation amounts from 1961 to 1989 and 1990 to 2018 were taken into account. According to the information already presented in the manuscript, the precipitation datasets were collected from a total number of 32 meteorological stations located around the study area. In the first step, the MFI values were calculated for each individual meteorological station using Equation (2). In order to represent the spatial variability of MFI values across the study area, the residual Kriging method was used [100]. Thus, the first stage of the application of the residual Kriging method was the exploration of the relationship between the MFI values and the altitude of the meteorological stations. The analysis of this relation was performed through ordinary least squares (OLS) regression [101]. In the applied OLS regression model, the independent variable is the altitude of meteorological stations, while the dependent variable is the MFI value [96]. Finally, the OLS equation, in which the MFI values (Y) will vary according to the altitude (X), were derived for each of the two periods. These OLS equations were further used in the Map Algebra (Raster Calculator tool) function of ArcGIS 10.5, where the X parameter was replaced with the DEM and the output of this procedure was the theoretical values of MFI in a raster format. The residual of the OLS regression was also calculated for each point represented by the meteorological stations and then were interpolated at a 30 × 30-m spatial resolution. The final estimated MFI values across the study area were obtained through the addition of the theoretical MFI values in raster format with the rasters resulting from the interpolation of the residual values [102] (Figure 5c,d). Once obtained, the MFI values were divided into four classes, taking into account the existing classification in the international literature [96].

The land uses/land covers for 1989 and 2019 (Figure 7a,b, respectively) were acquired via the classification of the Landsat images (USGS). This last factor represents a major driving force for the surface runoff as it influences different processes, such as evapotranspiration or interception of the pluvial water, as well as producing a specific roughness coefficient. Water runoff will be diminished above the forest vegetation land use/land cover, while on areas having a low roughness coefficient [103], such as pastures, bare rocks, or built-up areas, the runoff will increase.

3.3.3. Description and Configuration of Multilayer Perceptron (MLP) Model for FFPI Computation

In order to avoid the noisy data within the machine learning models and also in order to reduce the redundant information used in the training process, the predictive ability of the flash-flood predictors will be assessed. The predictive ability was estimated through the linear support vector machine, implemented in Weka 3.9 software developed by the Univerisity of Wakaito, Hamilton, New Zealand.

The multilayer perceptron (MLP) is a type of neural network that contains many interconnected nodes whose main aim is to solve the complex problems regarding the spatial relationships between influencing variables and the presence of a phenomenon. In the present case study, the influencing variables were flash-flood predictors, while the flash-flood location represented the dependent variable. A detailed description of the MLP background was carried out by Costache and Tien Bui [86] and Costache et al. [93]. The MLP architecture is composed of three main elements: (i) the input layer that contains the input neurons represented by flash-flood conditioning factors; (ii) the hidden layer, which has a various number of hidden neurons that have a crucial role in the MLP training procedure; and (iii) the output layer, which in the present case study, contained two neurons represented by flash-flood locations and non-flash-flood locations. An extremely important role in the MLP algorithm is held by the number of hidden neurons within the hidden layer. In this case, the number of hidden neurons was established based on the following formula [85,104]: 2 × N + 1, where N is the number of flash-flood predictors. Thus, 21 hidden neurons were included in the hidden layer (Figure 6).

The training sample of flood and non-flood locations for 1965–1989 and 1990–2019, together with the 10 flash-flood predictors, was used in order to train the MLP model. It should be mentioned that in the case of flash-flood predictors used in the MLP, the classes/categories were assigned the frequency ratio values calculated using Equation (3) [105]:

$$FR=\left(\frac{N\left(FXi\right)}{{{\displaystyle \sum }}_{i=1}^{m}N\left(FXi\right)}\right)/\left(\frac{N\left(Xj\right)}{{{\displaystyle \sum }}_{j=1}^{n}N\left(Xj\right)}\right)$$

(3)

where FR is the value of frequency ratio assigned to a class/category i of parameter j, N(FXi) is the number of flash-flood locations in a class i of a variable X, N(Xj) is the number of pixels inside a variable Xj, m is the number of classes contained by a flash-flood predictor Xi, and n is the number of flood-influencing factors inside the research zone.

Furthermore, the FR coefficients, obtained for the two periods, were normalized between 0.1 and 0.9 using the Equation (4):

$$v=\frac{\left(a-\min \left(r\right)\right)\times \left(\max \left(l\right)-\min \left(l\right)\right)}{\max \left(r\right)-min\left(r\right)}+min\left(l\right)$$

(4)

where v represents the standardized value of a, a represents the current value of the variable, r represents the limits of the range value, and l represents the limits of the standardization range.

The normalized values of FR were used as input data in MLP algorithm for determining FFPI₁₉₈₉ and FFPI₂₀₁₉. Thus, the MLP model was trained by using a maximum of 500 epochs and 30 validation thresholds. The values of the root mean square error (RMSE) for both years were also determined in order to assess the model efficiency. By applying the MLP algorithm using Weka 3.9 software, the importance of each flash-flood predictor was calculated. By multiplying the values of the MLP importance with the FR coefficients in ArcGIS 10.5, the FFPI₁₉₈₉ and FFPI₂₀₁₉ were derived. In order to compare the values and the changes produced between the years 1989 and 2019, the values of the two indices were normalized between 0.1 and 0.9 (Equation (4)) and were then grouped into the same five intervals.

3.3.4. FFPI Results Validation

The assessment of the results reliability is a mandatory stage for their use in the subsequent steps performed in the present analysis. One of the most used method for results validation is the ROC curve with its statistical indicator, the area under curve (AUC). A detailed explanation regarding the manner in which the ROC curve is constructed and how AUC values can be calculated has been done before by Costache and Tien Bui [85] and Costache [105].

3.3.5. FFPI Differences

The next stage of the study was the calculation of FFPI differences between 2019 and 1989. In order to correlate the variation of FFPI between 1989 and 2019 with land-use/land-cover change for the same period, the mean values of FFPI for both years were calculated at the same grid-cell level of 1 km², which was also used for ${\mathrm{TRD}}_{\mathrm{SDLUI}}$. This procedure was performed by using the Zonal Statistics form of the Spatial Analyst extension of ArcGIS for Desktop 10.5. Subsequently, FFPI values calculated for 1989 and 2019 were integrated into a relative evolution, according to Equation (5):

$${\mathrm{RD}}_{\mathrm{FFPI}\left(1989-2019\right)}=\frac{|{\mathrm{FFPI}}_2-{\mathrm{FFPI}}_1|}{{\mathrm{FFPI}}_1}\times 100(\%)$$

(5)

where ${\mathrm{RD}}_{\mathrm{FFPI}}$—relative difference for the flash flood potential index

${\mathrm{FFPI}}_1$—flash flood potential index at the starting date (1989)

${\mathrm{FFPI}}_2$—flash flood potential index at the end date (2019)

3.4. Geo-Statistical Analysis

Due to the complexity of our hypothesis and the heterogeneity of the studied area, the statistical correlation between ${\mathrm{TRD}}_{\mathrm{SDLUI}}$ and ${\mathrm{RD}}_{\mathrm{FFPI}}$ indicators was explored using a geographically weighted regression [106]. This method was chosen instead of using ordinary least squares (OLS), in spite of the reduced number of variables (only two), since the value of the Akaike information criterion [107] for GWR was inferior to OLS, suggesting a better suitability for the GWR [108]. When performing the GWR analysis, a moving window including 30 neighbours was applied, with the correlation coefficient being separately calculated for each of the 627 cells of the grid, taking into consideration only 30 neighbour cells. The spatial distribution of Pearson’s correlation coefficient and two-tailed probabilities for each class of values was obtained by means of the GWR. Subsequently, Moran’s I for the spatial autocorrelation [109] was performed for the residuals.

4. Results

4.1. Results of the Imagery Classification

By appliyng the supervised classification described in Section 3.1, the land use/land cover for 1989 (Figure 7a) and 2019 (Figure 7b) was derived.

Furthermore, the classification accuracy was assessd through the confusion matrix created with the help of data from the testing sample. Thus, based on the values from the confusion matrix for 1989 and 2019 (

Table 1. Confusion matrix for classification accuracy assessment.

Year.	Overall Acc. (%)	Kappa Index	Class	Ground Truth Samples (Pixels)							T.C. Pixels	User Acc. (%)
				B.A.	A.Z.	F.T.	P.	F.	T.W.	W.B.
1989	93.7	0.92	B.A.	347	0	2	0	2	0	7	358	96.9
			A.Z.	3	179	1	3	3	4	3	196	91.3
			F.T.	0	0	125	0	0	0	4	129	96.9
			P.	1	4	1	271	18	7	6	308	88.0
			F.	0	3	6	5	783	10	1	808	96.9
			T.W.	0	0	0	0	29	252	0	281	89.7
			W.B.	7	0	3	0	2	4	114	130	87.7
			T.G.T. Pixels	358	186	138	279	837	277	135	2210
			Prod. Acc. (%)	96.9	96.2	90.6	97.1	93.5	91.0	84.4
2019	95.23	0.939	B.A.	419	0	0	0	1	0	6	426	98.4
			A.Z.	2	245	7	0	0	0	0	254	96.5
			F.T.	0	2	56	0	0	0	0	58	96.6
			P.	10	0	0	173	9	0	6	198	87.4
			F.	0	5	1	7	706	25	8	752	93.9
			T.W.	2	0	0	7	0	229	0	238	96.2
			W.B.	2	0	0	0	0	0	172	174	98.9
			T.G.T. Pixels	435	252	64	187	716	254	192	2100
			Prod. Acc. (%)	96.3	97.2	87.5	92.5	98.6	90.2	90.2

B.A.—built-up areas; A.Z.—agriculture zones; F.T.—fruit trees; P.—pastures; F.—forests; T.W.—transitional wood-land; W.B.—water bodies; T.G.T. Pixels—total ground truth pixels; T.C. Pixels—total classified pixels; Prod. Acc.—producer accuracy.

), the overall classification accuracy and kappa index were calculated. For 1989, the overall classification accuracy was 93.7%, while the Kappa index was 0.92. Meanwhile, for 2019, the overall classification accuracy was 95.23%, while the Kappa index was 0.939. It should also be mentioned that for both analyzed years, the user accuracy and producer accuracy values were determined. Thus, for 1989, the user accuracy ranged between 87.7% (water bodies) and 96.9% (built-up areas), while the producer accuracy ranged between 84.4% (water bodies) and 97.1% (pastures). If we bring into discussion the year 2019, we observe that the user accuracy ranged between 87.4% (pastures) and 98.9% (water bodies), while the producer accuracy was situated between 87.5% (fruit trees) and 98.6% (forests).

4.2. Result of Land-Use/Land-Cover Changes

The Markov matrix (

Table 2. Markov matrix for the period 1989–2019.

1989\2019	Agricultural Areas	Transitional Woodland	Built-Up Areas	Forests	Pastures	Water Bodies	Fruit Trees	Losses (ha)
Agricultural areas	-	2.91	12.06	55.92	318.61	0	0	389.5
Transitional woodland	21.99	-	2.62	1554.83	97.87	0.27	12.43	1677.58
Built-up areas	564.13	0	-	34.39	3293.4	17.49	42.42	3909.41
Forests	53.73	4088.73	104.8	-	4088.73	23.3	0	8359.29
Pastures	906.6	666.66	220.01	1090.5	-	16.91	100.32	3000.9
Water bodies	7.45	16.58	51.37	48.75	32.91	-	0	157.06
Fruit trees	0	4.51	12.54	0	4.23	0	-	21.28
Gain (ha)	1553.9	4779.39	403.4	2784.39	7835.75	57.97	155.17	17,569.97

) was used in order to identify the direction and area of each land-use/land-cover transition. According to the matrix, the most dynamic classes were: agriculture areas, transitional woodlands, forests, and pastures. Agriculture areas faced an overall extension, especially by means of pasture conversion, but certain losses also occured as a consequence of the abandonment of several agricultural areas. The general balance revealed an increase of 1164.4 ha of agriculture zones during the studied period. The evolution of the transitional woodlands followed the inverse track of the forested areas, where forest areas decreased almost 5574.9 ha, which was 14.8% of the total forest area in 1989 (37,660 ha). According to the Markov matrix, the surface of built-up areas decreased approximately 3600 ha. This fact was due to the depopulation process that occurred after the political changes in 1989. Numerous types of transitions occurred, including conversions to agriculture zones, transitional woodlands, built-up areas, and pastures. The most important transition was from forests to pastures (4088.73 ha), mostly due to deforestation.

These transitions highlight major and complex land-use/land-cover changes that affected the Zăbala catchment between 1989 and 2019 involving all land-use/land-cover classes, a situation that strengthens our hypothesis that significat changes in runoff and flash flood potential occurred too.

In order to quantize the land-use/land-cover changes at the appropriate analysis scale, ${\mathrm{TRD}}_{\mathrm{SDLUI}}$ was calculated at the grid-cell level (1 km²) (Figure 10a). The spatial distribution of this indicator highlights the differences concerning the intensity of land-use/land-cover changes within the study area. The general pattern shows a west to east gradient, which suggests an increase in intensity of land-use/land-cover change that was simultaneous with the transition from the mountainous areas to the sub-Carpathian area due to the increase in accessibility offered by topographic conditions and by proximity to human communities.

The null values of this index, where no change occurred in land use/land cover, represented only 1% of the total area and composed an uneven spatial distribution, formed by several areas dispersed within the Carpathian zone. Low and intermediate values ranging from 0.1% to 16.85% had the greatest frequency (79% of the grid cells) (Figure 10a) and were distributed over the entire area: in the Carpathian zone, where shrubbery to forest changes mostly occurred due to the natural evolution of secondary vegetation installed after several old forests, and in the sub-Carpathian zone, especially along the valleys (Zăbala, Năruja, Petic), where different types of conversions took place: forest/pasture to built-up areas/agriculture areas. A total of 17% of the area recorded changes between 16.86% and 29.75%; this class formed less compacted areas, except several groups located at the limit between the Carpathians and sub-Carpathians, affected by deforestations and other types of conversion in the proximity of human settlements. The last class, where ${\mathrm{TRD}}_{\mathrm{SDLUI}}$ varied from 29.76% to 48.6%, characterized only 3% of the grid cells (Figure 10a) and was especially located in the sub-Carpathian zone where certain punctual and less extended changes occurred.

4.3. Run-off Risk Assessment (FFPI)

4.3.1. Flash-Flood Predictor Selection

It should be mentioned that the predictive ability was estimated for the years of the analysis. For 1989, the results of linear support vector machine (LSVM), in average merit (AM) values (

Table 3. Predictive ability of flash-flood predictors.

Flood Predictor	AM1989	AM2019
Slope	0.87	0.91
TPI	0.28	0.39
TWI	0.61	0.46
Land use/land cover	0.47	0.73
Lithology	0.52	0.59
Profile curvature	0.41	0.22
Aspect	0.17	0.13
Convergence index	0.35	0.32
Hydrological soil groups	0.23	0.26
MFI	0.68	0.55

), show that the slope had the highest predictive ability (AM = 0.87), followed by MFI₁₉₈₉ (AM = 0.68), TWI (AM = 0.61), lithology (AM = 0.52), land use/land cover₁₉₈₉ (AM = 0.47), profile curvature (AM = 0.41), convergence index (AM = 0.35), TPI (AM = 0.28), hydrological soil groups (AM = 0.23), and aspect (AM = 0.17). In terms of 2019, the highest predictive ability was also assigned to slope (AM = 0.91), followed by land use/land cover₂₀₁₉ (AM = 0.73), lithology (AM = 0.59), MFI₂₀₁₈ (AM = 0.55), TWI (AM = 0.46), TPI (AM = 0.39), convergence index (AM = 0.32), hydrological soil group (AM = 0.26), profile curvature (AM = 0.22), and aspect (AM = 0.13).

4.3.2. Application of the MLP Model for FFPI Computation

The first step in MLP training and in calculating FFPI₁₉₈₉ and FFPI₂₀₁₉ was done through the determination of FR coefficients. For 1989, the highest FR values were obtained by the following categories of land use/land cover: fruit trees (FR = 3.46), water bodies (FR = 2.45), shrub (FR = 2.21), and built-up areas (FR = 2.13). Hydrological soil group D (FR = 2.1) and slopes between 15 and 25° (FR = 1.98) also obtained high FR coefficients. For 2019, the highest FR values were achieved by the following class/categories of land use/land cover and lithology: water bodies (FR = 3.7), shrubs (FR = 3.16), sandstone and conglomerates (FR = 3), sandstone and tuffs (FR = 2.79), sandstone and marls (FR = 2.79), and slopes between 25 and 45° (FR = 2.72) (

Table 4. Frequency ratio and multilayer perceptron weights.

Factor	Class	FR1989	FR1989 N	FR2019	FR2019 N	MLP1989 Weight	MLP2019 Weight
Slope	0–3°	0.00	0.10	0.33	0.20	0.373	0.404
	3–7°	0.20	0.18	0.00	0.10
	7–15°	0.43	0.27	0.25	0.17
	15–25°	1.98	0.90	1.93	0.67
	25–45°	1.20	0.58	2.72	0.90
TPI	−25.1 to −4.74	1.47	0.90	1.35	0.90	0.097	0.103
	−4.73 to −1.31	0.93	0.33	1.21	0.75
	−1.3 to 1.52	1.01	0.41	1.03	0.56
	1.53–5.15	0.98	0.38	0.70	0.21
	5.16–26.32	0.72	0.10	0.60	0.10
TWI	3.09–6.04	1.27	0.90	1.14	0.66	0.104	0.194
	6.05–7.74	1.13	0.81	1.09	0.64
	7.75–10.1	0.14	0.10	0.37	0.28
	10.11–14.75	0.90	0.64	1.62	0.90
	14.76–24.64	0.67	0.48	0.00	0.10
Land use/land cover	Built-up areas	2.13	0.59	2.21	0.56	0.287	0.236
	Agriculture zone	0.28	0.15	1.01	0.29
	Shrub	2.21	0.60	3.16	0.78
	Fruit trees	3.46	0.90	2.58	0.65
	Forests	0.08	0.10	0.16	0.10
	Pastures	0.64	0.23	0.37	0.15
	Water bodies	2.45	0.66	3.70	0.90
Lithology	Sandstone flysch	0.91	0.65	0.58	0.17	0.146	0.128
	Gravels, sands	0.00	0.10	1.02	0.30
	Clay with blocks	0.90	0.64	0.79	0.23
	Sandstone, shales	1.12	0.77	1.54	0.46
	Sandstone, marls	0.47	0.38	2.79	0.84
	Sandstone, tuffs	0.70	0.52	2.79	0.84
	Sandstone, conglomerates	1.09	0.75	3.00	0.90
	Sandstone-shale	1.34	0.90	0.35	0.10
Profile curvature	0.9–1.4	1.52	0.90	1.79	0.90	0.059	0.045
	0–0.9	0.86	0.10	1.11	0.40
	−1.6 to 0	0.99	0.26	0.71	0.10
Aspect	Flat surfaces	0.00	0.10	0.00	0.10	0.051	0.032
	North	0.98	0.63	1.09	0.86
	North-East	1.32	0.81	1.12	0.88
	East	1.25	0.77	0.92	0.74
	South-East	0.72	0.49	0.72	0.60
	South	0.79	0.53	1.15	0.90
	South-West	0.67	0.46	1.00	0.80
	West	0.56	0.40	0.96	0.77
	North-East	1.48	0.90	1.06	0.84
Convergence index	−99.3 to −3	0.76	0.10	0.40	0.14	0.05	0.084
	−3 to −2	0.79	0.16	0.34	0.10
	−2 to −1	1.13	0.90	1.04	0.54
	−1–0	1.10	0.85	1.61	0.90
	0–100	1.09	0.82	1.21	0.65
(Hydrological Soil Group (HSG)	A	1.70	0.67	1.70	0.10	0.043	0.067
	B	0.70	0.10	1.90	0.29
	C	1.51	0.57	2.10	0.47
	D	2.10	0.90	2.56	0.90
MFI	<60	1.01	0.10	1.34	0.10	0.122	0.289
	60–90	1.24	0.43	1.52	0.52
	90–120	1.56	0.90	1.68	0.90
	>120	1.32	0.55	1.48	0.43

).

The RMSE for 1989 was equal to 0.0183, while for 2019, the RMSE was equal to 0.0157. Both values denote a very good performance for the MLP algorithm.

Furthermore, through MLP, the importance of each flash-flood predictor was estimated for each of the two years. Thus, for 1989, the most important factor was slope angle (0.373), followed by land use/land cover₁₉₈₉ (0.287), lithology (0.146), MFI₁₉₈₉ (0.122), TWI (0.104), TPI (0.097), profile curvature (0.059), aspect (0.051), convergence index (0.05), and hydrological soil groups (0.043). In can be noted that for 2019, the slope angle also had the highest importance (0.404), followed by MFI₂₀₁₉ (0.289), land use/land cover (0.236), TWI (0.194), lithology (0.128), TPI (0.103), convergence index (0.084), hydrological soil group (0.067), profile curvature (0.045), and aspect (0.032).

Through the multiplication of MLP weights with the FR normalized coefficient, the FFPI values were obtained. The highest FFPI values were concentrated in the eastern part of the Zăbala catchment for both 1989 and 2019 (Figure 8a and Figure 8b, respectively). Several areas with a high FFPI occurred in the western part of the study area, on high-altitude pastures, having a declivity higher than 15° and developed on soils with a clay texture. For 1989, the high and very high FFPI values covered approximately 34% of the total area (Figure 8a), while in 2019, they represented 46% of the Zăbala catchment (Figure 8b). Intermediate values for FFPI for both years were evenly distributed across 29–30% of the total area. The low and very low values of FFPI₁₉₈₉, ranging from 0.1 to 0.31, were associated with 35% of the catchment, while the same values of FFPI₂₀₁₉ were distributed over 24% of the Zăbala river basin.

4.3.3. FFPI Results Validation

As was presented in Section 3.3.5, the results validation and model’s performance evaluation were done through the ROC curve and AUC. It can be observed that in terms of the success rate, both the MLP models obtained very good results. Thus, the AUC for MLP₁₉₈₉ was equal to 0.895, while for MLP₂₀₁₉, the AUC was equal to 0.874 (Figure 9a). Very good results were also achieved after the application of the prediction rate. Thus, the AUC for MLP₁₉₈₉ was 0.864, while for MLP₂₀₁₉, it was 0.837 (Figure 9b). These results highlight a high degree of accuracy associated with the FFPI maps. Therefore, the results obtained were used for further analysis.

4.3.4. FFPI Differences

Regarding the RD_FFPI spatial variations, the largest surfaces, which accounted for approximately 59% of the study area, was covered by the changes that occurred in a percentage between 0.01% and 3.13%. The class of changes between 3.14% and 5.2% was found for 27% of the study area. The null changes and the changes of between 8.3% and 13.96% occurred on approximately 2% of the Zăbala river catchment (Figure 10b).

4.4. Statistical Analysis for Correlating $TR{D}_{SDLUI}$ and $R{D}_{FFPI}$

By analyzing Figure 11, apart from the obvious overall correlation between the two variables, several particular spatial patterns can be observed. In the northern part of the study area (Năruja valley), serious land-use/land-cover changes occurred without clearly influencing the flash-flood potential. This situation can be explained by the conversion between land uses/land covers having almost similar impacts on flash-flood potential (forest, fruit trees, or transitional wood-lands). Several areas faced intermediate values for the correlation coefficient (0.33–0.48) and were mainly distributed at the contact between the sub-Carpathian and Carpathian zones. The dominant class had correlation coefficients ranging from 0.67 to 0.94 and was distributed in areas where both phenomena (land-use/land-cover change and flash-flood potential) increased during the studied period, such as: the periphery of rural settlements (Paltin, Spulber), high mountain peaks, and several valley sectors (Zăbala, Lapos). The correlation was finally validated using a spatial autocorrelation test, i.e., Moran’s Index, which highlighted a random distribution for the residuals of GWR (Figure 11).

5. Discussion

The image classification procedure using the remote sensing techniques had a crucial role in the present paper because it allowed us to assess the changes that occurred in land use/land cover during the last 30 years, and further, to analyze the correlation between these changes and those that occurred in flash-flood potential. The importance of this remote sensing technique it is even greater as the results obtained had a very high accuracy. According to the results achieved by Pal and Mather [110], the accuracy obtained by the maximum likelihood method (82.9%) was close to the accuracy of the artificial neural network (85.1%) and higher than the accuracy of the support vector machine classification (79.73%). The fact that the maximum likelihood achieved accuracies almost equal to or higher than more advanced techniques in some situations is a solid argument that this method does not have relevant limitations against the quality of the outputs. The reliability of this method was also highlighted by the accuracy of the classification undertaken in the present study. Thus, according to the results presented above, the maximum likelihood supervised classification method performed very well on the Zăbala river catchment and provided very precise results with an overall accuracy higher than 93% and the Kappa index higher than 0.92. The accuracy of these results was higher than the accuracy of 81.4% achieved by Asamoah et al. [111], who used the maximum likelihood algorithm to classify a Landsat scene in a study carried out on a region of Ghana. An accuracy of 91.34%, closer to that obtained in the present study, resulted from a study carried out by Ali et al. [112], who classified a Landsat scene from Egypt through the maximum likelihood method. Ajaj et al. [113] managed to extract the land use/land cover from Landsat 7 Enhanced Thematic Mapper Plus (ETM+) imagery with an overall accuracy of 94.25% and a value of Kappa index of 0.921. In the present research, the results achieved through remote sensing techniques were analyzed with the help of the Markov matrix, which revealed that the forest surfaces decreased by more than 5500 hectares. Due to the fact that this area represented approximately 13.1% from the total surface covered by forest in 1989, we consider that deforestation was an important factor that contributed to the increase of flash-flood potential across the Zăbala catchment. The situation of these areas is all the more worrying as the deforestations occurred in areas with high slopes and a high degree of hydrographic convergence that favours, due the absence of the forest, a very high potential for flash-flood genesis.

In this regard, there is no doubt of the fact that the assessment of flash-flood potential, especially within the mountain and hilly catchments, is of a crucial importance for future measurements that the authorities should take in order to mitigate the flash-flood damages. This is also the case for the present study, in which the FFPI was estimated through the multilayer perceptron method. It should be remarked that the use of the MLP model provided accurate results with AUC values higher than 0.8. The MLP model was also involved in the determination of the FFPI values across the Prahova river catchment from Romania [87], where this machine learning algorithm also provided very accurate results (AUC > 0.8). An optimization of multilayer neural networks was tried by Ngo et al. [49], who attempted to estimate the flash-flood susceptibility across a mountain catchment from Vietnam. In this study, the AUC values of the optimized neural networks were also situated above 0.8, which indicates the high accuracy of machine learning model [114,115,116,117,118,119,120,121].

Nevertheless, the present research did not stop with a simple assessment of the flash-flood potential and was further developed with the analysis of the correlation between the changes in the land use/land cover and flash-flood potential between 1989 and 2019. To our best knowledge, the results of the present research represent a first for the international literature because a study regarding the correlation between the land-use/land-cover changes and the bi-temporal evolution of flash-flood potential index (FFPI) has not been made before. Our research highlights that the land-use/land-cover changes between 1989 and 2019 determined the increase in the flash-flood potential across the Zăbala catchment. Numerous studies in the literature have analyzed the quantitative influence of land-use/land-cover changes in the runoff characteristics using the SWAT model [122,123,124,125,126]. The vast majority of these studies, in which the influence of land-use/land-cover changes on the runoff characteristics is explored, indicate that during the last few decades, the aggressiveness of run-off and flash-flood frequency significantly increased, mainly due to deforestation activities [127,128,129].

6. Conclusions

The present study, concerning the connection between land-use/land-cover changes between 1989 and 2019 and flash-flood potential evolution in the Zăbala catchment holds great practical importance in the context of the multiplication of flash floods events that represent hydrological risk phenomena, especially affecting the sub-Carpathian zone of the basin.

The applied methodology was based on the correlation between the values of TRDSDLUI, which highlighted the land-use/land-cover changes between 1989 and 2019, and RD_FFPI values, which highlighted the changes within the flash-flood potential index. The values of TRDSDLUI were derived on 1 km² grid cell size with the help of remote sensing and GIS techniques, while for computing the RD_FFPI values on the same grid cell size, the use of machine learning (MLP) had a crucial role. The high efficiency of the MLP algorithm, in the case of FFPI mapping, was demonstrated by the AUC values of the ROC curve, which were between 0.837 and 0.895.

Having obtained the highly accurate results for both indicators of land-use/land-cover changes and flash-flood potential changes, the correlation between these two processes were assessed through the geographically weighted regression. The analysis tools provided by geostatistics were efficient as they were able: (1) to confirm the overall correlation between the two phenomena; and (2) to identify the local particularities of each analyzed spatial unit, which were highlighted by the spatial variability of the Pearson correlation coefficient.

Given the results presented above, we can state that the changes that occurred in the land use/land cover were correlated to a high degree with the changes that occurred in the FFPI, especially in the eastern, southern, and southwestern part of the Zăbala river catchment. Therefore, it can be assumed that, within the study area, the land-use/land-cover changes had a crucial role in the intensification of flash-flood events, which were more and more frequent within the study area.

The main novelty of the present study is the fact that, to our best knowledge, it is the first study that addresses the problem of land-use/land-cover changes in relation to the changes within the flash-flood potential index at international level. Also, taking into account that, nowadays, the subject of the effects of land-use/land-cover changes regarding the intensification of floods and flash-floods is a very highly discussed issue in the literature, the present study can be a very solid benchmark for future studies and new research directions.

The results of the present study can also be used in order to efficiently plan the use of territory from the flood risk management point of view.