Assessing Hydrological Impact of Forested Area Change: A Remote Sensing Case Study

The interaction between precipitation and vegetation plays a significant role in the formation of runoff, and it is still a widely discussed issue in hydrology. The difficulty lies in estimating the degree to which a forest influences runoff generation, especially flood peaks, on the one hand, due to the small amount of information regarding the evolution of the forest area and density, and, on the other hand, the correlations between these indicators and the runoff and precipitation values. The analysis focuses on a small basin in the mountain region of Romania, the upper basin of the Ruscova River located in northwestern Romania. In this river basin, there is no significant anthropic influence, other than the intense deforestation and reforestation actions. Using satellite images captured by Landsat missions 5, 7 and 8 for the period 1985–2019, the forest canopy density vegetation index was extracted. Using a gridded precipitation dataset, a hydrological model was calibrated based on three scenarios to assess the impact of forest vegetation on the runoff. Analysis of the results of these models conducted on scenarios allowed us to deduce a simple equation for estimating the influence of forest area on maximum river flows. The analysis showed that even small differences in the forest surface area exert an influence on the peak flow, varying between −5.28% and 8.09%.


Introduction
Land use and land cover change (LULCC) represents one of the most important factors that alters the natural runoff regime. Andréassian (2004) [1], citing Pliny the Elder, from the first century AD, already identified a correlation between forest removal and the subsequent increase in streamflow. At the global level, many researchers have tried to estimate/evaluate the impacts of LULCC on the evolution and changes in runoff regimes [1][2][3][4].
Forests play an important role in the natural water balance, influencing runoff and generally regions' water resources. Anthropogenic deforestation and reforestation are among the most important land use and land cover changes worldwide [5]. In forested catchments, the climate change impact is significant since a rising temperature can reduce trees' vitality even though precipitation amounts do not change [6]. A recent study found that, in some cases, changes in land use influence the surface runoff response more significantly than climate change [7].
According the PannEx White Book, it is particularly vital to study the effect of deforestation on catchment runoff [8]. Understanding the influence of LULCC on river In the study area, a snow-dominated hydrologic regime is characteristic. In spring, there are high waters dominated by melting snow and high amounts of precipitations; in summer and autumn, waters are generally low due to the global evaporation which is interrupted by floods due to torrential rains. At the end of autumn and in winter, there are also low waters: in autumn due to the lack of precipitation, and in the cold season because a part of the water is retained in the ice formations and because the precipitations are in a solid form.
Concerning the altitude of the basin, the upper basin of the Ruscova River is located between 584 and 1944 m. The topography is characterized by steep slopes (the average slope of the basin is 42.52%).
From the total surface of the upper basin of the Ruscova River, 69% (in the year 2018) is covered by forests ( Figure 2). The predominant types of forests are mixed forest, formed by coniferous and broad-leaved forests (26.5%), followed closely by coniferous (25.2%) forests and, with a smaller weight, broad-leaved forests (17.3%). Among the tree species, fir and spruce are dominant, but soft-deciduous mixed forests are present too. In the study area, a snow-dominated hydrologic regime is characteristic. In spring, there are high waters dominated by melting snow and high amounts of precipitations; in summer and autumn, waters are generally low due to the global evaporation which is interrupted by floods due to torrential rains. At the end of autumn and in winter, there are also low waters: in autumn due to the lack of precipitation, and in the cold season because a part of the water is retained in the ice formations and because the precipitations are in a solid form.
Concerning the altitude of the basin, the upper basin of the Ruscova River is located between 584 and 1944 m. The topography is characterized by steep slopes (the average slope of the basin is 42.52%).
From the total surface of the upper basin of the Ruscova River, 69% (in the year 2018) is covered by forests ( Figure 2). The predominant types of forests are mixed forest, formed by coniferous and broad-leaved forests (26.5%), followed closely by coniferous (25.2%) forests and, with a smaller weight, broad-leaved forests (17.3%). Among the tree species, fir and spruce are dominant, but soft-deciduous mixed forests are present too.
The basic criteria for choosing this river basin were the following: • The surface of the river basin must be covered as much as possible by forest; • There are no human settlements or elements of infrastructure that would lead to the permanent reduction in forested areas; • Afforestation and deforestation actions were detected in the focus area over recent decades.
The upper basin of the Ruscova River is characterized by large areas covered by forests. Over time (especially in the period 1995-2010, in the focus region, more or less legal deforestation has been carried out, or numerous destructions due to natural causes (storms or erosions) have occurred. There were also several campaigns to reforest the area by ROMSILVA (Romania's National Forest Administration is the Romanian state company in charge of protecting and managing the Romanian public property forests). The basic criteria for choosing this river basin were the following: • The surface of the river basin must be covered as much as possible by forest; • There are no human settlements or elements of infrastructure that would lead to the permanent reduction in forested areas; • Afforestation and deforestation actions were detected in the focus area over recent decades.
The upper basin of the Ruscova River is characterized by large areas covered by forests. Over time (especially in the period 1995-2010, in the focus region, more or less legal deforestation has been carried out, or numerous destructions due to natural causes (storms or erosions) have occurred. There were also several campaigns to reforest the area by ROMSILVA (Romania's National Forest Administration is the Romanian state company in charge of protecting and managing the Romanian public property forests).
Out of the total area of the focus region, 70.1% is covered by acid brown soils. These types of soil are formed in cold and humid climates, which causes the organic matter to be only partially decomposed by microorganisms, and the very slow process leads to an increased acidification of the soil. They are poor in nutrients with a high-permeability texture and are generally covered with forest vegetation. Along with acid brown soils, lithosols (1.3%) and feriiluvial brown soils (podzolic) (28.6%) are present. Out of the total area of the focus region, 70.1% is covered by acid brown soils. These types of soil are formed in cold and humid climates, which causes the organic matter to be only partially decomposed by microorganisms, and the very slow process leads to an increased acidification of the soil. They are poor in nutrients with a highpermeability texture and are generally covered with forest vegetation. Along with acid brown soils, lithosols (1.3%) and feriiluvial brown soils (podzolic) (28.6%) are present.

Precipitation
Precipitation data were extracted from the daily gridded dataset provided by the European Climate Assessment and Dataset project (EOBS) [25,26]. The EOBS daily gridded database is an ensemble dataset developed at the European level for the period 1950-2020 (updated year by year). The spatial resolution of the grid is 0.
The comparative analysis between the EOBS gridded dataset and the observed values for the NW area of Romania recommended their use for this type of analysis [27].

Runoff Data
For the runoff data, we used one gauging station which is part of the Romanian national network of hydrology and water management: Luhei Hydrometric Station on the Ruscova River (Figure 1), and the data were daily mean discharge values for the entire study period.

Satellite Images
Analysis of the evolution of the forest area was made based on the Landsat Multispectral Scanner (MSS) Surface Reflectance images freely provided by the United States Geological Survey (USGS) for the period 1985-2019 (Landsat 4-8) [28]. All used images are cloud-free.
The analysis of the method's accuracy was performed based on the orthophoto images provided by the Romanian National Cadastral Agency and Real Estate Advertising service [29].

Methods
The forest plays a significant role in the soil-water-atmosphere system [30]. Thus, precipitation water partially seeps into the looser soil of the forest, another part is retained by the canopy of trees and subjected to the process of evapotranspiration and perspiration, another amount is consumed by tree roots and transformed into sap and another significant part flows as surface runoff [1,10,30]. Therefore, in the analysis of the influence exerted by the forest on the runoff, the following information and data are needed:

•
The amount of precipitation falling in a certain area; • The evolution of the forest area over a certain period of time (afforestation and deforestation actions) in the focus region; • The evolution of potential evapotranspiration; • The amount of water drained (flow) at the closing section of the analyzed basin.

Estimating the Amount of Precipitation
Correct determination of the amount of precipitation is an essential step in this analysis. Using the geographical limits of the basin, all the daily data series contained in the EOBS gridded dataset for the period 1985-2019 were extracted. Then, the extracted values for each grid were mediated to determine the average amount of precipitation falling in the watershed.

Evolution of Forest Area and Density
In the analyzed area, for the period after 1985, no first-hand reliable information was found regarding the evolution of forested areas. To follow the evolution of the afforested surfaces in the study area, several vegetation indices were extracted and analyzed from the Landsat images for the period 1985-2019. The following issues were taken into account:

1.
The image collection period should be in the interval between the end of May and the start of October so that the phenophase of the forest is somewhat similar, fully developed canopy; 2.
The analyzed surfaces should not be covered by clouds.
Since the research area is located in a mountain region characterized by high cloudiness, it was quite difficult to identify images without significant cloud cover for each year after 1985. Table 1 shows the data and characteristics of satellite images considered in this analysis; all images have a 30 m spatial resolution and are from the T1 collection category.
Landsat Surface Reflectance images provided by the U.S. Geological Survey were used [28]. These images measure the fraction of incoming solar radiation reflected from Earth's surface to the Landsat sensor, for Landsat missions 5, 7 and 8.
All images are converted within an 8-bit range.
To choose the best method, we tested several vegetation indices (VI) obtained based on the analysis of the spectral bands of Landsat satellite images: leaf area index (LAI), normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), a soiladjusted vegetation index (SAVI), forest canopy density (FCD). In this evaluation, we used predefined indices in different programs, such as the Raster Function module of ArcGIS PRO software (produced by ESRI US) for NDVI and LAI indices, or we analyzed the satellite images using the software function Sentinel Application Platform (SNAP) produced by the European Space Agency (ESA) for SAVI and EVI indices. The best results were obtained using the method proposed in [31][32][33] for the determination of the FCD index. It is defined as the proportion of coverage of the forest floor by the vertical projection of the peaks of trees, and it is also known as crown coverage or canopy coverage [34,35]. The FCD index is achieved by composing four other indices: advanced vegetation index (AVI), canopy shadow index (SI), bare soil index (BI), thermal index (TI).
Using the Model Builder implemented in the ArcGIS desktop program, a workflow was implemented (see Supplementary Material Figure S1) to obtain the FCD map, as described in Figure 3.  Image Normalization. The first step in obtaining the four already mentioned indicators was the normalization of the spectral bands. This stage involves the elimination of distortions caused by the angle of the Sun with the vertical of the place at the zenith, the average Earth-Sun distance when taking the image, atmospheric conditions, etc. [36].
The Landsat TM bands (except for band 6) were normalized using linear transformation (Equations (1)-(3)): where 1 = − 2 ; 2 = + 2 ; 1 = 20; 2 = 220; M-mean; S-standard deviation; X-original data; and Y-normalization data. Image Normalization. The first step in obtaining the four already mentioned indicators was the normalization of the spectral bands. This stage involves the elimination of distortions caused by the angle of the Sun with the vertical of the place at the zenith, the average Earth-Sun distance when taking the image, atmospheric conditions, etc. [36].

Obtaining the Four Vegetation Indices.
According to the methodology proposed by Rikimaru [31,32], we employed the following vegetation indices, derived from satellite images.
Advanced Vegetation Index (AVI). This method uses the combination of red and near-infrared (NIR) bands. By raising the power, the spectral signature of chlorophyll is amplified. This way, the index is able to capture the factors and type of vegetation (it is able to identify vegetation inversions, but especially the anthropogenic influence of reforestation). Equation (4) shows how to obtain this indicator: where: B4-near-infrared band; B3-red band.

Bare Soil Index (BI).
Adding an index to identify bare soil surfaces, in combination with vegetation indices, results in a stronger contrast in the analysis [31][32][33]. The basic logic of this approach is to build on the high reciprocity between the state of the bare soil and the state of vegetation. By combining vegetation and bare soil indices in the analysis, one can assess the state of forests on a continuum ranging from high vegetation conditions to exposed soil conditions. BI ( Figure 5) was calculated by using Equation (5): where: B5-shortwave infrared band; B1-blue band.

Bare Soil Index (BI).
Adding an index to identify bare soil surfaces, in combination with vegetation indices, results in a stronger contrast in the analysis [31][32][33]. The basic logic of this approach is to build on the high reciprocity between the state of the bare soil and the state of vegetation. By combining vegetation and bare soil indices in the analysis, one can assess the state of forests on a continuum ranging from high vegetation conditions to exposed soil conditions. BI ( Figure 5) was calculated by using Equation (5): Atmosphere 2021, 12, 817 8 of 24 to exposed soil conditions. BI ( Figure 5) was calculated by using Equation (5): where: B5-shortwave infrared band; B1-blue band.
(a) (b) Shadow Index (SI). The use of this index highlights the three-dimensional characteristics of the forest by placing the focus on the shaded areas of the forest canopy. It is based on the low brightness of the visible bands. Through this process, younger tree structures where: B5-shortwave infrared band; B1-blue band. Shadow Index (SI). The use of this index highlights the three-dimensional characteristics of the forest by placing the focus on the shaded areas of the forest canopy. It is based on the low brightness of the visible bands. Through this process, younger tree structures can be identified [31][32][33]. These structures usually have a lower brightness than more mature trees. At the same time, this index allows the delineation of the surfaces covered by bushes and shrubs from the surfaces with trees.
The shadow index is calculated using Equation (6): where: B2-green band. The SI value is directly proportional to the height of the canopy and the density of the forest ( Figure 6). can be identified [31][32][33]. These structures usually have a lower brightness than more mature trees. At the same time, this index allows the delineation of the surfaces covered by bushes and shrubs from the surfaces with trees. The shadow index is calculated using Equation (6): where: B2-green band. The SI value is directly proportional to the height of the canopy and the density of the forest ( Figure 6). Thermal Index (TI). This index was developed considering that bare surfaces (without vegetation) have a higher temperature at ground level than shaded surfaces in areas covered by forests [33]. This indicator allows distinguishing between the shaded slopes with forests and the shaded slopes without forests.
The first step in obtaining the brightness temperature is to convert the reflectance values into radiant values by using relation 7: Thermal Index (TI). This index was developed considering that bare surfaces (without vegetation) have a higher temperature at ground level than shaded surfaces in areas covered by forests [33]. This indicator allows distinguishing between the shaded slopes with forests and the shaded slopes without forests.
The first step in obtaining the brightness temperature is to convert the reflectance values into radiant values by using relation 7: where: L-spectral radiance at the pixel level (meter squared * ster* µm): DN-digital number or spectral value of the pixel. By employing the spectral radiance and some calibration coefficients (K1, K2), the temperature at ground level is derived, as presented in Equation (8): For Landsat missions 5 and 7, K1 = 666.09 watts (m2/sr/µm), and K2 = 1282.71 K; for Landsat 8, K1 = 774.8853 watts (m2/sr/µm), and K2 = 1321.0789 K.
In Figure 7, images with the TI for the years 1989 (a) and 2017 (b) are presented. Combining Indices to Achieve the FCD Model.
The composition of the FCD model implies the derivation from the initially obtained indices of two other indicators: the vegetation density (VD) and the scaled shadow index (SSI).
VD is obtained by synthesizing the AVI and BI indices by the principal component analysis method [31][32][33]. Next, each pixel value is divided by the maximum value and multiplied by 100 to obtain the percentage value.
The SI index is an index whose value depends on the low irradiant area of each visible band. The use of the SI indicator to determine the SSI may present false information due to the possibility of delimiting including the shaded slopes. As bare soils have a higher temperature than the shady ones (due to their higher capacity to absorb solar radiation), the SI was combined with the TI index. Thus, the possible errors generated by the exposure of the slopes are avoided [31][32][33]. In the next step, this indicator was transformed by relating each pixel to the absolute maximum value. The SSI value is scaled from 0 to 100%: the value of SSI 100 corresponds to a forest with the highest shadow value, whereas the value of 0 indicates a forest with the lowest shadow value. Thus, in areas where the SSI is VD is obtained by synthesizing the AVI and BI indices by the principal component analysis method [31][32][33]. Next, each pixel value is divided by the maximum value and multiplied by 100 to obtain the percentage value.
The SI index is an index whose value depends on the low irradiant area of each visible band. The use of the SI indicator to determine the SSI may present false information due to the possibility of delimiting including the shaded slopes. As bare soils have a higher temperature than the shady ones (due to their higher capacity to absorb solar radiation), the SI was combined with the TI index. Thus, the possible errors generated by the exposure of the slopes are avoided [31][32][33]. In the next step, this indicator was transformed by relating each pixel to the absolute maximum value. The SSI value is scaled from 0 to 100%: the value of SSI 100 corresponds to a forest with the highest shadow value, whereas the value of 0 indicates a forest with the lowest shadow value. Thus, in areas where the SSI is 100, forest is present.
FCD. This step involves combining the VD and SSI (Equation (9)) [31][32][33]. As the values of the two parameters represent the percentage of values, the FDC ranges from 0 to 100:

Evaluation of the Method for Determining the FCD Index
For evaluating the method, we used the kappa coefficient (κ), which is employed to measure the accuracy of a method by evaluating the actual values, derived from observations and the estimated ones, determined by calculation [37].
To perform this analysis, we created 20 random points that we intersected with the FCD layers made for each analyzed year. For intersections with forest areas, they were quantified as 1; for intersections with pasture areas, we used 0. The comparison of these estimated values with the observed values was conducted by comparing them with aerial photographs from 2018, where the same scores of 1 and 0 were given.

Analysis of Precipitation Values and Flow Hydrographs Related to the Evolution of Forest Surfaces
The entire amount of precipitation that falls in a region does not reach the hydrographic network. Part of this amount is returned to the atmosphere by evapotranspiration, another part is absorbed by vegetation at the interception level and yet another part infiltrates and forms underground runoff and/or hypodermic runoff.
The comparative analysis of the discharge values with the surface covered by forest or a comparative analysis between runoff coefficients with the forest surface values can give us an image of the influence of the forest surfaces on the river discharge in the analyzed area. However, this comparative analysis is not sufficient to determine the impact that the forest has on the runoff. For such an analysis, a hydrological model able to simulate hydrological processes is required. One of the setting parameters of the model should be the vegetation, or an index derived from it.
For this study, we used the well-known unit hydrograph method (UHM) from the Mike Hydro River modeling system produced by DHI [38,39]. This type of model cannot simulate events in which snow influences the flow of the watercourse. That is why this model is not recommended to be used for those events in which snow or melting snow influences the runoff. In the UHM module, the excess rain is calculated assuming that the losses to infiltration can be described by the SCS general-

Evaluation of the Method for Determining the FCD Index
For evaluating the method, we used the kappa coefficient (κ), which is employed to measure the accuracy of a method by evaluating the actual values, derived from observations and the estimated ones, determined by calculation [37].
To perform this analysis, we created 20 random points that we intersected with the FCD layers made for each analyzed year. For intersections with forest areas, they were quantified as 1; for intersections with pasture areas, we used 0. The comparison of these estimated values with the observed values was conducted by comparing them with aerial photographs from 2018, where the same scores of 1 and 0 were given.

Analysis of Precipitation Values and Flow Hydrographs Related to the Evolution of Forest Surfaces
The entire amount of precipitation that falls in a region does not reach the hydrographic network. Part of this amount is returned to the atmosphere by evapotranspiration, another part is absorbed by vegetation at the interception level and yet another part infiltrates and forms underground runoff and/or hypodermic runoff.
The comparative analysis of the discharge values with the surface covered by forest or a comparative analysis between runoff coefficients with the forest surface values can give us an image of the influence of the forest surfaces on the river discharge in the analyzed area. However, this comparative analysis is not sufficient to determine the impact that the forest has on the runoff. For such an analysis, a hydrological model able to simulate hydrological processes is required. One of the setting parameters of the model should be the vegetation, or an index derived from it.
For this study, we used the well-known unit hydrograph method (UHM) from the Mike Hydro River modeling system produced by DHI [38,39]. This type of model cannot simulate events in which snow influences the flow of the watercourse. That is why this model is not recommended to be used for those events in which snow or melting snow influences the runoff. In the UHM module, the excess rain is calculated assuming that the losses to infiltration can be described by the SCS generalized CN (Soil Conservation Service curve number) method [39], which is a simple and efficient method for determining the approximate amount of runoff from a rainfall event in a particular area. The curve number is a dimensionless index, which can take values from 0 to 100. It is based on both land use and the hydrologic soil group and reflects the potential to generate water runoff on different surfaces. By processing the DEM using the spatial analysis tools of the ArcGIS Pro software (ESRI, US), a series of parameters required by the model was extracted: slope angle ( Figure  9b), flow direction, flow length, sub-basin delineation.
The data source referring to the soil texture, used in the model to determine the value of the CN index (Figure 9c), is the National Research-Development Institute for Pedology, Agrochemistry and Environmental Protection Bucharest (ICPA). The spatial resolution of the soil map is 1:25,000. By processing the DEM using the spatial analysis tools of the ArcGIS Pro software (ESRI, US), a series of parameters required by the model was extracted: slope angle (Figure 9b), flow direction, flow length, sub-basin delineation.
The data source referring to the soil texture, used in the model to determine the value of the CN index (Figure 9c), is the National Research-Development Institute for Pedology, Agrochemistry and Environmental Protection Bucharest (ICPA). The spatial resolution of the soil map is 1:25,000.
Based on the river basin catchment area, we extracted the related precipitation data series. The spatial representation of the extracted values is points located in the center of each grid cell. Further on, based on the grid pixel centers, Thiessen polygons were derived, and their overlap in the basin area allowed us to calculate a weighted average of precipitation values for the considered area (Figure 9d).

Soil Conservation Service (SCS) Generalized (CN) Loss Method
As mentioned above, the CN is based on the area's hydrologic soil group, the type of land use and the hydrologic condition. The U.S. Soil Conservation Service developed this method for computing losses from storm rainfall in 1972 [37][38][39].
The general equation for the SCS CN method is as in Equation (10): where: Q-runoff(m 3 /s); P-precipitation (mm); S-potential maximum retention after runoff begins; I a -user-defined initial loss (mm). According to the DHI Mike 11 Reference Manual [38,39], an empirical relation was developed after testing several small experimental watersheds (Equation (11)).
The potential maximum retention S is calculated from a dimensionless CN using the empirical formula derived by SCS based on rainfall runoff analyses of a large number of catchments. The value of S can be calculated with the formula (Equation (12)) In this analysis, we used the SCS generalized method. The main difference from the SCS CN method is that the SCS generalized method does not make use of the concept of an antecedent moisture content (AMC) but applies a storage initial abstraction depth [38][39][40].

Establishing the Soil Hydrological Groups
As presented in Equation (11), the maximum retention potential is based on the CN, which is determined based on the hydrological group of soils and land use at the level of the river basin. From this point of view, the soils were classified into 4 hydrological groups (A, B, C and D), depending on the water infiltration rate. Chendes [41] adapted the soil hydrological groups classification to the Romanian soil classification system. Here, we used the same approach. The classification is based on the texture of the soil, assessing the percentage of sand, silt and clay in each type of soil, which provides a characteristic from permeable to impermeable soil or from light (sandy) to heavy soil (clays) (see Supplementary Material Table S1).
The analyzed area does not have a great diversity of soil textures. Only three types of textures were identified: sandy-loam-loamy (75.7%), sandy-loam (23.8%) and sandy-loamclay-loam (0.5%). The high percentage of the sandy-loamy soil texture is due to the fact that most of the focus region is forest land, and the main species are coniferous (Figure 9c). In order to determine a single value of the CN index, a weighted average value was computed using the corresponding area for every land use type (Equation (13)): where: (x)-area weighted average CN; c i -certain land use area; x i -the values x 1 , x 2 . . . x n represent the CN values from the original reference table [39,41]; A-total area of the river basin. The above procedure is necessary to measure the contribution of the land use type (in our case, forested areas) to the runoff. In essence, by modifying the surface area of the forest and the pastures, we automatically obtain a different CN, which influences the simulation of the runoff during a certain event.

Model Calibration and Event Isolation
UHM is used to estimate the runoff of a single storm event using the well-known unit hydrograph technique.
The models divide storm rainfall into excess rainfall (or runoff) and water loss (or infiltration). From a conceptual point of view, this model represents the average response of a unitary surface of the hydrographic basin to the impulse generated by the excess of precipitations. Therefore, for each calculation time step, the model determines the response of the basin to excess rainfall, as generated by the loss model during the time step, and cumulates this response with the flow values formed in the previous time steps [38].
The following parameters were used for model calibration: lag time, hydraulic length, slope of the catchment, baseflow, catchment area and the CN.
Most of these parameters (except for the CN value) were extracted from the DEM, using the spatial analysis tools of the ArcGIS Pro program. The values used to calibrate the model are presented in Table 3. The CN values are variable values depending on the land use area obtained for each year. Its values will be mentioned in the Results section. The lag time can either be specified directly or may be calculated from catchment characteristics using the standard SCS formula as in (Equation (14)) [39]: where: t l -lag time in hours; L-hydraulic length of the catchment (km); CN-SCS curve number; Y-average catchment slope (%). For each analyzed year, an event was isolated to show the evolution of runoff in relation to the precipitation ( Figure 10).  The average altitude of the basin and its location favor a long snow cover duration, usually from October to May, but sometimes it extends even into June. The high peak flow of March-April is well above the value of the simulated flow, due to the large amount of snow accumulated during the winter and which begins to melt in April. The situation is somewhat opposite starting from autumn, when the simulated flow values are higher than the observed ones. This situation may be explained by the large amount of precipitation falling and accumulated as a snow layer in the high mountains. The climatic parameter temperature plays an extremely important role in this case; therefore, the model cannot correctly simulate the events that took place during the winter season. Under these circumstances, for model calibration, we selected events only from the warm season (June to September) when the influence of snow melting and deposition is absent.
A difference can be observed between the simulated and the observed values ( Figure  10). The causes can be multiple, but the most important seem to be the uncertainties regarding the spatial and temporal variation in the precipitation value, as well as the initial humidity state of the basin; further, the structure of the model can be suspected as a cause, which is too simple. The generalized SCS model used is the only model that captures the influence of vegetation on runoff, using only two parameters for calibration (lag time, The average altitude of the basin and its location favor a long snow cover duration, usually from October to May, but sometimes it extends even into June. The high peak flow of March-April is well above the value of the simulated flow, due to the large amount of snow accumulated during the winter and which begins to melt in April. The situation is somewhat opposite starting from autumn, when the simulated flow values are higher than the observed ones. This situation may be explained by the large amount of precipitation falling and accumulated as a snow layer in the high mountains. The climatic parameter temperature plays an extremely important role in this case; therefore, the model cannot correctly simulate the events that took place during the winter season. Under these circumstances, for model calibration, we selected events only from the warm season (June to September) when the influence of snow melting and deposition is absent.
A difference can be observed between the simulated and the observed values ( Figure 10). The causes can be multiple, but the most important seem to be the uncertainties regarding the spatial and temporal variation in the precipitation value, as well as the initial humidity state of the basin; further, the structure of the model can be suspected as a cause, which is too simple. The generalized SCS model used is the only model that captures the influence of vegetation on runoff, using only two parameters for calibration (lag time, CN).
The daily precipitation value contained in the gridded dataset does not allow the detailed capture of the rainfall-runoff response. Although the shape of the simulated discharge hydrograph does not precisely follow the observed one, the overall accuracy of the amplitude and phase of the discharge hydrograph fall under the accepted and recommended thresholds, which confirms the reliability of the model for testing the proposed methodology. The main purpose of this analysis is to calculate the influence of the forested area on the runoff and not a rigorous analysis of the precipitation-runoff phenomenon.
In choosing the events under analysis, we did not necessarily take into account certain characteristics of the phenomenon to look for each year (amount of precipitation, similar duration, etc.) because the pre-event moisture conditions are different from case to case (e.g., the moisture condition prior to the occurrence of the rain event with a direct influence on soil permeability, air temperature, wind speed, etc.). Thus, the events themselves cannot be compared among them in terms of total rainfall amount over a period of time. What we followed in choosing the events was to identify events in which the simulated values are as close as possible to those observed.

Simulated Land Use Change Scenarios
In order to evaluate and understand the influence of reforestation and deforestation actions, it is necessary, first of all, to observe the effects that this change in land use has on the discharge values over time. As the precipitation amount is not the same each year, it is quite difficult to take a reference situation only based on precipitation values and forested areas (identifying two years with similar land use values and observing whether the discharge values are similar). However, because the volume of precipitation events is different, the runoff response or the volume of drained water is expected to be different.
We appreciate that the best method is to take into consideration three scenarios: scenario 1 (S1) is the initial one and by which we calibrate the model for each year event; scenario 2 (S2) assumes the hypothesis that the forested area is at the maximum potential area calculated by the FCD index method (the forested area determined in 1986 was 162 km 2 ); scenario 3 (S3), opposite to S2, refers to the minimum area calculated based on the FCD index and corresponds to that identified for the year 2003 (135 km 2 ) ( Figure 11). In choosing the events under analysis, we did not necessarily take into account certain characteristics of the phenomenon to look for each year (amount of precipitation, similar duration, etc.) because the pre-event moisture conditions are different from case to case (e.g., the moisture condition prior to the occurrence of the rain event with a direct influence on soil permeability, air temperature, wind speed, etc.). Thus, the events themselves cannot be compared among them in terms of total rainfall amount over a period of time. What we followed in choosing the events was to identify events in which the simulated values are as close as possible to those observed.

Simulated Land Use Change Scenarios
In order to evaluate and understand the influence of reforestation and deforestation actions, it is necessary, first of all, to observe the effects that this change in land use has on the discharge values over time. As the precipitation amount is not the same each year, it is quite difficult to take a reference situation only based on precipitation values and forested areas (identifying two years with similar land use values and observing whether the discharge values are similar). However, because the volume of precipitation events is different, the runoff response or the volume of drained water is expected to be different.
We appreciate that the best method is to take into consideration three scenarios: scenario 1 (S1) is the initial one and by which we calibrate the model for each year event; scenario 2 (S2) assumes the hypothesis that the forested area is at the maximum potential area calculated by the FCD index method (the forested area determined in 1986 was 162 km 2 ); scenario 3 (S3), opposite to S2, refers to the minimum area calculated based on the FCD index and corresponds to that identified for the year 2003 (135 km 2 ) ( Figure 11).

Delimitation of Forest-Covered Areas Using the FCD Index
The values of the FCD index for the analyzed years are presented in detail in Figure  12.

Delimitation of Forest-Covered Areas Using the FCD Index
The values of the FCD index for the analyzed years are presented in detail in Figure 12.

Evaluation of the Method for Determining the FCD Index
Analyzing Table 4, the low value of the kappa coefficient for 1985 is due to the much larger forest area in the period before 1989. Therefore, some points identified as being areas covered with pasture in 2018 were identified as areas covered with forests in 1985. A quite similar situation occurred for 2015, but this time, the error is also for forests found as pastures and pastures as forests.
The assessment of the method's accuracy was conducted by comparing the extension of the calculated forest area with the values observed on the orthophoto images from 2018.
The result of this analysis is presented in Table 4.   Table 4, the low value of the kappa coefficient for 1985 is due to the much larger forest area in the period before 1989. Therefore, some points identified as being areas covered with pasture in 2018 were identified as areas covered with forests in 1985. A quite similar situation occurred for 2015, but this time, the error is also for forests found as pastures and pastures as forests. The assessment of the method's accuracy was conducted by comparing the extension of the calculated forest area with the values observed on the orthophoto images from 2018.
The result of this analysis is presented in Table 4.

The Impact of Changes in Forested Areas on Peak Flood Discharges
As described in the methodology chapter, we used three scenarios to assess the impact of forest vegetation on runoff.
By modifying the values of the forested areas, according to S2 and S3, the values of the average CN will also change. The change in average CN values leads to changes in the estimated volume and flow values.
The results of the simulations for the three scenarios are presented in Figure 13 and Table 5. Once the model is set and calibrated, according to the methodology described in Section 2, the runoff hydrograph is analyzed for each year event ( Figure 14).     The results of the simulations on the three scenarios show an important influence of forests in reducing the peak flows. This influence is preserved up to a fairly significant share of the drained volume during floods (Table 5 and Figure 14). Even if the differences between the forest surfaces determined for different years and the values S2 and S3 are not very high, analyzing Table 5, we can observe an influence exerted on the peak flow, varying between −5.28% (S2) and 8.09% (S3).
The influence of forests in attenuating flood peaks is visible. The weight of this influence can be determined by drawing the correlation graph based on the surface differences of the three scenarios with the cumulated runoff volumes on each scenario: Forest area S1-forest area S2 compared with total volume S1-total volume S2 ( Figure 14a); Forest area S1-forest area S3 compared with total volume S1-total volume S3 (Figure 14b). The result of this analysis is shown in the two graphs of Figure 15. The results of the simulations on the three scenarios show an important influence of forests in reducing the peak flows. This influence is preserved up to a fairly significant share of the drained volume during floods (Table 5 and Figure 14). Even if the differences between the forest surfaces determined for different years and the values S2 and S3 are not very high, analyzing Table 5, we can observe an influence exerted on the peak flow, varying between −5.28% (S2) and 8.09% (S3).
The influence of forests in attenuating flood peaks is visible. The weight of this influence can be determined by drawing the correlation graph based on the surface differences of the three scenarios with the cumulated runoff volumes on each scenario: Forest area S1-forest area S2 compared with total volume S1-total volume S2 ( Figure  14a); Forest area S1-forest area S3 compared with total volume S1-total volume S3 ( Figure  14b).
The result of this analysis is shown in the two graphs of Figure 15.
(a) (b) Figure 15. Comparative analysis of the total volume of water drained during the events and the differences in forested areas, between the three scenarios S1, S2 and S3. (a) S1 compared with S2; (b) S1 compared with S3.
In the focus region, the catchment response to the increase in the forested surface was the reduction in the total runoff volumes for the selected event (Figure 15a). In the case of a decrease in the forested area (scenario S3), the total runoff increased during the analyzed event (Figure 15b). Figure 15. Comparative analysis of the total volume of water drained during the events and the differences in forested areas, between the three scenarios S1, S2 and S3. (a) S1 compared with S2; (b) S1 compared with S3.
In the focus region, the catchment response to the increase in the forested surface was the reduction in the total runoff volumes for the selected event (Figure 15a). In the case of a decrease in the forested area (scenario S3), the total runoff increased during the analyzed event (Figure 15b).
A similar situation was detected when comparative analysis was performed for the maximum runoff during the considered events ( Figure 16). In this case, the analyses were conducted considering the following differences: -Difference between forest area S1 and forest area S2 was compared with the difference between maximum runoff S1 and maximum runoff S2 ( Figure 16a); -Difference between forest area S1 and forest area S3 was compared with the difference between maximum runoff S1 and maximum runoff S3 (Figure 16b). A similar situation was detected when comparative analysis was performed for the maximum runoff during the considered events ( Figure 16). In this case, the analyses were conducted considering the following differences: -Difference between forest area S1 and forest area S2 was compared with the difference between maximum runoff S1 and maximum runoff S2 ( Figure 16a); -Difference between forest area S1 and forest area S3 was compared with the difference between maximum runoff S1 and maximum runoff S3 (Figure 16b).
(a) (b) Figure 16. Comparative analysis for the maximum value of water flow during events and differences in forested areas, between the three scenarios S1, S2 and S3. (a) S1 compared with S2; (b) S1 compared with S3.
A negative value of differences was obtained in all the graphs presented in Figures  15 and 16, where an increase in the forested area corresponds to a decrease in the maximum flow and vice versa. Additionally, the Pearson squared correlation index R 2 in both Figure 16. Comparative analysis for the maximum value of water flow during events and differences in forested areas, between the three scenarios S1, S2 and S3. (a) S1 compared with S2; (b) S1 compared with S3.
A negative value of differences was obtained in all the graphs presented in Figures 15 and 16, where an increase in the forested area corresponds to a decrease in the maximum flow and vice versa. Additionally, the Pearson squared correlation index R 2 in both Figures 15a and 16a (i.e., the relation between S1 and S2) indicated a stronger correlation compared to the relations S1-S3 (Figures 15b and 16b). Given the similar direction for both types of relations, S1-S2 and S1-S3, regardless of whether they are about the volume or flows, one can choose the best correlation. The correlation index presented in Figure 16a has the highest value, which is why we can use it as an equation in estimating the share of influence of forests on runoff (Equation (15): where: y-the estimate percentage of the maximum flow; x-the percentage of the forested area related to the reference value S1; 0.3717-slope of the regression line; −0.001-Y axis intercept.

Limitations
The analysis of the vegetation indices in order to determine the forest surface was conducted on satellite images, which have a low spatial resolution (30 m), and this can influence the results. However, the accuracy analysis of the method was performed on the orthophoto images from 2018, which have a 5 cm/pixel resolution. An analysis conducted on satellite images with a better spatial accuracy would lead to a much more accurate result, but in this paper, we aimed to find the order of magnitude of this influence.
Another limitation of the method is the quality and time step of precipitation data. We used daily gridded datasets with a spatial resolution of 0.1 • . For this type of analysis on small catchments, it is recommended to use hourly precipitation because of the short concentration time and to perform a comparative analysis of the data contained within the observed values [27] to determine the quality of the data contained in the gridded dataset. However, hourly data are not available for the area under study.
The model used in these analyses is designed to simulate events in which the snow parameter has no influence. Therefore, in choosing the event for simulation, it is essential that the snow or melting snow does not exert any influence.

Discussion
LULCC represents one of the most important factors that alter natural runoff regimes. The uncertainties in many studies were linked to the type of land use change. Among all types of land use, forests seem to play a role of crucial importance in the soil-water-atmosphere system [42]. Previous studies developed worldwide considered that the most important influence must be attributed to forested areas. Thus, a recent study in a mountain region in Brazil stated that "an increase in forest area is expected to decrease surface runoff and the associated water yield, and may mitigate the impacts of extreme precipitation events" [7]. Several other studies concluded that land use plays an essential role in determining the amount of percolating water, surface runoff and the lost water to the atmosphere through evapotranspiration [42,43].
In this research, we tried to capture and quantify the importance of forests. To achieve this objective, the study area was chosen taking into account the idea that the catchment should be small and controlled hydrometrically and have significant areal afforestation and deforestation change over time.
In fact, our correlations underlined that the forested area decrease contributed to the growth in the runoff water volume and discharge. Moreover, our results emphasize the importance of forests in the presented phenomena, where a forested area increase corresponds to a decrease in maximum flow, in agreement with the current literature findings, which reveal that deforestation leads to an immediate increase in the total discharge and the maximum flow [43].
In Romania, the areas covered with forests and their density have evolved over the years. Before 1989 (the year when the communist system collapsed in Eastern Europe), more attention was paid to forests to increase their economic efficiency; after 1989, there was a series of deforestation actions, both natural and anthropogenic (most of them illegal) [35].
After 2003, a new forest law was implemented, which stopped the advance of illegal deforestation, and after 2015, with the implementation of the "forest radar" application, this measure was intensified. Throughout the period, however, the administrator of the forest ROMSILVA tried numerous actions to reforest the deforested slopes. The results of these administrative measures are slowly beginning to show their results through a steady increase in forested areas.

Conclusions
Forested areas in river basins play an imperative role in attenuating flood peaks. One of the questions that this article wanted to answer is the quantification of this influence. This objective is covered by the PannEx project Flagship Question 4 [7]. According to the PannEx White Book, it is crucial to assess all anthropogenic influences which could alter the original natural water cycle, which, in part, represents the only true image regarding the available water resources [44].
The use of satellite images can be a valuable source of information regarding the extent of forested areas. Especially in mountainous areas, but not only in these areas, it is vital to consider that the entire land is not suitable for afforestation, either due to the vertical vegetation steps (alpine vegetation, areas of juniper and dwarf trees, moorland, etc.) or due to steep slopes or bare terrain. This type of analysis, based on successive time series assessments, can help to determine the maximum area that is suitable for afforestation.
Another advantage of this method is the use of freely available gridded datasets, meaning they can be used by any authority in order to evaluate the possible impact the extension or reduction in the forest surface would have.
The delineation of the forested area for certain regions (and at a certain date) can often be difficult to achieve for various reasons. Often the plans or maps held by the forest administrator do not provide a complete or real picture of the forested area, or the magnitude of deforestation or afforestation. Extraction of the FCD vegetation index based on Landsat images can be an important and relevant source of information regarding the delineation of the forested area. The results obtained in this research are encouraging and recommend the use of such a method for delimiting forested areas.
The analysis of the model results, performed on three scenarios, allows us to deduce a simple equation for estimating the influence of the forest surface on the maximum river flows. The use of such an equation allows us to hypothetically determine the percentage by which the maximum flow rate increases or decreases depending on the decrease or increase in the forest area.
In Romania, over recent decades, a lot of money has been invested in large-scale water works regarding structural measures for the attenuation of flood waves by dam constructions, accumulations, non-permanent accumulations, etc. This type of analysis revealed that for certain areas of the catchment, reforestation measures are enough to reduce the flood peaks, or to remove the negative effect of floods. Thus, by using this method, a lot of money and effort can be saved and can be invested in the quality of people's lives.
As further work, to test this method for areas with lower slopes, or hill or plain areas, identifying where the forest can play an even more important role is a priority. Furthermore, using the satellite images provided by the Copernicus Sentinel missions, with a much better resolution, could improve the results.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/atmos12070817/s1. Ref. [41] is also cited in supplementary materials. Data Availability Statement: The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.