The Impact of Land Use Changes on Soil Erosion in the River Basin of Miocki Potok, Montenegro

Land use change in all river basins leads to changes in hydrologic response, soil erosion, and sediment dynamics characteristics. Those changes are often viewed as the main cause of accelerated erosion rates. We studied the impact of land use changes on soil erosion processes in one of the watersheds in Montenegro: the Miocki Potok, using this watershed as a pilot river basin for this area. We simulated responses of soil erosion processes by using a process-oriented soil erosion Intensity of Erosion and Outflow (IntErO) model, with different settings of land use for the years 1970, 1980, 1990, 2000, 2010, and 2020. The model provides fast, effective, and affordable insight into the effects of land use change on soil erosion processes. Testing of the applied procedures was important for the further establishment of watershed management methodologies at the national level, for the other 300 river basins of Montenegro. For the current state of land use, calculated peak discharge for the Miocki Potok was 364 m3 s−1 (2020)–372 m3 s−1 (1970) for the incidence of 100 years, and there is a possibility for large flood waves to appear in the studied basin. Real soil losses, Gyear, were calculated at 13680 m3 year−1 (2020) and specific 333 m3 km−2 year−1 (2020). A Z coefficient value of 0.439 (2020) indicated that the river basin belongs to destruction category III. The strength of the erosion process was medium, and according to the erosion type, it was mixed erosion. According to our analysis, the land use changes in the last 50 years influenced a decrease in the soil erosion intensity for 14% in the Miocki Potok River Basin. Further studies should be focused on the detailed analysis of the land use changes trends with the other river basins at the national level, closely following responses of soil erosion to the changed land use structure, and effects of plant-and-soil interaction on soil erosion and sediment dynamics.


Introduction
Soil is one of the most important and most complex natural resources, but current developments (climate change, soil erosion, and urbanization) increasingly threaten this valuable resource. Soil degradation caused by erosion and rapid population increase is ranked among the most important environmental problems in the world [1][2][3]. Erosion is a key driver of land degradation, heavily affecting sustainable land management in various environments worldwide [4][5][6][7][8]. Human activities may have an important impact on the Earth's surface processes and contribute to changing the Earth's surface significantly since soil erosion and sediment dynamics are expected to change accordingly.
The negative effects of soil erosion on agricultural productivity and ecosystem are almost invisible over a shorter period of time [9,10]. Under normal conditions, soil erosion is a constant, progressing phenomenon, and, if no accelerated land degradation process such as gullying or landsliding occurs, we commonly are aware of the consequences only when it could be too late or too expensive to solve it. Conditions of increasing soil erosion not only affect agriculture but may dramatically reduce the life expectancy of reservoirs. In the long term, such processes may lead to changes in the provision of ecosystem services and to deplete fundamental natural resources.
The analysis based on gauging data can contribute to identifying the effects of land use changes to soil erosion and sediment dynamics, but these data are available only for a small number of medium to large catchments, whereas for small basins and a variety of landscapes, this kind of information is still missing. Thus, quantification of the responses to soil erosion and sediment yield to land use variations can only be achieved by numerical simulation of soil erosion processes [11].
In their 2009 study, Volk [12] encouraged researchers to conduct such analysis highlighting the importance of achieving the aims of the European Water Framework Directive. According to this directive, land use and land management options are to be used as tools for water quantity and quality control in order to achieve and maintain ecologically stable and productive water bodies. This idea on policy issues has been well received among researchers in Montenegro and Serbia, taking into consideration the current EU accession agenda of these two countries.
Numerous semi-quantitative methods are at the moment in use for assessing sediment yield at the catchment scale. These methods take into account parameters such as geology, soil, geomorphology, topography, climate, rainfall erosivity, and land cover [13][14][15].
The objective of this research was to simulate responses of soil erosion processes with the different settings of the land use in the River Basin of Miocki Potok in Montenegro over a 50 year time interval (1970-2020) by using a process-oriented soil erosion Intensity of Erosion and Outflow (IntErO) model [13]. The model provides fast, effective, and affordable insight into the effects of land use change on soil erosion processes and sediment yield. Simulation results under different scenarios of land use change were used to analyze the response characteristics of soil erosion processes to the variations in land use with the aim of optimizing current and future land use options and to predict or optimize future land and water management options in the wider region.
It is well known that soil erosion processes can be, to some extent, controlled by appropriate land management. This requires the collection of field data and the use of predictive models for the evaluation of different management scenarios for the protection of soils. Field measurements of erosion and sedimentation using classical techniques are time-consuming and expensive; therefore, modeling of the erosion process represents a valid alternative to adequately predict both the runoff and soil loss.
According to the various reports of the European Soil Bureau Institute for Environment and Sustainability and the European Environment Agency, the Universal Soil Loss Equation model (USLE) is extensively used in the following European countries: Austria, Bosnia, and Herzegovina (including the Republic of Srpska), Bulgaria, Italy, Hungary, Norway, Romania, Slovakia, Finland, Czech Republic, Spain, and Switzerland. The Revised Universal Soil Loss Equation (RUSLE) is used in Belgium; the UK, Germany, and France are using their domestic/national models.
Assessing the sustainability of water supply is of key importance, especially in Montenegro, the wider Balkans, and the Mediterranean. Montenegro has been identified as a hot spot in regards to the impact of climate change and as a region in which water demand is continuously growing for population needs, tourism development, and irrigated areas expansion.
Climate projections in some scenarios predict an increase in temperature of up to 2 • C and a decrease in daily precipitation by 20%. In this context, there is a great need for experts in the Balkans to assess future changes in water distribution capacities in the basins, taking into account the newly emerging climate change, management of reservoirs and dams. It is necessary to work in a broader context on the preparation of scenarios for the use of available water.
This research attempts to establish one of the sustainable forms of modeling that would be calibrated and validated in the Polimlje region of the River Lim Basin and further used to assess how water resources can meet water requirements for local, agricultural, and eco sectors.

