Ecological-Economic (Eco-Eco) Modelling in the River Basins of Mountainous Regions: Impact of Land Cover Changes on Sediment Yield in the Velicka Rijeka, Montenegro

This paper presents an Ecological-Economic (Eco-Eco) modelling using the Intensity of Erosion and Outflow (IntErO) model for calculation of sediment yield and runoff assessing the impacts of different land covers on soil erosion intensity. Calculations have been made for the Velicka River basin, which is one of 57 sub-basins of the Lim River in the Northeast Montenegro. Several different land use scenarios were then simulated in the model in order to find the optimal scenario of land use for intensive seed potato production. The results of Ecological (Eco-) analysis shown that the real soil loss under current conditions is 18148 m3yr. If seed potato production is introduced, the model calculated a soil loss of 20834 m3yr as sediment yield. In order to balance the damage caused by the introduction of seed potato production we considered also the ecological measure of afforestation to reduce soil loss caused by seed potato production. The model calculated that afforestation would result in a decrease of sediment yield to 17886 m3yr. The results of Economic (-Eco) analysis revealed that the investment of €3,385 per ha for the establishment of the seed potato production will provide the income for the farmers of €15,000 per hectare annually. In parallel, we proposed the investment for the protection of the area (258 ha) with afforestation that amounts to €330,608 (€1,281 per ha), for the period of two years, with no other costs in the next decade. The research results demonstrate that the application of the Eco-Eco modelling, by using the IntErO model for studying the effect of soil erosion and possible land use for intensive seed potato production in the Velicka River Basin provides cost effective solutions for the benefit of the local population.