Study Area
Northeastern Montenegro is mainly mountainous, with the presence of deep valleys incised into limestone ranges. Other parts of the northeast are rather hilly and underlain by Palaeozoic rocks. In this region, the highest peaks of Montenegro are found, including Komovi (2487 m a.s.l.) and Zla Kolata (2535 m a.s.l.) in the Prokletije Mountains. Nevertheless, this area is quite densely populated and includes the towns of Plav, Andrijevica, Berane, and Bijelo Polje. The rivers in this region drain to the Black Sea, and some of them form deep canyons crossing limestone formations. Further downstream, they form broader valleys flowing through softer Palaeozoic material [13].
This study was conducted in the area of the river basin of Miocki Potok ( Figure 1 Figure 1). This watershed encompasses an area of 41 km 2 , and the natural length of the main watercourse, Lv, is 6.4 km.
A flat area occurs around the village of Kanje in the lower reach, and steep slopes make up the upper part of the watershed. The average gradient of the catchment, Isr, is 31.6% and indicates that very steep slopes are present in the river basin. The average river basin altitude, Hsr, is 917.83 m a.s.l; the average elevation difference of the basin, D, is 382.83 m. The drainage density is low indicating a rather permeable substrate.
On the basis of more than 70 years (1948-2020) of temperature and precipitation data, the study area is characterized by a mountain Mediterranean climate with rainy autumns and springs, cold winters, and a deficit of precipitation in the summer months. Basic data on the area needed for the calculation of soil erosion intensity and runoff are presented in Table 1.  The absolute maximum air temperature ever recorded was 39.2 • C. Winters are severe, with negative temperatures as low as −27.6 • C. The average annual air temperature, t0, was 8.9 • C. The average annual precipitation, Hyear, was 873 mm. The temperature coefficient for the region, T, was calculated at 0.99. The torrential rain, hb, was calculated at 84.7 mm.

The Geological Structure and Soils of the Area
Montenegro is a part of the Dinaric Alps, which are, in turn, parts of the complex thrust-and-fold system of the Mediterranean area. This area represents a proxy of long-lasting interactions between Eurasia and Gondwana, resulting in a system of fold-and-thrust belts and associated foreland and back-arc basins. The system cannot be interpreted as the end product of one single Alpine orogeny as the major suture zones result from various tectonic events that closed different oceanic basins of the former Tethyan Ocean. The Dinarides-Albanides-Hellenides orogenic belt is caused by a Tertiary collision between the Adriatic promontory and the Serbo-Macedonian-Rhodope blocks. The belt is bordered to the west by a foreland basin in the Eastern Adriatic basin filled with Eocene-Quaternary deep marine sediments. The study area at large consists of various types of sediment, magmatic and metamorphic rocks generated in the long, Palaeozoic to Quaternary, interval. Most of the terrain is underlain by Mesozoic formations of carbonate composition, while magmatic and silico-clastic rocks are substantially less present. Paleozoic geological formations consist of sedimentary and metamorphic, silico-clastic rocks found mostly in the north-eastern parts of Montenegro, while Cainozoic rocks of carbonate and clastic composition occur almost in all regions of Montenegro.
The main rocks outcropping in the area are clastic and subordinate carbonate rocks from the Paleozoic, Triassic clastites, volcanites, tuffs, limestone and dolomites, Jurassic clastic rocks with diabasic effusions and metamorphic rocks and Quaternary, mainly alluvial and colluvial deposits.
In order to define the permeability of the rocks of the study area, we used the Geological Atlas of Serbia [37] and extracted a geological map of the Miocki Potok watershed from the Geological map of Montenegro [38]. The basin of the Miocki Potok consists of Devonian-Carboniferous (D+C) and Permian (P) phyllites, argyllo-phyllites, metasandstones, and conglomerates.
The coefficient of the region's permeability, S1, was calculated to be 0.91. Within the study basin, the area underlain by limestone is very permeable (class fp, 10%), and the rest has poor permeability (class fo, 90%).
According to the results of the field visits and supplementary laboratory analysis, but also using the previous research data of the project Soils of Montenegro  carried out by the team of the Biotechnical institute of the University of Montenegro [39] and Spalevic [13], the most common soil types in the study basin were: Dystric Cambisols and Fluvisols and Colluvial Fluvisols in the lower alluvial plain.

Methods
Morphometric methods were used to determine slope gradient, specific lengths, aspect, and shape. The depth of the erosion base and the density of erosion rills were identified during field visits. Google Earth was used to further investigate the morphological features of the watershed with a realistic 3D view. In that manner, some minor sites or phenomena that were not reached during the fieldwork were also taken into consideration. According to Chang et al. [40], Google Earth is considered more than adequate to identify and to characterize soil erosion processes.
Some pedological profiles were excavated, and soil samples were collected for physical and chemical analysis. The grain-size distribution of the soils was determined by the pipette method [41]; the soil samples were air-dried at 105 • C and dispersed using sodium pyrophosphate. The soil reaction (pH in H 2 O and nKCl) was determined with a potentiometer. Total carbonates were determined by the volumetric Scheibler method [42]; the content of the total organic matter was determined by the Kotzman method [43]; easily accessible phosphorous and potassium were determined by the Al-method [44], and the adsorptive complex (y1, S, T, V) was determined by the Kappen [45] method.
Blinkov and Kostadinov [46] evaluated the applicability of various erosion risk assessment methods for engineering purposes. The factors to be considered depends on the scale, various erosion tasks, and local specific needs. According to their findings, the Erosion Potential Method-EPM [16] was the most suitable at the catchment level for watershed management in southeastern European countries. The EPM is distinguished by its high degree of reliability in calculating sediment yields, suspended sediment yield, and reservoir sedimentation [47].
In Montenegro, the EPM proved to return good predictions of runoff and soil erosion intensity, specifically in the Region of Polimlje [13,[48][49][50]. In this study, the Intensity of Erosion and Outflow (IntErO) program package [13] was used to obtain estimates of soil erosion intensity in the river basin of Miocki Potok. The method has clearly defined procedures, and subjective evaluations are reduced to a minimum. The EPM is embedded in the algorithm of the IntErO model.

The IntErO Model
The IntErO model uses the Erosion Potential Method in its algorithm background [13,26]. The IntErO, an upgrading of the River Basins [48,51] and the Surface and Distance Measuring [48] programs, is simple in handling and can be used to calculate a large number of data with the processing of 22 input parameters, returning, after the processing, 26 result parameters (Coefficient of the river basin form, A; Coefficient of the watershed development, m; Average river basin width, B; (A)symmetry of the river basin, a; Density of the river network of the basin, G; Coefficient of the river basin tortuousness, K; Average river basin altitude, Hsr; Average elevation difference of the river basin, D; Average river basin decline, Isr; The height of the local erosion base of the river basin, Hleb; Coefficient of the erosion energy of the river basin's relief, Er; Coefficient of the region's permeability, S1; Coefficient of the vegetation cover, S2; Analytical presentation of the water retention in inflow, W; Energetic potential of water flow during torrent rains, 2×gDFˆ1 2 ; Maximal outflow from the river basin, Qmax; Temperature coefficient of the region, T; Coefficient of the river basin erosion, Z; Production of erosion material in the river basin, Wyear; Coefficient of the deposit retention, Ru; Real soil losses, Gsp; Real soil losses per km 2 . The model considers six factors related to lithology (rocks permeability by percentage: fp, permeable; fpp, semipermeable; fo, low permeability) and soil type (erodibility coefficient, Y), topographic and relief data (I coefficient), monthly mean and annual precipitation (P coefficient), temperatures annual averages (t coefficient), land cover data (X coefficient), the state of erosion patterns, and development of the watercourse network (Φ coefficient). The IntErO model can be characterized as semi-quantitative because it is based on a combination of descriptive and quantitative procedures. Compared to other semi-quantitative methods, this is the most quantitative because it uses descriptive evaluation for three parameters only: soil erodibility, soil protection, and the extent of erosion in the catchment. All other parameters are quantitative catchment descriptors. IntErO flowchart is presented in Figure 2. The annual volume of soil detached due to soil erosion (Wyear) was calculated in m 3 km −2 y −1 by the following Equation (1): P is the annual average rainfall in mm.
T is the temperature coefficient, calculated from the following equation: t represents the mean annual temperature in °C, which is extracted from the meteorological data collected from the State Hydrometeorological Institute of Montenegro and the Biotechnical faculty of the University of Montenegro.
The Z coefficient describes the intensity of the erosive process, can be classified according to the degree of erosion (Table 2); it is calculated using Equation (3).

= . . ( + )
(3) The annual volume of soil detached due to soil erosion (Wyear) was calculated in m 3 km −2 y −1 by the following Equation (1): P is the annual average rainfall in mm.
T is the temperature coefficient, calculated from the following equation: t represents the mean annual temperature in • C, which is extracted from the meteorological data collected from the State Hydrometeorological Institute of Montenegro and the Biotechnical faculty of the University of Montenegro.
The Z coefficient describes the intensity of the erosive process, can be classified according to the degree of erosion (Table 2); it is calculated using Equation (3).
X represents the coefficient of soil protection (Table 3), it is a dimensionless parameter based on the vegetal cover and catchment's land use, and it varies from 0.05 to 1; values close to 0 indicate low soil protection while values close to 1 mean high soil protection [52]. Y is the coefficient of soil erodibility, which depends on the pedological and lithological characteristics of the watershed and indicates the resistance of soils to erosion. The values of this dimensionless factor can be determined either by laboratory experiments or by field measurements [48,52] and usually ranges from 0.25 to 2. Values close to 0.2 indicate low erodibility, whereas values close to 2 represent strongly erodible formations [52,53]. In this study, we used the soil map data and the geological map of the study area to evaluate the Y factor. Subsequently, coefficient values were assigned to each soil type based on the EPM guidelines and previous studies [52,[54][55][56][57].
ϕ is a coefficient that depends on the active erosion and on the degree of extension of linear erosion and mass movements. It is a dimensionless factor whose values range from 0.1 to 1 [58,59]. Each type of erosion form has taken on a value of Φ according to the guidelines of the EPM method [52,54,57]. In calculating Z, Y, Xa, and ϕ, the weighted average value was calculated all the polygons with various values [26].
Isr represents the average mean basin slope of the studied area expressed as a percentage. This parameter is extracted from the Digital Terrain Model of the catchment (DTM-ASTER-GDEM) with a horizontal resolution of 30 m and a vertical resolution of 20 m (NASA earth data), while F is a watershed area in km 2 .
The total volume of sediment produced in the different areas of the watershed does not fully reach downstream. A portion is redeposited in streams or other areas of the basin; therefore, it is essential to calculate the specific real sediment production (Gsp) in m 3 km −2 y −1 by the following equation: This is possible by multiplying the average annual production of sediments Wyear by a delivery ratio coefficient Ru, as calculated in Equation (5).
Lv represents the length of the main river in Km, O is the perimeter of the catchment area in Km, and D is the average height of the basin in relation to the closing section [26].

IntEro Model Verification
In order to support the idea of using the model for watershed management planning in the study region, it was necessary to perform model validation and verification. To validate the sediment yield calculations for the Polimlje region, which includes the Miocki Potok basin, bathymetric measurements were carried out in the Potpec reservoir on the river Lim in 2017 to calculate the accumulation of sediment in the 1967-1976; 1967-1981; 1967-1986; 1967-1991; 1967-1999; 1967-2006; 1967-2012; 2012-2017 intervals.
The "Potpec" reservoir is about 17 km long, with a maximum depth of about 30 m close to the river dam. Over the years, it has been receiving large amounts of sediment that have reduces its useful volume. In order to determine the quantities of deposited sediment, a survey of the bottom of the Lim riverbed was organized in the period of 11-26 April 2017.
Field measurements of the accumulation in 2012 were performed with professional hydrographic equipment. We used a GPS rover, which provides measurement accuracy of 10 mm horizontally, and 20 mm vertically, and a Trimble R6 base with an internal radio modem, while the depth in the accumulation was measured with a single-frequency portable echo sounder Oda Hydrotrac which provides measurement accuracy of 1 cm. The equipment is installed on a boat that is structurally adapted for the mentioned type of such research work. The trigonometric points closest to the accumulation on the left and right banks were used as the geodetic basis for the survey.
After the field recording, we developed a digital terrain model of the bottom of the accumulation. From the digital terrain model, transverse profiles are drawn whose position is defined by coordinates.

Rainfall Change
Rainfall and overland flow are key factors in soil particle detachment, entrainment, and transport. Overland flow depth depends on several factors, including soil field capacity, slope gradient, and the volume and intensity of precipitation. The rainfall factor is, in fact, included in almost all the most renowned models to predict soil erosion and sediment yield, such as the Modified Universal Soil Loss Equation (MUSLE), the Unit Sediment graph (USG), the Water Erosion Prediction Project (WEPP), etc.. Rainfall is the soil erosion triggering factor and, after land use/cover, it is the most variable, especially over short timescales.
In the study area, long time series of rainfall data (1950-2016) is available only for two meteo-stations: Berane and Bijelo Polje. Soil erosion is not influenced by the rainfall amount and mostly by rainfall intensity. For this study, unfortunately, only daily precipitation, for Bijelo Polje (1966-2019), and maximum snow thickness and duration data, For Berane, Bijelo Polje, and Plav, are available. Langbein and Schumm [60] were the first to investigate the relationship between mean annual precipitation and sediment yield. Since that pioneering paper, very many other studies on the prediction of sediment yield followed, and many of them used annual precipitation as an entry parameter [61,62]. The variability of annual precipitation for Berane and Bijelo Polje is reported in Figure 3. Berane does not show any trend, whereas the Bijelo Polje diagram is characterized by a moderate increasing trend of 99 mm across a 66-year interval. From the daily rainfall data set two time series were extracted: the mean monthly maximum daily intensity and the absolute highest daily intensity of each year. These two time series are reported in Figure 4, which shows clearly that there was no trend for the average monthly maximum daily intensity and a very moderate increase of only 9 mm in 24 hours across the 54−year interval for the annual peak values of daily intensity. The highest precipitation of 158 mm in one day was recorded in 1992; other peaks around 100 mm in 24 hours were recorded in 1972, 1986, and 1996. Runoff generation and volume are largely influenced by rainfall intensity; however, the threshold for soil particles detachment and beginning of sediment supply to the drainage network varies for different environments according to the climatic conditions and vegetation cover type and density. Though in his plot experiment, sediment yield was not measured [63], showed that cycles of intermittent rain with an intensity of 10 mm/h for periods of 3-25 minutes are able to generate an appreciable runoff of 1-7 mm depth.
Some photos in Dunkerley's paper [63] show evidence of soil saturation, and soil particle entrain met should be expected. Based on Dunkerley's experiments and on the author's experience of field observations during rainstorms, a daily intensity of 20 mm/24h was considered as an appropriate threshold for the generation of suspended sediment transport events in the study catchment. In order to surrogate the information about rainfall intensity at an hourly scale, we built a time series of the number of days with daily rainfall higher than 20 mm and investigated the occurrence of any trend ( Figure 5). From the daily rainfall data set two time series were extracted: the mean monthly maximum daily intensity and the absolute highest daily intensity of each year. These two time series are reported in Figure 4, which shows clearly that there was no trend for the average monthly maximum daily intensity and a very moderate increase of only 9 mm in 24 hours across the 54-year interval for the annual peak values of daily intensity. The highest precipitation of 158 mm in one day was recorded in 1992; other peaks around 100 mm in 24 hours were recorded in 1972, 1986, and 1996.  From the daily rainfall data set two time series were extracted: the mean monthly maximum daily intensity and the absolute highest daily intensity of each year. These two time series are reported in Figure 4, which shows clearly that there was no trend for the average monthly maximum daily intensity and a very moderate increase of only 9 mm in 24 hours across the 54−year interval for the annual peak values of daily intensity. The highest precipitation of 158 mm in one day was recorded in 1992; other peaks around 100 mm in 24 hours were recorded in 1972, 1986, and 1996. Runoff generation and volume are largely influenced by rainfall intensity; however, the threshold for soil particles detachment and beginning of sediment supply to the drainage network varies for different environments according to the climatic conditions and vegetation cover type and density. Though in his plot experiment, sediment yield was not measured [63], showed that cycles of intermittent rain with an intensity of 10 mm/h for periods of 3-25 minutes are able to generate an appreciable runoff of 1-7 mm depth.
Some photos in Dunkerley's paper [63] show evidence of soil saturation, and soil particle entrain met should be expected. Based on Dunkerley's experiments and on the author's experience of field observations during rainstorms, a daily intensity of 20 mm/24h was considered as an appropriate threshold for the generation of suspended sediment transport events in the study catchment. In order to surrogate the information about rainfall intensity at an hourly scale, we built a time series of the number of days with daily rainfall higher than 20 mm and investigated the occurrence of any trend ( Figure 5). Runoff generation and volume are largely influenced by rainfall intensity; however, the threshold for soil particles detachment and beginning of sediment supply to the drainage network varies for different environments according to the climatic conditions and vegetation cover type and density. Though in his plot experiment, sediment yield was not measured [63], showed that cycles of intermittent rain with an intensity of 10 mm/h for periods of 3-25 minutes are able to generate an appreciable runoff of 1-7 mm depth.
Some photos in Dunkerley's paper [63] show evidence of soil saturation, and soil particle entrain met should be expected. Based on Dunkerley's experiments and on the author's experience of field observations during rainstorms, a daily intensity of 20 mm/24 h was considered as an appropriate threshold for the generation of suspended sediment transport events in the study catchment. In order to surrogate the information about rainfall intensity at an hourly scale, we built a time series of the number of days with daily rainfall higher than 20 mm and investigated the occurrence of any trend ( Figure 5). The data indicated no trend, but the diagram in Figure 5 shows that in the mid-late 1970s, more than 10 days in a year with rainfall depth of 200 mm or more were rather common. From 1990 to 2006, the highest number of days (19) and several years with 13-16 days are recorded. In the last decade (2007-2019), most of the data plotted below ten days and lower than average sediment yield values should be expected. In order to further investigate the role of annual precipitation on sediment yield of the study catchment, the Modified Fournier Index-MFI [64] was also calculated as: In which p is the amount of rain in a month (mm), and P is the annual precipitation (mm). The MFI has been developed and used as a surrogate rainfall erosivity index, and it is classified into five classes in which the values less than 90 indicate low erosivity; between 90 and 16 indicate moderate erosivity, and values beyond 160 indicate high and very high erosivity.
The time series of the MFI calculated for Berane and Bijelo Polje across the 1950-2016 interval are reported in Figure 6. This diagram shows that most of the MFI values were included in the moderate rainfall erosivity class and that no significant trend was present. High erosivity peaks occurred only in 1974, 1985, and in 1992. It is worth noticing that the latter two peaks corresponded to the peaks of absolute maximum daily intensity ( Figure 6). This index is used to represent the rainfall erosivity. The dashed lines indicate the transition between low to moderate and from moderate to high erosivity. The solid black lines are the trend lines; their interpolating equations are also reported. The data indicated no trend, but the diagram in Figure 5 shows that in the mid-late 1970s, more than 10 days in a year with rainfall depth of 200 mm or more were rather common. From 1990 to 2006, the highest number of days (19) and several years with 13-16 days are recorded. In the last decade (2007-2019), most of the data plotted below ten days and lower than average sediment yield values should be expected. In order to further investigate the role of annual precipitation on sediment yield of the study catchment, the Modified Fournier Index-MFI [64] was also calculated as: In which p is the amount of rain in a month (mm), and P is the annual precipitation (mm). The MFI has been developed and used as a surrogate rainfall erosivity index, and it is classified into five classes in which the values less than 90 indicate low erosivity; between 90 and 16 indicate moderate erosivity, and values beyond 160 indicate high and very high erosivity.
The time series of the MFI calculated for Berane and Bijelo Polje across the 1950-2016 interval are reported in Figure 6. This diagram shows that most of the MFI values were included in the moderate rainfall erosivity class and that no significant trend was present. High erosivity peaks occurred only in 1974, 1985, and in 1992. It is worth noticing that the latter two peaks corresponded to the peaks of absolute maximum daily intensity ( Figure 6).
Water 2020, 12, x FOR PEER REVIEW 11 of 28 Figure 5. Time series of the number of days with rainfall higher than 20 mm in a year. The red line discriminates the years with more than ten days with a daily intensity of 20 mm or more.
The data indicated no trend, but the diagram in Figure 5 shows that in the mid-late 1970s, more than 10 days in a year with rainfall depth of 200 mm or more were rather common. From 1990 to 2006, the highest number of days (19) and several years with 13-16 days are recorded. In the last decade (2007-2019), most of the data plotted below ten days and lower than average sediment yield values should be expected. In order to further investigate the role of annual precipitation on sediment yield of the study catchment, the Modified Fournier Index-MFI [64] was also calculated as: In which p is the amount of rain in a month (mm), and P is the annual precipitation (mm). The MFI has been developed and used as a surrogate rainfall erosivity index, and it is classified into five classes in which the values less than 90 indicate low erosivity; between 90 and 16 indicate moderate erosivity, and values beyond 160 indicate high and very high erosivity.
The time series of the MFI calculated for Berane and Bijelo Polje across the 1950-2016 interval are reported in Figure 6. This diagram shows that most of the MFI values were included in the moderate rainfall erosivity class and that no significant trend was present. High erosivity peaks occurred only in 1974, 1985, and in 1992. It is worth noticing that the latter two peaks corresponded to the peaks of absolute maximum daily intensity ( Figure 6). This index is used to represent the rainfall erosivity. The dashed lines indicate the transition between low to moderate and from moderate to high erosivity. The solid black lines are the trend lines; their interpolating equations are also reported. This index is used to represent the rainfall erosivity. The dashed lines indicate the transition between low to moderate and from moderate to high erosivity. The solid black lines are the trend lines; their interpolating equations are also reported.

Demographic Changes in the River Basin of Miocki Potok
According to the 2011 census, the municipality of Bijelo Polje, covers an area of 923 km 2 , with 13,199 households, 137 settlements, and 46,051 inhabitants, which became 42,191 in 2020 (Figure 7). The borders of the Municipality of Bijelo Polje are also largely the borders of Montenegro to Serbia, and they extend to the tops of the mountains. In the middle of the past century, 1267 inhabitants lived in the catchment area of the Miocki Potok that is part of the Bijelo Polje municipality (Figure 7). The average density of the population was 30.9 inhabitants per km 2 , which was slightly less than the density of the population of the municipality of Bijelo Polje (39.8).
Until the 1970s, population growth was registered in all settlements of Miocki Potok basin. Index rates ranged from 101.  The population decline since the 1970s was caused by industrialization, which initially had positive effects, as the villages were relieved of latent redundancies. Workers moved primarily went to the city center of Bijelo Polje and later to other cities in Montenegro and beyond (especially to other republics of former Yugoslavia). However, as the industrialization process intensified, the countryside was increasingly neglected, and the rural exodus, instead of declining, increased over time, leading to a significant depopulation of all rural settlements, not only in the municipality of Bijelo Polje but also throughout the Northern region of Montenegro.
Migration regained intensity at the beginning of the 21st century. In many villages, the birth rate decreased to a minimum or to zero, which led to a lack of labor force and economic decline of the local communities. The emigration regarded mainly the young who left behind villages consisting of elderly households. With such a marked change in the age composition of the villages and equivalent marked change in the land use/cover was inevitable.

Vegetation and Land Use
In Montenegro, vegetation distribution is strongly influenced by (i) the distance to the Mediterranean Sea (continental or maritime climate) and (ii) the altitude (temperate or mountain climate). Due to Montenegro's rugged topography, there are large differences in vegetation over short distances [65]. The Miocki Potok river basin is located in the Dinarid Province of the middle-south-eastern European mountainous biogeographical region. Forests dominate this river basin accounting for 60% of the total vegetation cover.
Most of the river basin is covered by degraded beech forests (Fagetum montanum). The degraded forests are located near settlements and roads because of the common practice of firewood harvesting. These forests are characterized by a terminate canopy and by a large number of species of ground flora, shrubs, and lower trees. They differ from beech forests in the inner parts of the basin, especially from the sub-association Fagetum montanum typicum, which is characterized by a dense canopy. Near the villages of Milovo and Dobrakovo, especially on southern exposures, there are forests of Sessile oak and Turkish oak (Quercetum petraeae cerridis Lak.). In the lower part of the basin, a narrow belt along the river channel is covered with hydrophilic forest (Alnetea glutinosae, Salicetea herbacea).
The State-owned forest portions have a much better structure than the private forests. We measured seven clusters within the basin. Position of clusters in the river basin study area presented in the Figure 9. consisting of elderly households. With such a marked change in the age composition of the villages and equivalent marked change in the land use/cover was inevitable.

Vegetation and Land Use
In Montenegro, vegetation distribution is strongly influenced by (i) the distance to the Mediterranean Sea (continental or maritime climate) and (ii) the altitude (temperate or mountain climate). Due to Montenegro's rugged topography, there are large differences in vegetation over short distances [65]. The Miocki Potok river basin is located in the Dinarid Province of the middle-south-eastern European mountainous biogeographical region. Forests dominate this river basin accounting for 60% of the total vegetation cover.
Most of the river basin is covered by degraded beech forests (Fagetum montanum). The degraded forests are located near settlements and roads because of the common practice of firewood harvesting. These forests are characterized by a terminate canopy and by a large number of species of ground flora, shrubs, and lower trees. They differ from beech forests in the inner parts of the basin, especially from the sub-association Fagetum montanum typicum, which is characterized by a dense canopy. Near the villages of Milovo and Dobrakovo, especially on southern exposures, there are forests of Sessile oak and Turkish oak (Quercetum petraeae cerridis Lak.). In the lower part of the basin, a narrow belt along the river channel is covered with hydrophilic forest (Alnetea glutinosae, Salicetea herbacea).
The State-owned forest portions have a much better structure than the private forests. We measured seven clusters within the basin. Position of clusters in the river basin study area presented in the Figure 9. Data obtained by field measurements showed the prevalence of degraded forests. The average timber volume in the forests is 111 m 3 ha −1, with an average of 327 trees per hectare. In five representative clusters, the dominant species is beech. Beech wood volume varies from 81.5 m 3 ha −1 to 331 m 3 ha −1 . In the Miocki Potok basin, we also recorded Betula verrucosa, Quercus cerris, Carpinus betulus, Acer campestre, Quercus petraea, and Prunus avium.
In recent years, the migration of the rural population to towns resulted in an increase in forest areas with decreasing processes of forest degradation and devastation.
In order to describe the impact of land use on soil erosion in the Miocki Potok catchment, we used a series of data from The First National Forest Inventory of Montenegro [66], Forest Management Plans, Cadastre, Landsat images, and Statistical Yearbooks to quantify variations in land use across six intervals : 1970, 1980, 1990, 2000, 2010, and 2020.
The data indicate that forests cover about half of the study basin and that the forest area has been slightly increasing during the last decades ( Figure 10). The findings of this study are in line with observations of Nyssen [67], who stated that there is marked forest regrowth in the former Yugoslavia. Overall, forest cover has been naturally increasing in this region. Data obtained by field measurements showed the prevalence of degraded forests. The average timber volume in the forests is 111 m 3 ha −1, with an average of 327 trees per hectare. In five representative clusters, the dominant species is beech. Beech wood volume varies from 81.5 m 3 ha −1 to 331 m 3 ha −1 . In the Miocki Potok basin, we also recorded Betula verrucosa, Quercus cerris, Carpinus betulus, Acer campestre, Quercus petraea, and Prunus avium.
In recent years, the migration of the rural population to towns resulted in an increase in forest areas with decreasing processes of forest degradation and devastation.
In order to describe the impact of land use on soil erosion in the Miocki Potok catchment, we used a series of data from The First National Forest Inventory of Montenegro [66], Forest Management Plans, Cadastre, Landsat images, and Statistical Yearbooks to quantify variations in land use across six intervals : 1970, 1980, 1990, 2000, 2010, and 2020.
The data indicate that forests cover about half of the study basin and that the forest area has been slightly increasing during the last decades ( Figure 10). The findings of this study are in line with observations of Nyssen [67], who stated that there is marked forest regrowth in the former Yugoslavia. Overall, forest cover has been naturally increasing in this region. This denser vegetation has led to higher water infiltration into the soil. The process of forest regrowth may be further accelerated by applying good silviculture measures and forest management practices. Meadows, pastures, and orchards cover around 40% of the river basin, but this land cover type has been decreasing in recent years. According to Grimes [68], some marginal lands were cultivated again in the first decade of this century, as a consequence of the economic crisis in the region.
This trend was also confirmed by our research on land use in the small Miocki Potok basin. For the period of 2000-2020, a decrease in the area of arable and cultivated land was recorded. Land use changes in the Miocki Potok basin are presented in Figure 11. This denser vegetation has led to higher water infiltration into the soil. The process of forest regrowth may be further accelerated by applying good silviculture measures and forest management practices. Meadows, pastures, and orchards cover around 40% of the river basin, but this land cover type has been decreasing in recent years. According to Grimes [68], some marginal lands were cultivated again in the first decade of this century, as a consequence of the economic crisis in the region.
This trend was also confirmed by our research on land use in the small Miocki Potok basin. For the period of 2000-2020, a decrease in the area of arable and cultivated land was recorded. Land use changes in the Miocki Potok basin are presented in Figure 11.  basin (1970, 1980, 1990, 2000, 2010, and 2020). Figure 11. Land use in the Miocki Potok river basin (1970,1980,1990,2000,2010, and 2020).
Within the basin area, degraded forests are the most widespread form of vegetation cover (29.0%, in 1990) and well-developed forests (28.3 in 2020). Other types are as follows: mountain pastures (21.5%, on average for the period 1970-2020), meadows (16.1% on average for the period 1970-2020), plow-lands (12.6% on average for the period 1970-2020), and orchards (less than 1%, on average for the period 1970-2020) ( Table 4). The coefficient of the vegetation cover, S2, was calculated to range from 0.72 to 0.74, the coefficient of the river basin planning, Xa, from 0.45 to 0.50 Figures 12 and 13.

Current Erosion and the Impact of Land Use on Soil Erosion Intensity
The dominant erosion form in the study area is from surface runoff, including severe forms of erosion such as rills, gullies, and ravines. Surface erosion has affected all soils and slopes, and it is most pronounced on steep slopes with scarce or degraded vegetation cover (Figure 14).  Within the basin area, degraded forests are the most widespread form of vegetation cover (29.0%, in 1990) and well-developed forests (28.3 in 2020). Other types are as follows: mountain pastures (21.5%, on average for the period 1970-2020), meadows (16.1% on average for the period 1970-2020), plow-lands (12.6% on average for the period 1970-2020), and orchards (less than 1%, on average for the period 1970-2020) ( Table 4). The coefficient of vegetation cover, S2, was calculated to range from 0.72 to 0.74, the coefficient of the river basin planning, Xa, from 0.45 to 0.50 Figures 12 and 13.

Current Erosion and the Impact of Land Use on Soil Erosion Intensity
The dominant erosion form in the study area is from surface runoff, including severe forms of erosion such as rills, gullies, and ravines. Surface erosion has affected all soils and slopes, and it is most pronounced on steep slopes with scarce or degraded vegetation cover (Figure 14).

Current Erosion and the Impact of Land Use on Soil Erosion Intensity
The dominant erosion form in the study area is from surface runoff, including severe forms of erosion such as rills, gullies, and ravines. Surface erosion has affected all soils and slopes, and it is most pronounced on steep slopes with scarce or degraded vegetation cover (Figure 14). Within the basin area, degraded forests are the most widespread form of vegetation cover (29.0%, in 1990) and well-developed forests (28.3 in 2020). Other types are as follows: mountain pastures (21.5%, on average for the period 1970-2020), meadows (16.1% on average for the period 1970-2020), plow-lands (12.6% on average for the period 1970-2020), and orchards (less than 1%, on average for the period 1970-2020) ( Table 4). The coefficient of vegetation cover, S2, was calculated to range from 0.72 to 0.74, the coefficient of the river basin planning, Xa, from 0.45 to 0.50 Figures 12 and 13.

Current Erosion and the Impact of Land Use on Soil Erosion Intensity
The dominant erosion form in the study area is from surface runoff, including severe forms of erosion such as rills, gullies, and ravines. Surface erosion has affected all soils and slopes, and it is most pronounced on steep slopes with scarce or degraded vegetation cover (Figure 14).  The IntErO model was used to calculate the soil erosion intensity and the peak discharge for the Miocki Potok river basin. The results of the impact of land use changes on soil erosion intensity in the Miocki Potok river basin in six time profiles (years: 1970, 1980, 1990, 2000, 2010, and 2020) are presented in the Tables 5 and 6 and Figures 15 and 16. Calculated soil losses in 57 river basins of the Polimlje region are presented in Table 7, as well. The IntErO model was used to calculate the soil erosion intensity and the peak discharge for the Miocki Potok river basin. The results of the impact of land use changes on soil erosion intensity in the Miocki Potok river basin in six time profiles (years: 1970, 1980, 1990, 2000, 2010, and 2020) are presented in the Tables 5 and 6 and Figures 15 and 16. Calculated soil losses in 57 river basins of the Polimlje region are presented in Table 7, as well.     Table 7, as well.     Source: [69][70][71][72][73]; * Data for the Studied Miocki Potok Basin [13].
Asymmetry of the river basin, a, is calculated at 0.55, and, according to that result, there is a possibility for large flood waves to appear in the Miocki Potok basin. Soil loss, G year , was calculated as 13680 m 3 year −1 and the specific soil loss is 333 m 3 km −2 year −1 (240 t km −2 year −1 ). The value of the Z coefficient of 0.439 indicates that the river basin belongs to destruction category III. The strength of the erosion process was moderate, and, according to the erosion type, was mixed erosion.
According to our analysis, land use changes, specifically the change of the forest structure (fs), as well as the slight decrease of the plow-lands (fg), influenced the decrease of the soil erosion intensity of 14% in the Miocki Potok River Basin.
In parallel, the value of peak discharge, Qmax in the 1970-2020 interval has decreased by 3.5% ( Figure 16).

Lim River Reservoir Sedimentation Survey
From the bathymetry of Potpec reservoir, it can be seen that the backfilling of the reservoir shows a declining trend compared to previous measurement periods ( Table 8). The small amount of sediment in the reservoir can also be affected by the occasional flush of water through the lower dam gates. In order to carry out the IntErO model verification for the study area, the average sediment yields of the 57 river basins of the Polimlje region in Montenegro were calculated as 347,273 m 3 year −1 [13].      The Potpec reservoir undertakes the Polimje region, which includes the Miocki Potok basin. The difference between measured accumulation and predicted sediment yield is less than 1%. This result indicates that the use of the IntErO model is acceptable for the preparation of land use scenarios on erosion intensity and runoff in the Polimlje basins and similar basins in the region; though factors such as the trap efficiency of the Potpec reservoir and the bottom gates flush routing should also be considered for more accurate predictions. The proportion of the incoming sediment that actually accumulates in a reservoir depends on many factors, including the reservoir trap efficiency, reservoir capacity, inflow discharge, and reservoir length.
Brune [74] and Churchill [75] presented curves relating to the percentage of sediment that is likely to be retained in a reservoir to the ratio of the reservoir capacity to the mean annual runoff.
Roberts [76] modified Churchill's curve to determine the trap efficiency of a reservoir with respect to a sedimentation index, in which the reservoir length is an important factor. The Potpec reservoir is very long (about 17 km) and very narrow (width ranging between 100 and 200 m). The longitudinal profiles showed that only in the very vicinity of the dam was there evidence of bottom-gate sediment flushing. On the basis of these considerations, though we do not have data of sediment transport beyond the dam, we can assume that likely 80 to 90% of the incoming sediment was retained in the reservoir. The sedimentation rate of the Potpec reservoir, obtained from bathymetric measurements, was 340.000 m 3 yr −1 , so the actual sediment input can be assessed to be about 390-400,000 m 3 yr −1 . The average sediment input predicted by the IntEro model was slightly smaller (347,273 m 3 yr −1 ), but given the size of the Lim river catchment with respect to that of the individual 57 small catchments used for the model calibration, the model predictions could be considered accurate and reliable, and the IntErO model could be used to quantify the sediment yield variations induced by the land use changes associated with the demographic changes in the Polimje region and the Miocki Potok basin as well.

Discussion
Over the last six decades, anthropogenic factors have significantly increased the pressure on agricultural and forest land in Montenegro, degrading the vegetation cover, which eventually results in serious degradation and loss of fertile soil. In Montenegro, water erosion is the most important erosion type. Water erosion is primarily caused by precipitation and consecutive runoff, but also by fluvial erosion in streams [77]. Given the extreme precipitation values in some parts of the country (the highest of Europe), the influence of this erosion type on the landscape is enormous.
According to Nyssen [67], land-use changes in Montenegro show that, on average, the area covered by dense vegetation increased from about 35% in the early-20th century to 56% of today. In the mountain region, the wooded areas increased slightly. The share of agricultural land, represented by meadow and farmland, remained constant. This pattern is similar to our findings, whereby the forests cover increased from 46% in 1970 to 52% in 2020. Our results for forest land were about four percent less than those of Nyssen [67] because the Miocki Potok study area had a higher than average number of villages and population density.
Depending on available equipment and accessibility of the terrain, different verification methods can be applied. With the idea to compare the data processed by the model measurements of sediment yield on erosion plots were conducted at the Rokava (Dragonja) River basin in Slovenia by Zorn and Komac in 2011 [23]. In North Macedonia Milevski in 2008, [21] found a very good correspondence between the results obtained for the Bregalnica river basin, using the Erosion Potential Method and on-site measurements. Bagherzadeh, in 2011 [78], verified model outputs by a field survey using GPS and a visual comparison of areas characterized as areas with moderate and heavy annual sediment yields. Kouhpeima, in 2011 [33], analyzed five different catchments in Semnan Province in Iran (Amrovan, Atary, Ali Abad, Ebrahim Abad, Royan catchments) and used comparison results to measure sediment deposits in the reservoir as a verification method. Haghizadeh in 2009 [79] used a comparison of model output results with field observations and a GLASGOD (Global Assessment of Soil Degradation) map as a verification method. In the Prescudin catchment in Italy, Bemporad in 1997 [80] measured the sediment deposits, showing minimum deviation between predicted and measured sediment yield values.
Nuclear probes for suspended-load measurements were used at the Esino and Musone river basin in Italy by Tazioli [22]. The measurements exhibited some deviations in comparison with the overall sediment yield production estimated with the EPM method, but, overall, a good order-of-magnitude correspondence was obtained concerning yearly sediment yield. It was concluded that further measurements were necessary because the EPM considers total sediment load, whereas the measurements were conducted to take into account suspended load only.
Other verification methods include the use of a personal digital assistant (PDA) device and on-site observations [18], and certain verification methods remain unspecified in papers [81,82] but provide poor overall ratings for the EPM method, which is said to overestimate sediment yield [81].
Surface erosion has affected all soils and slopes, and it is most pronounced on steep slopes with scarce or degraded vegetation cover. A similar situation has been described for Brazilian watersheds, where higher slopes [83][84][85] and poorly designed rural roads [85] are the main reason for increasing surface erosion. Under the same conditions of soil, slope, and vegetation cover, Avanzi [86] reported that a more concave slope could concentrate the runoff flow leading to a greater erosion rate. These changes can significantly affect the further purpose and use of land, especially in the future agricultural and overall economic development [87].
The calculated erosion rates were rather low at international standards, and that explained the nature of the catchment Table (Tables 6 and 7). Besides, Babic Mladenovic [88] calculated real soil losses for the total Lim River basin at 350 m 3 km −2 year −1 . By using the IntErO model, we found the average value was 331.78 m 3 km −2 year −1 for the 57 river basins of Polimlje in Montenegro [13], and in the range of 394 m 3 km −2 yr −1 (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) to 333 m 3 km −2 yr −1 (2020) for the Miocki Potok river basin ( Figure 15). This correspondence suggests that the results of the assessment obtained by the IntErO model were appropriate for the study area.
In parallel, the value of peak discharge, Qmax in the 1970-2020 interval, decreased by 3.5% ( Figure 16). This is in line with the findings of Nyssen [67] that in Montenegro, in recent times, runoff has decreased due to an increase in vegetation cover.
The calculated erosion rates are relatively low compared to other Mediterranean regions. That is likely related to the large occurrence of forest in the Polimje region and especially in the Miocki Potok basin.
The demographic changes observed in the study area are similar and with similar effects to those recorded in other Mediterranean areas. In the northern Adriatic rivers, for instance, the sediment flux to the sea was much reduced because of the reforestation, both spontaneous and controlled by the forest authorities, following the abandonment of cultivated areas, especially at in the mountain areas for the fast industrialization of Italy and the demand for a larger workforce [89].

Conclusions
The scientific literature about the impact of land use/cover change on the hydrological systems in the Mediterranean is rather abundant. Typically, studies focussed on the influence of different cultivations on runoff generation and soil loss. Others investigated the role of deforestation in increasing the risk of flooding and/or land degradation. Very few papers report about the effect of vegetation regeneration on sediment yield and hydrology at the catchment scale. This study used the IntErO model to quantify the sediment yield and runoff response of a small watershed in Montenegro to past variations of land use and natural forest regeneration, driven by rural population migration. The results of this study, though obtained for a small catchment in a mountainous area of Montenegro, provide an insight on the interconnection between social issues and hydrological processes that can be used to depict realistic scenarios of soil erosion and hydrological changes in other regions with a similar bio-physical setting. The approach proposed in this study proved to be able to indicate with a good level of accuracy the sediment yield and hydrological changes that may be expected if vegetation is allowed to regrow.
In the study area, the migration of rural population to towns proved to be a crucial factor. In fact, the resulting abandonment of farmland in recent years resulted in an increase in the forested area. In the Miocki Potok watershed, whose biophysical characteristics can be considered representative of the whole region of the upper Lim River Basin, the forests increased from 46% in 1970 to 52% in 2020. Such an increase in the forested area resulted in a decrease of 14% and 3.5% for soil loss and peak discharge, respectively, which is in line with the findings of Nyssen et al. (2014) showing that in Montenegro, recent runoff decreases are due to vegetation increase. During the 1970-2020 interval, in fact, no significant change in annual precipitation was found. Yet, no relevant change in both the daily precipitation maximum intensity and the number of days with rainfall higher than 20 mm in a year was observed. These results proved that during the 1970-2020 interval, rainfall had a secondary role compared to the forest regrowth in the predicted hydrological changes. The IntErO model accuracy in predicting sediment yield was tested against the sediment accumulation within the Potpec reservoir on the Lim River. The IntErO model predicted an average annual sediment yield of 330,000 m 3 , whereas from repeated topographic filed surveys, an average annual sedimentation rate in the reservoir of 340,000 m 3 . The Potpec reservoir is very long and narrow, and the average volume of the annually accumulated sediment was assessed by Roberts [76] equation around 400,000 m 3 yr −1 . This result and the IntErO model prediction were very close, thus confirming the usefulness of the IntErO model to calculate runoff and sediment yield at a watershed scale in southeastern Europe, and areas with physical characteristics similar to the Polimlje Region of Montenegro. In the study area, the calculated erosion rates are relatively low compared to other Mediterranean regions. This result can be accounted for by the large proportion of forested areas in the Polimje region and especially in the Miocki Potok basin.
Seven longitudinal profiles of the Potpec reservoir bottom were surveyed in 1967, 1976, 1981, 1999, 2000, 2012, and 2017. The largest sediment accumulation rates occurred from 1967 to 1981, then a period of decreasing sedimentation followed with very low rates observed in the last two decades. This pattern, though with a decade delay, very well reflects the migration phases of the rural population and the consequent abandonment of their farmlands.
The fast industrialization of Montenegro and the growing need for a larger workforce produced deep demographic changes, especially in the rural, mountain areas. The consequent abandonment of the cultivated areas resulted in a natural process of reforestation whose hydrological implications are well predicted by the IntErO model.
Though the approach proposed in this study proved to be successful for the Balkan region, other studies in other countries are necessary to define the environmental constraints of the model application.