Introduction
In mountain areas such as the most part of Montenegro, watersheds are often affected by natural disasters including overflows, floods and inundations, erosion and landslides. Soil erosion is one of the most widespread and a major environmental threat which decreases agricultural productivity and affects water quality (Poesen et al., 1997;Weltzin et al., 2003;Nearing et al., 2005). There are several stages/types of water erosion, including splash, sheet, rill, gully and stream bank erosion (Toy et al., 2002;Poesen et al., 2003;Khaledi Darvishan et al., 2014 and. Their uncertainties regarding the relative significance of each characteristic (Bockstael et al., 1995), sediment yield can be used to show how land cover changes can affect soil erosion and, therefore, the ecological and economic conditions of a watershed. The importance of land cover changes on runoff, soil erosion and the economic conditions of the watershed has been rarely considered as one of the main anthropogenic effects. Land cover management is one of the main issues of sustainable development to design and implement reclamation measures especially in degraded areas (Ballesteros Cánovas et al., 2017).
The objective of this paper is to apply ecologicaleconomic modeling by using the IntErO computer graphic model to predict soil erosion scenarios associated with intensive seed potato production at the Velicka Rijeka River Basin in the Northeast Mountainous region in Montenegro. This study aims also to emphasize how ecological-economic modelling approach may contribute to provide policy makers with scientifically based information for best practice decision making in various fields. Namely, it examines ecological aspects of soil erosion and economic development of areas affected by soil erosion processes and proposes measures to overcome these issues.

Research area
The river basin of the Velicka Rijeka is the right-hand tributary of the river Lim in the Northeast Montenegro. This basin covers a surface area of 32.2 km 2 encompassing the villages of Velika, Volujak and Radevic. It is located 5.3 km north of Plav; 9 km south of Sekular, settlement of Spalevici; 15 km south-east of Andrijevica (Fig. 1).
The catchment slopes are very steep all the way from the Prijedolska glava (Hmax, 2077 m asl) to the confluence of the Velicka Rijeka into the Lim River (H min 879 m asl) for a distance of 5.20 km. The length of the main watercourse is 6.9 km. The shortest distance between the source and the mouth is 5.4 km; length of the basin, measured with a series of parallel lines, Lb, is 8.9 km. Average altitude of the catchment is 1455.83 m asl, and the average height difference is 576.83 m.
Satellite imagery, available from the Google Earth and Google Maps, were used to estimate standard morphometric methods (Spalevic et al., 2014b) and analyze the erosion rills density and the depth of the erosion base, to measure specific lengths such as the natural length of the main watercourse and tributaries of 1 st and 2 nd class, the length of the watershed and other physical-geographical characteristics.
negative effects on soil thickness and fertility, plant cover, runoff coefficient and flood risk may be remarkable, hence, soil erosion and sediment yield studies are of great interest in the world (especially in arid and semi-arid regions where soil and water resources are highly vulnerable. The widespread environmental impacts of soil erosion and loss are often not enough faced by the governments (Behzadfar et al., 2014).
Thus, the accurate understanding and quantification of soil erosion at watershed scale is essential to address many environmental problems influenced by the amount of sediment transported and deposited out of the basins. On the other hand, multi-years measurements of sediment transport at the watershed outlet may represent the soil loss in a watershed (Tazioli, 2009).
Multi-years measurements of sediment load transport at the watershed outlet represents the soil loss in a watershed (Tazioli, 2009) and can be used to calibrate soil erosion models (Tazioli, 2009;Khaledi Darvishan et al., 2010).
Modelling, when calibrated, are useful tools to test hypothesis and to evaluate the amount of discharge and erosion in a watershed, especially when hydrometric data are not available (Behzadfar et al., 2014). Therefore, erosion models have been developed to assess soil erosion and sediment yield, based on some simple empirical equations such as the Universal Soil Loss Equation (USLE), or some modified and updated versions of it Smith, 1965, 1978).
Since it is difficult to accurately measure soil erosion in the field, also the performance assessment of soil erosion models is difficult (Conoscenti et al., 2008, Rawat et al., 2011. By contrast, sediment yield models are easier to apply and to be tested, because the data for these models can be measured at the watershed outlet (Kinnell and Riss 1998;Erskine et al., 2002;Kinnell, 2010).
Among several models, Erosion Potential Method -EPM, originally developed in Yugoslavia by Gavrilovic (1972), was in recent times repeatedly applied to watersheds in the Apennines and in the Balkan Peninsula (Stefanovic, 2004;Zorn and Komac, 2009;Milevski et al., 2008;Tazioli, 2009;Blinkov and Kostadinov, 2010;Ristic et al., 2012;Kostadinov et al., 2006Kostadinov et al., , 2014Spalevic et al., 2014a), but also in other regions in the world, for example in arid and semi-arid areas of the south-western USA (Gavrilovic, 1988), Saudi Arabia (Al-Turki et al., 2015). The method considers the main factors affecting erosion in a catchment, i.e. temperature, mean annual rainfall, soil use, geological properties and some other minor features of the catchment.
The Intensity of Erosion and Outflow model -IntErO program package (Spalevic, 2011), was developed to predict the intensity of soil erosion and the runoff peak discharge in a watershed. It is a computer-graphic method based on the Erosion Potential Method -EPM, which is embedded in its algorithm (Spalevic et al., 2013a).
The functions and processes of ecosystems are not easily characterized because of their complexity interrelation. There are some ecosystem characteristics which can be used to evaluate stress, such as runoff and sediment yield. Some ecosystems characteristics have obvious and immediate economic and/or human significance, while others appear important in a longer term or in more global sense (Bockstael et al., 1995). Though there are some considerable 603 The research part related to geology and soil is based on previous geological (Zivaljevic, 1989) and pedological studies (Fustic and Djuretic, 2000), who analyzed all geological formations and soils of Montenegro. Furthermore, we collected some soil samples for chemical and physical analysis. The grain size composition of the soil was determined by the pipette method. The soil samples were air-dried at 105 °C sifted through 2 mm sieve and dispersed using sodium pyrophosphate. Total carbonates were determined by the volumetric Scheibler method; the soil reaction (pH in H2O and nKCl) was determined with a potentiometer; the content of the total organic matter was determined by the Kotzman method; easily accessible phosphorous and potassium were determined by the Almethod and the adsorptive complex (y1, S, T, V) was determined by the Kappen method (Spalevic et al., 2013b). We used IntErO model (www.agricultforest.ac. me/Spalevic/IntErO), based on Erosion Potential Method -EPM (Gavrilovic, 1972) and designed to assess annual erosion rates: km -2 ]; lp -Length of the principal waterway [km]; la -Cumulated length of secondary waterways [km]; L -Cumulated length of the principal and secondary waterways [km]; Gy -Actual sediment yield [m 3 yr -1 ].
EPM is widely used in the Balkan Region because of its relative simplicity and it is preferred as a local model for calculation of soil erosion intensity for the territory of Ex-Yugoslavia (Spalevic et al., 2013c). The use of EPM, including the River Basin model, has been used widely in Montenegro, especially in the Region of Polimlje (Spalevic et al., 2014b), representing a standardized approach.
In comparison with some other procedures, the EPM/IntErO model does not explore the physics of erosion processes; therefore, it is suitable for areas where basic data are available, or where there is a lack of previous erosion research. As such, the model can provide not only the amount of sediment yield, but also the erosion intensity as a preliminary result and indications or areas of potential erosion threats (Dragicevic et al., 2016).

Climatic characteristics
The effects of climate have the high impacts on land degradation; rainfalls and torrential rains are amongst the main triggers of soil erosion processes. Global warming is leading to a more vigorous hydrological cycle, including higher total precipitation and more frequent high intensity rainfall events. Rainfall amounts and intensities increased around the world during the 20 th century, and according to climate change models they are expected to continue to increase during the 21 st century (Nearing et al., 2004). These rainfall changes, along with expected changes in temperature have significant impacts on soil erosion rates.
Montenegro is experiencing increasing temperatures and evapotranspiration, most notably in the northern mountainous region. The 2001-2010 decade was the warmest since records began, with the most prominent changes in the northern mountainous region of +1.4 °C and with a decrease in the number of frost days and very cold days and nights. Changing rainfall pattern is also forecasted in the near future (more precipitations in winter, less in summer) increasing erosion, flood risk (winter) and water stress (summer). The analysis of the climatic patterns undertaken confirmed that climate in Montenegro has already changed and that the main impacts foreseen for temperatures and extreme events are confirmed (Froslev Christensen and Spalevic, 2017).
Regarding rainfall there has been no significant reduction in the average annual precipitation: rainfall has increased in autumn whereas it has decreased during spring, summer and winter. However, there has been a damaging and significant increase in the number of extreme weather events. This pertains especially to heat waves, that are increasingly frequent, and their length shows a high year-toyear variability. Secondly, but equally important, storms have become more frequent and more intense since 1998, resulting in huge amounts of precipitation and high flooding.
Extreme whether events (e.g. droughts, flooding and heatwaves) are increasingly impacting natural resources (soils, water bodies, pastures, others). Moreover, heavy snowfalls and flash floods are becoming more common.  (Spalevic et al., 2013a). The Velicka Rijeka catchment has a continental climate, with rainy autumns and springs and cold winters. The absolute air temperatures recorded range from a minimum of -29.8 o C up to 35 °C. On the basis of the available data, the average annual air temperature, t0, is 8.1°C; average annual precipitation is 1182.3 mm. Calculated temperature coefficient for this area, T, is 0.95; the torrential rain, hb, is calculated to be 89.4 mm.

Geological structure
Mountains in Montenegro are part of the Dinaric Alps. The terrain around the study area consists of various types of sedimentary, magmatic and metamorphic rocks, ranging in age from Palaeozoic to Quaternary. Most of the study area is underlain by Mesozoic formations of carbonate composition, while magmatic and clastic silicate rocks are significantly less present. Using the Geological map of Montenegro (Zivaljevic, 1989), permeability of the rocks of the study river basin has been defined. The coefficient of the region permeability, S1, is calculated to be 0.9. Rocks of poor permeability (class fo) cover 79% of the study river basin area, very permeable rocks (class fp) cover 13% of the territory and semi-permeable rocks (class fpp) only 8%.

Soils of the area
Prevailing soils in Montenegro are characterized by limited to low fertility (90% of soils), acid reaction (95% of soils in Montenegro are naturally acidic), often skeletal and shallow, with small retention capacity for moisture and nutrients. On the basis of both previous pedological results (Fustic and Djuretic, 2000) and original field and laboratory research, the most important types of soil in the basin are listed according to the percentage distribution: Dystric Cambisols (83.68%), Eutric Cambisols (12.60%), Kalkomelanosols (3.01%), and traces of Fluvisols and Colluvial Fluvisols (0.71%).

Vegetation
The study area belongs to the Dinaridi Province of the Middle-Southern-East European mountainous biogeographical region (Dees et al., 2013). Some wooded areas are found at higher altitudes in the Velicka river basin and on the top of Cakor mountain and its slopes, where forests of spruce (Picea abies), one of the most important tree species from Europe, are found. At Djevojacki krs, towards to the Sabova glava, the basin is covered by forests of Macedonian pine (Pinus peuce), and on the steep slopes towards the village of Velika, the forests of fir (Abies alba) and spruce (Picea abies) are prevailing. To the right side of the basin, around the rural community of Radevici, the beech forests (Fagetum montanum) have been gradually replaced by pastures. In the lower parts of the study river basin, at the confluence of the Velicka Rijeka into the Lim River, mixed forests of beech (Fagetum montanum) and oak (Quercus cerris) are present. In the lower parts, near the river bed, we recorded some hydrophilic forest (Alnetea glutinosae, Salicetea herbacea) and, subordinately, some Betula verrucosa, Quercus cerris, Quercus petraea and Prunus avium.

Model parameters calculation
The coefficient of vegetation cover, S2, was calculated to be 0.73; the coefficient of the river basin planning, Xa, 0.42. According to the analysis of available data, meadows, pastures and orchards cover around 59% of the study river basin. Mountain pastures are the most widespread vegetation cover (44.82%). Meadows cover 12.95%, orchards 1%. Forests are covering 39%: well-constituted forests (30.18%) and degraded forests (9.06%). Arable and cultivated land less than 2%.

Results and Discussion
In order to carry out model verification for the study area, sediment yields were calculated for all the tributaries of the Lim river basins, which include the Velicka river basin. The model results were then compared with the measurements obtained at the Potpec reservoir. Using the Model the sediment yield was calculated to be 347 273 m 3 year -1 ; while actual geodetically performed measurements were 350 000 m 3 year -1 . This validates calculations of the results for sediment yield obtained by the model. This leads to a conclusion that the model is applicable for the observed area (Spalevic et al., 2016).
In Italy, using the same methodology, Tazioli (2009) found that this model corresponds well concerning annual sediment yield using nuclear probes for suspended-load measurements on Musone and Esino watersheds. Similar studies were applied earlier at the Prescudin catchment in Italy, (Bemporad et al., 1997) recording a minimum deviation between predicted and measured sediment yield values. At the Bregalnica basin in Macedonia (Milevski et al., 2008), a very good match has been achieved between the results obtained using the EPM method and onsite measurements. It should be highlighted that the EPM/IntErO model considers the total sediment load, whereas most of the measurements conducted in the studies cited take into account suspended load only.

Establishment of seed potato production and afforestation as a measure of soil conservation
We have introduced eco-eco modeling in this research with the main idea to guide the farmers on how to get economic benefits, while respecting at the same time sustainable river basin management. Economic part of the main research hypothesis is that commercial production of seed potato will increase income of farmers, income of the state (by increase of collected taxes) and improve the state balance of payments by reducing the import of food (Montenegro produces approximately 18% of food that is consumed in the country). The increase of income will improve living standard and quality of life in the study catchment, thus reducing migration from this area.
The study river basin is suitable for growing potatoes. Most farmers produce vegetables for home consumption, because of higher value crop opportunities which exist due to the particular comparative advantage of this area. Farmland at +800 m above sea level is ideal for producing high quality seed potato given the clean, virus free, low disease conditions. In an environment where no fruit can be cultivated above 1,300 m asl, highland varieties of potato are an important resource for certain households in the study area -on a commercial basis. For the establishment of intensive agricultural production it is important to improve land preparation, together with fertilizers and manure usage that will increase yields by one-third.
At the initial phase we proposed an expansion of the potato cultivated area that brings increased income and also incurs various production costs as well as cost of nature conservation of the study basin. For the purposes of this research we have taken into account theoretically maximum possible area of 2.58 km 2 (258 ha) for establishment of an intensive agricultural production at this territory. This area extends on a 5 km long course within the territory of village Velika, the valley across the Papratista and Lijeva Rijeka to the village Radevic. In real terms it would be hardly feasible for the agricultural production to be implemented on the entire territory mentioned above, due to property issues, labor shortages and other issues. We decided to prove our concept using simple calculation in accordance with potentials of the nature in the observed area.
However, establishing intensive potato production causes environmental damage by increasing the soil erosion intensity, which can be balanced by the introduction of conservation pathways. In order to calculate these economic benefits, we used the IntErO model to calculate soil erosion and its repercussion on the potato farming economic system. Calculation of seed potato production costs per 1 hectare for the region of Velicka river basin is presented in the Table 1. Potting of potatoes 20 2 Loading -unloading fertilizers 10 3 Spreading of fertilizers 10 4 Loading -unloading of potatoes 20 5 Seeding of potato 40 7 Weeding weeds 20 8 Treatment against diseases and pests 20 9 Extraction of potatoes 400 10 Loading -unloading of potatoes 100 Total 640 In order to calculate costs and income of seed potato production we have quantified the following costs: production operations costs, raw material costs and labor costs. In agreement with Article 24 of Montenegrin VAT Law, VAT rate for seed potato is 7%; however, we did not examine in detail taxation issues in order to not overcomplicate the model.
Taking into account that the study area is 258 ha and the calculated profit per hectare (see table 1d) is €10,629, this adds to the total profit of €2,742,282 annually. This amounts to a total nominal profit of €27,422,820 for the next decade (without taking into account price changes and time value of money), which is the basic time frame for calculation of afforestation costs in this research paper ( Table 2).
The model has taken into account all of the above parameters (27 input data and 22 results) ( Table 4) and it was calculated that the real soil loss under current conditions is 18,148 m³yr -1 . If seed potato production is introduced, the model calculated a soil loss of 20,834 m³yr -1 as sediment yield. In order to balance the damage caused by the introduction of seed potato production we considered also the ecological measure of afforestation to reduce soil loss caused by seed potato production. The model calculated that afforestation would result in a decrease of sediment yield to 17,886 m³yr -1 . The outcomes of model calculation are presented in Fig. 2.
From these results the twofold benefit of introducing the potato production coupled with afforestation is evident. On the one side, we have an increase in the farmer income and on the other side with planting the new areas under the forest a significant environmental improvement is achieved.  (1), under the seed potato production (2), under the seed potato production after the afforestation (3)

Conclusions
Several different land use scenarios were simulated using the IntErO model in order to find the optimal scenario of land use for intensive seed potato production in the Velicka Rijeka in Montenegro. The main points of research are the interactions between seed potato production, as an intense agricultural/economic activity, and afforestation, as a wellknown soil and water conservation activity. It can be simply concluded that agriculture can be used more conservatively with afforestation and, as a kind of agro-forestry, has the ability to decrease even soil loss in the study area. The effects of afforestation on various variables such as soil physical, chemical and biological characteristics, water infiltration, runoff and soil loss and etc. have been already proven in previous studies but the interaction effect between afforestation and seed potato production especially in soil loss was studied and proved in present research. The abilities of IntErO to simply investigate the effect of various scenarios of land use change made lots of such studies possible without doing any changes in the actual land uses. This advantage can be a very important tool on the working table of managers and policy planners. The quantitative results of soil loss in studied scenarios also showed that 608 afforestation can compensate the increasing effect of seed potato production on soil loss. In other word, seed potato production increased soil loss by 15%, while seed potato production after afforestation decreased it again by even more than the same rate.
There is another important point about the type of agricultural activities as well as the species of trees used for afforestation in each study area and should be taken into account. Technical notes of planting the selected tree species such as the distance and network dimensions and similar can play a very important role in the final ecologicaleconomic (Eco-Eco) modelling and results.