Effects of Erosion Control Works: Case Study–Reservoir Celije, Rasina River Basin, the Zapadna Morava River (Serbia)

Abstract

The aim of this research was to analyze the impact of implemented erosion control works (ECW) on soil erosion intensity in the watershed of the Ćelije reservoir (Rasina River) in the period between 1968 and 2022. The Erosion Potential Method was used to calculate the annual gross erosion (W), sediment transport (G), and erosion coefficient (Z) in the study area. As a result of the performed ECW there was a general decreasing trend in the intensity of soil erosion processes in the last 54 years. The specific annual gross erosion was 1189.12 m³/km⁻²/year⁻¹ in 1968, while in 2022 it was 554.20 m³/km⁻²/year⁻¹. The specific sediment transport was 540.18 m³/km⁻²/year⁻¹ in 1968 and 253.55 m³/km⁻²/year⁻¹ in 2022. Due to the changes in the intensity of erosion processes, the specific annual gross erosion decreased by 634.92 m³/km⁻²/year⁻¹ and the specific sediment transport decreased by 286.63 m³/km⁻²/year⁻¹. The erosion coefficient was reduced from Z = 0.62 to Z = 0.35. A dependence between the slope of siltation and the natural bed slope was defined. The results show a significant correlation between erosion intensity and performed ECW, providing a basis for future watershed management and defining a strategy for soil erosion control in the Ćelije reservoir watershed.

1. Introduction

Soil erosion is the most widespread form of land degradation. Due to climate change, the intensity of erosion processes is expected to increase. Compared to other forms of erosion, water erosion causes the greatest losses of soil. The rate of soil erosion is primarily affected by topography, precipitation, soil type, land use and land management practices [1,2,3,4].

Disturbances in watersheds due to deforestation, mining activities, construction of facilities and other infrastructure projects lead to dramatic erosion processes, increased soil loss and sediment transport into watercourses [5]. Sediment yield in river basins and sediment transport in watercourses represent two components of a global natural process that, due to its consequences, have significant ecological and hydrological importance.

The territory of Serbia is threatened by soil erosion of various intensities, with 36% of the territory affected by erosion processes categorized as excessive, intensive or medium intensity [6].

The Ćelije reservoir was formed in 1979 by constructing a 55 m high dam in the Zlatarska Gorge (Rasina, right tributary of the Zapadna Morava). The initial purpose of the reservoir was flood protection, sediment retention, and irrigation, but later it began to be used for water supply. The Ćelije reservoir is a strategic part of the Rasina-Pomoravlje regional water supply system, which includes the downstream part of the Zapadna Morava River and the upstream part of the Velika Morava River. According to measurements from 1926 to 1975, the average annual flow of Rasina river was 5.9 m³/s, so the water retention time in the reservoir in 1979 was about 80 days. Due to global warming, the flow through the reservoir has tended to decrease. In the period between 2001 and 2015, the average annual flow was 4.5 m³/s, changing the retention time of the reservoir to about 105 days. Due to the influence of the reservoir within the watershed, there are no major issues with flash floods. The management of Ćelije reservoir is overseen by the Public Enterprise Srbijavode, while the dam is under the jurisdiction of the Public Enterprise Vodovod Kruševac. There has not been any reclamation work conducted, such as dredging, resulting in a reduction of the lake’s length by 1 km due to sediment accumulation at the entrance to the reservoir. The projected functioning of the Ćelije reservoir was until 2030, but this has been extended by reducing sediment production in the watershed due to implementation of erosion control measures.

The successful functioning of the reservoir over time has been significantly endangered by the reduction of the storage capacity due to sediment siltation [7]. In addition to reducing the storage capacity and water volume in the reservoir, erosion processes contribute to chemical pollution of water with materials of both natural and anthropogenic origin. Chemical substances of natural origin, resulting from the decay of rocks, are transported into the hydrographic network through sediment erosion. On the other hand, anthropogenic chemical pollution occurs due to the runoff of artificial fertilizers and pesticides from agricultural areas into river basins [8]. Due to erosion processes, many reservoirs in Serbia have been rapidly silted up with sediment and practically converted into sediment reservoirs. A typical example is the Zvornik reservoir on the Drina River, as well as reservoirs on the Zapadna Morava River, and so on [9].

Soil erosion is the main cause of reservoir sedimentation, meaning that sedimentation is also a complex process that depends on numerous natural and anthropogenic factors [10,11]. Several research studies have been conducted on erosion processes, sediment yield and transport in reservoir basins, indicating the significance of sediment issues [12,13,14,15,16,17,18,19].

Erosion and torrent control works in Serbia started at the end of 19th century, but not until 1907 as organized activity. Since then, a significant amount of erosion control works have been performed in the territory of Serbia [20].

Properly designed and performed erosion control works cause the decrease of soil erosion intensity. In this respect, technical and biological works have different activities and effects. Biological works such as afforestation, reseeding, and orchard establishment on terraces are carried out on watershed slopes. Technical works can be longitudinal, for the protection of stream banks against erosion and sloughing and transversal for prevention of torrent beds erosion and retaining the sediment (primarily bedload) [21,22].

Multiple studies on the effects of erosion control works at the national, regional [23,24,25,26,27,28] and global levels indicate their significance in reducing soil erosion intensity [29,30,31,32,33]. Technical works provide the desired results in a shorter time period, while biological and biotechnical works provide a long-lasting solution to erosion and sedimentation issues within the watershed.

The effects of check dams are significant and studying them enables improvement and rationalization of erosion control works. By constructing check dams in their accumulation area, sediment is retained, forming a siltation behind the structure. The slope of the upper surface of siltation represents the slope of siltation, and acts as a practical representation of specific cases of the equilibrium bed slope [34,35,36,37,38,39]. Field research on the effects of check dams in Serbia shows that the best results in calculating the slope of siltation are based on the regional analytical dependence of the slope of siltation on natural bed slope.

The main aim of this research is to determine the spatial and temporal changes in soil erosion intensity and determine the effects of erosion control works on sediment yield, sediment transport and spatial distribution of soil erosion processes in the watershed of the Ćelije reservoir based on field research conducted across two time periods (1968 and 2022). The analysis of the effects of erosion control works is performed by comparing soil erosion intensity, sediment yield and transport before and after the works, as well as by calculating the amount of sediment retained upstream of the check dams and defining the slope of siltation.

2. Materials and Methods

2.1. Study Area

The Rasina River is a right tributary of the Zapadna Morava River, with a total watershed area of 981 km². The Rasina springs on the slopes of the Goč and Željin mountains at an elevation of 1340 m above sea level, formed by the merging of the Velika and Burmanska rivers (Figure 1). The lowest elevation of the Rasina is 134 m above sea level, at the confluence in Zapadna Morava, five kilometers downstream from Kruševac. Up to the village of Razbojna, it flows through a canyon-like valley with small erosion enlargements, such as the one near Brus. Between Razbojna and Brus, it flows through a shallow and wide channel through the Dobroljubac basin before entering the Zlatarska Gorge where it forms several entrenched meanders. The watershed is elongated with a developed hydrographic network. Apart from the Zagrža river, which is the left tributary, all other tributaries are from the right. All these rivers flow through canyon-like valleys with steep slopes transporting sediment. The largest among them are the Graševac River, coming from Kopaonik, and the Blatašnica River, from Jastrebac.

The Rasina River drains the slopes of the Jastrebac and Kopaonik mountains, reaching elevations of 1500 m above sea level and 1900 m above sea level, respectively. According to the topographic map, the watershed area upstream of the dam is 609.15 km². The perimeter of the watershed is 184.47 km. The river length is 69.95 km, and the length of the watershed is 49.39 km. The lowest elevation in the watershed is 239 m above sea level at the dam profile, while the highest elevation is 1936 m above sea level. The average elevation in the watershed is 695 m above sea level, with a mean elevation difference of 447 m. The average slope of the watershed is 34.65%. Physical characteristics of the Celije Reservoir watershed were derived from a topographic map, scale 1:25,000 and digital elevation model (SRTM) (

Table 1. Physical characteristics of the Celije Reservoir watershed.

Parameter	Mark	Value	Unit
Watershed area	A	609.15	km2
Perimeter	O	184.47	km
River length	L	69.95	km
Watershed length	Ls	49.39	km
Peak point of the watershed	Hmax	1936	m
Confluence point of the watershed	Hmin	239	m
Mean slope of the riverbed	Iu	1.62	%
Mean slope of terrain	Isr	34.65	%
Mean altitude	Hsr	695	m
Mean altitude difference	D	447	m

).

Slope angles up to 10 degrees cover 29% of the watershed area. Dominant slopes are between 10–25 degrees, covering 60% of the surface. Areas with slopes greater than 25 degrees constitute 11% of the surface (Figure 2a). Analysis of hypsometric characteristics of the Ćelije reservoir watershed area shows that 23% of the surface is located between 300 and 500 m above sea level, 58% of the surface is between 500 and 1000 m above sea level and 18% is located above 1000 m. Only 1% is below 300 m (Figure 2b). A slope map and a hypsometric map were produced using a topographic map at a scale of 1:25,000 and digital elevation model (SRTM).

These natural characteristics indicate a significant predisposition of the watershed area for the development of erosion processes, sediment yield and transport.

The geological characteristics are reviewed on the basis of the Main Geological Map of SFR Yugoslavia, scale 1:100,000. The geological composition of the Ćelije reservoir watershed is predominantly characterized by the following geological formations: Paleozoic bedrock complex (crystalline schists), Mesozoic bedrock complex (flysch), Paleogene bedrock complex, Neogene bedrock complex and Quaternary deposits. Legend of the geological map is reduced, rocks presented on area smaller than 3 km² are not shown (Figure 3).

According to the resistance to erosion processes, bedrocks are categorized into four categories: very solid rocks, conditionally solid rocks, conditionally erodible rocks and very erodible rocks. Conditionally solid are present on 13.75% of the area, conditionally erodible rocks on 38.27% and very erodible rocks on 44.88%, meaning that more than 90% of the area is potentially threatened by erosion. A map of the rock resistance to erosion is derived from a geological map of the study area (Figure 4b).

The dominant soil types in the watershed are of weaker structure. The most prevalent soil types are eutric cambisols (35%), followed by distric cambisols (12.8%), ranker eutric soils (8.7%), ranker distric soils (6.6%) and regosol (7.5%), while other soils are present on smaller areas. Soil characteristics are reviewed on the basis of the Main Soil Map of SFR Yugoslavia, scale 1:50,000 (Figure 4a).

The lowland part of the Rasina River watershed is characterized by a moderately continental climate, while the hilly-mountainous area of the watershed has uneven elevation influenced by the continental climate.

To determining the climatic characteristics of the Ćelije reservoir watershed area for the period from 1990 to 2022 (

Table 2. Meteorological stations in the Celije reservoir watershed.

Meteorological Stations	Coordinates		Elevation	Average Annual Temperature	Average Annual Rainfall
	X	Y
			m a.s.l.	°C	mm
Kruševac	7,528,274.36	4,824,754.19	166	12.4	738
Brus	7,502,705.92	4,804,323.12	426	10.6	822
Blace	7,524,339.64	4,795,099.41	425	10.6	651
Goč	7,487,890.04	4,822,849.21	990	8.3	978
Kopaonik	7,484,294.00	4,793,337.00	1710	4.2	1035
Jastrebac	7,528,303.00	4,814,539.00	256	11.6	884

), data from meteorological stations in Kruševac, Brus, Blace, Goč, Kopaonik and Jastrebac were used. Isohyet and isoterme maps for the Study area were developed based on these data [40].

The annual precipitation during the period from 1946 to 1990 at the Kruševac meteorological station was 649 mm, while at the Brus meteorological station it was 678 mm, indicating that the precipitation amount was higher by 89 mm and 144 mm, respectively (Figure 5a).

Analyzing the data from the “Kruševac” meteorological station, the average annual temperature for the lower part of the watershed is 12.4 °C. For the mountainous part of the watershed, the average annual temperature is 4.2 °C, based on data from the “Kopaonik” meteorological station at 1710 m above sea level (Figure 5b).

Compared to the values for the period 1946–1990, the average annual temperature at the Kruševac meteorological station has increased by 1.4 °C.

The Rasina River exhibits all the characteristics of a torrential watercourse because its flow amplitudes are very pronounced. In April, it has 5.5 times more water than in August. The highest flow occurred in April 1958 and amounted to 342 m³/s.

According to the population census data from 2022 (Figure 6), the population in the area of the Ćelije reservoir watershed is 20,402 inhabitants, continuing the trend of population decline since 1953 when there were 34,631 inhabitants [41].

2.2. Methodology

There are several erosion risk assessment models in use today [42,43,44]. The selection of a specific model largely depends on the purpose of the research, data availability, time and finances required for the project. Most models used in soil erosion studies are empirical models. These models vary in complexity and incorporate different techniques and approaches to define erosion intensity, sediment yield and transport [45,46,47,48].

The Erosion Potential Method (EPM), also known as Gavrilović, was developed and calibrated using field research in the Juzna Morava River basin [49,50,51,52,53,54]. This method is widely used in Serbia and the former Yugoslav republics [55,56,57,58,59], as well as in the Mediterranean region, Central Europe, the Middle East, and other climates [60,61,62,63,64,65,66].

The method has provided reliable results for evaluating of intensity of soil erosion processes, mean annual soil erosion rate and sediment transport and implementing erosion control measures.

The biggest challenge lies in verifying the obtained data. At the regional and broader levels, reliable data for comparing estimated and actual soil loss practically does not exist. The method involves both qualitative and quantitative analyses and is adapted to GIS environments. The erosion mapping procedure includes processing various datasets and calculating numerical indicators.

The basic premise of the method is to reduce subjective errors in coefficient estimation to an acceptable level, and the development and improvement of the method are aimed at eliminating them.

The EPM involves classifying the intensity of erosion into five categories, from very slight (V category) to excessive (I category). The primary parameter used to classify the intensity and category of erosion in this method is the erosion coefficient (Z).

Parameters for determining the erosion coefficient include topographic, geological and pedological characteristics, land use practices within the watershed and the degree of erosion vulnerability of the watershed. Each of these parameters can be represented by a digital map and by overlaying these maps a separate composite map is generated, where areas with identical coefficient values are formed. For each of these areas, the average slope of the terrain is determined, and all these data represent input values for calculating the erosion coefficient (Z) according to the following Equation (1):

$$\mathrm{Z}=\mathrm{X}\times \mathrm{a}\times \mathrm{Y}\times (\mathsf{\phi }+\sqrt{\mathrm{I}})$$

(1)

where X represents soil protection coefficient, Y is the coefficient of soil resistance, φ is the erosion and stream network development coefficient and I is the average slope of the watershed.

The values of coefficients used in the EPM are shown in

Table 3. Basic characteristics of X, Y and φ coefficients used in EPM.

Additional Info	X Description—Soil Protection	X
Before implementation of anti-erosion measures	Completely bare, uncultivable land (bare land)	1.00
	Arable land with plowing up or down hill	1.00
	Orchards and vineyards, without ground vegetation	0.70
	Mountain pastures and drylands	0.60
	Meadows, fields and similar agricultural crops	0.40
	Degraded forests and thickets with eroded soil	0.60
	Forests and thickets with good structure and vegetation	0.05
After implementation of anti-erosion measures	Plows with contour plowing (isohypsis direction)	0.60
	Arable land well cared for and protected by mulching	0.50
	Contour strip cultivation with crop rotation (fields)	0.45
	Contour orchards and vineyards	0.30
	Terracing of arable land, terraces and tiers	0.35
	Weeding of bare land and melioration of pastures and drylands	0.30
	Construction of contour trenches of medium density
	Retardation waterways and micro-accumulations	0.25
	Afforestation with construction of tiers	0.10
	Channel regulation, dam construction and channelization	0.70
Y Description—Soil Erodibility		Y
Sand, gravel and loose soils		2.00
Loess, tuffs, salt marshes, steppe soils		1.60
Disintegrated limestones and marls		1.20
Serpentines, red sandstones, flysch deposits		1.10
Podzol soils and alike; decomposed shales: micaschist, gneiss slates, clay slates		1.00
Core limestones and shales, red rocks and humus-silicate soils		0.90
Cambisol and mountain lands		0.80
Vertisol, humogley and wetlands		0.60
Chernozem and alluvial soils of good structure		0.50
Bare compact eruptives (volcanic origin)		0.25
φ Description—Type and Extent of Erosion and Slumps		φ
Watershed completely under gully erosion and primordial processes		1.00
About 80% of the watershed is under furrow and gully erosion		0.90
About 50% of the watershed is under furrow and gully erosion		0.80
The entire watershed is subject to surface erosion		0.70
The entire watershed is under surface erosion, without deep processes		0.60
Land with 50% of the area covered by surface erosion		0.50
Land with 20% of the area covered by surface erosion		0.30
No visible signs of erosion, minor slips and slides in watercourses		0.20

.

Precipitation and temperature data were obtained from six meteorological stations for the period 1990 to 2022. The coefficient Y was determined based on a geological map at a scale of 1:100,000. Coefficients X and ϕ were determined using orthophotos, satellite imagery and field research of erosion conditions during 2022, based on which a land use map was created. Different land use categories were identified and harmonized for both periods.

Morphometric characteristics and slopes were determined using a topographic map at a scale of 1:25,000 and a digital elevation model (SRTM).

The calculated values for each delineated area form the basis for the classification of erosion processes. The coefficient values (Z) are categorized into five classes, which provide an optimal number for graphical representation of the erosion status and surface prevalence (

Table 4. Classification of erosion coefficient Z and specific annual gross erosion W.

Erosion Category	Erosion Intensity	Range of Z	Range of W (m3/km−2/year−1)
I	Excessive Erosion	1.01–1.50	>3000
II	Intensive Erosion	0.71–1.00	1200–3000
III	Medium Erosion	0.41–0.70	800–1200
IV	Weak Erosion	0.21–0.40	400–800
V	Very Weak Erosion	0.01–0.020	100–400

). For practical applications and calculations, each category has been assigned a qualitative name and corresponding mean value for the erosion coefficient (Z). This erosion map thus created serves as the foundation for further application in other areas.

The described procedure resulted in the creation of an erosion map for the entire analyzed area. Data on the surface prevalence of erosion, defined by the mentioned method, serve as the initial basis for calculating sediment yield and transport from the watershed area. These data have enabled the development of new methods for calculating sediment transport.

The calculation of specific annual gross erosion was performed according to the following Equation (2)

$${\mathrm{W}}_{\mathrm{s}\mathrm{p}}=\mathrm{T}\times \mathrm{H}\times \mathsf{\pi }\times \sqrt{{\mathrm{Z}}^3}$$

(2)

where, W_sp is specific annual erosion (m³/km⁻²/year⁻¹), T is the temperature coefficient, H is mean annual precipitation (mm), Z is the erosion coefficient and F is the basin area (km²).

Temperature coefficient (T) is calculated by Equation (3):

$$\mathrm{T}=\sqrt{\frac{\mathrm{t}}{10}+0.1}$$

(3)

where t is the mean annual temperature (°C).

The annual gross erosion for the entire watershed area (W_year) is calculated by multiplying the watershed area (F) by the specific annual gross erosion (W_sp) by Equation (4).

$${\mathrm{W}}_{\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}}=\mathrm{F}\times {\mathrm{W}}_{\mathrm{s}\mathrm{p}} {\mathrm{m}}^3/\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}$$

(4)

Having calculated the total annual soil erosion rates, we calculated the sediment delivery ratio, which is needed for the actual sediment yield calculation. Gavrilovic [52] suggested the following Equation (5) for the determination of the sediment delivery ratio Ru:

$${\mathrm{R}}_{\mathrm{U}}=\frac{{(\mathrm{O}\times \mathrm{D})}^{0.5}}{0.25\times (\mathrm{L}+10.0)}$$

(5)

where O is perimeter of watershed (km), D is average distance of elevation of the watershed (km) and L is length of the watershed (km).

Specific sediment transport is calculated as:

$$\mathrm{G}=\mathrm{W}\times {\mathrm{R}}_{\mathrm{U}} {\mathrm{m}}^3/{\mathrm{km}}^{-2}/{\mathrm{year}}^{-1}$$

(6)

The actual sediment transport is:

$${\mathrm{G}}_{\mathrm{G}}={\mathrm{W}}_{\mathrm{G}}\times {\mathrm{R}}_{\mathrm{U}} {\mathrm{m}}^3/\mathrm{year}$$

(7)

In this study, erosion mapping was based on the Erosion Potential Method successfully integrated into the GIS environment, using the open-source geographic information system QGIS [67]. To define the change in soil erosion intensity, GIS was used to overlay the erosion map from 1968 with the current erosion map from 2022. The modeling, digitalization and mapping were conducted using QGIS 3.32.0.

3. Results

3.1. Analisys of Performed Erosion Control Works (ECW) in the Study Area

3.1.1. Transverse Structures in Torrential Streams (Check Dams)

In the Rasina River watershed, technical erosion control works were actively undertaken to retain the bedload and stabilize torrential channels in the immediate vicinity of the reservoir as well as upstream of the reservoir. These works accompanied a period of dam construction and continued even after the dam was built and reservoir filled.

Between 1972 and 1990, over thirty check dams were constructed (Figure 7), primarily using concrete or stone in cement mortar. A significant number of these structures are still evident in the field. For the purposes of this research, 14 check dams built within the channels of 6 streams were analyzed, representing the entire basin by their characteristics. Ten of these barriers were made of stone in cement mortar (Blatašnica, Žunjska River, Zagrža and Rasina), while four were constructed from concrete (Koznička and Popovačka Rivers).

The check dams constructed within the channels of torrential streams serve multiple purposes: to retain sediment and to decrease bed slope. Reducing the bad slope of the watercourse also implies a reduction in the transport capacity for sediment (

Table 5. The Main Check Dams effects.

	River Basin	No	It	Iz	Iz/It	hk	m	W
			%	%		m	m	m3
1	Blatasnica River	1	1.6	0.28	0.18	4.0	19.5	11,818.18
2		2	1.6	0.24	0.15	3.0	25.0	8272.06
3		3	1.1	0.27	0.25	2.5	28.0	10,542.17
4	Popovacka River	1	3.45	1.43	0.41	1.5	14.5	807.55
5	Zunjska River	1	2.77	1.95	0.70	2.0	47.0	11,463.41
6		2	3.07	1.68	0.55	3.0	35.0	11,330.94
7		3	3.07	1.45	0.47	3.0	22.0	6111.11
8		4	4.22	1.72	0.41	5.0	44.0	22,000.00
9	Koznicka River	1	5.79	4.19	0.72	4.0	48.0	24,000.00
10		2	7.69	4.40	0.57	5.5	19.0	8734.80
11		3	7.42	5.45	0.73	3.0	22.0	5025.38
12	Zagrza	1	2.35	1.02	0.43	3.0	14.0	4736.84
13		2	2.50	1.15	0.46	3.5	15.0	6805.56
14	Rasina River	1	3.85	1.98	0.51	2.5	21.5	3592.91

).

The amount of sediment retained upstream of the check dams depends on the effective height of check dams, the valley geometry upstream of check dams and the slope of siltation. The quantity of the retained sediment is calculated using the Kitin 1975 formula [68]:

$$\mathrm{W}=\frac{1}{2}\mathrm{A}\times \mathrm{L}\Rightarrow \mathrm{W}=\frac{\mathrm{m}\times {\mathrm{h}}_{\mathrm{k}}^2}{2\times ({\mathrm{I}}_{\mathrm{t}}-{\mathrm{I}}_{\mathrm{z}})}$$

(8)

where W is the volume of the retained sediment in the silting-up area (m³), m is mean width of the silting-up area (m), h_k—dam effective height (m), I_t is natural torrent bed slope in decimal form and I_z is slope of siltation in decimal form.

The input parameters in the formula were measured for each check dam directly in the field.

The calculation results indicate that check dam 1 on the Koznicka River retained the highest amount of sediment, which is logical considering the geometric characteristics of the channel upstream of check dam. In comparison with this check dam, check dam 2, which has a greater effective height, retained much less sediment, primarily due to unfavorable geometric characteristics (large bed slope of the channel upstream of the check dam and narrow width).

Based on field research it was concluded that all check dams in the watershed have fulfilled their function: they retained the maximum amount of sediment, stabilized the riverbed and reduced the bed slope. The siltation upstream of the check dams is covered with vegetation.

Considering that these check dams were constructed with the purpose of protecting the Ćelije reservoir from sedimentation, which is a large-scale water management facility with significant investment value, there is no need for a separate economic analysis to conclude that these check dams have justified their construction costs. It should be noted that no technical works have been carried out since the 1990s. The existing structures are not adequately maintained, and many of them are damaged.

3.1.2. Biological Works

In addition to technical erosion control works, activities related to afforestation and grassing have been carried out every year.

During the period from 1960 to 1980, an annual afforestation of about 100 to 150 hectares was carried out on these surfaces, resulting in afforestation of approximately 2000 to 3000 hectares for this period. In the Žunjska River basin, 150 hectares were afforested with black pine. From 1980 to 1988, a total of 1221 hectares were afforested, and after the year 2000, an average of 50 hectares were afforested annually.

Based on available data on afforestation works, it can be estimated that 5000 hectares have been afforested in the last 54 years. The tree species used for afforestation include spruce, white pine, black pine, acacia, poplar, red oak, douglasia, larch and Scots pine (Figure 8).

In parallel with afforestation activities, distribution of grass seeds to the local population was carried out to establish artificial pastures, thus preventing soil erosion.

The grass mixture was suitable for these terrains in terms of both soil and altitude. In most cases, the mixture consisted of the following grass species: Dactylis glomerata (orchard grass), Festuca pratense (meadow fescue), Festuca rubra (red fescue), Lolium multiflorum (Italian ryegrass), Lolium perenne (perennial ryegrass), Phleum pratense (timothy grass) and legumes: Trifolium pratense (red clover) and Lotus corniculatus (bird’s-foot trefoil).

Based on available data, it can be estimated that around 4000 hectares were seeded with grass.

3.2. Land Use before and after ECW

According to data from 1968, before the construction of the dam and the formation of the reservoir, forests covered 33.05% of the watershed. Meadows and pastures accounted for 26.84% of the area (

Table 6. Land Use in Watershed of the Celije Reservoir in 1968 and 2022.

Year	1968		2022
Land Use	Area km2	%	Area km2	%
Forests	201.35	33.05	384.19	63.07
Barren land	6.15	1.01	1.88	0.31
Degraded forests	97.25	15.96	54.55	8.96
Meadows and pastures	163.45	26.84	66.80	10.97
Arable lands	115.27	18.92	83.96	13.78
Settlements	25.68	4.22	14.53	2.38
Reservoir	-	-	3.25	0.53
Total	609.15	100	609.15	100

). Arable land covered 18.92% of the watershed. Areas under degraded forest accounted for almost 16%. Settlements comprised 4.22% of the territory. Barren land made up 1% of the territory (Figure 9a).

According to data from the land use map from 2022, forests cover 63.07% of the watershed area. These are mostly deciduous forests (56%), while coniferous forests cover 5% of the area. Degraded forests account for 8.96%. Meadows and pastures make up almost 11% of the watershed area. Arable land covers 13.78% and is located along the Rasina River and at higher altitudes in settlements. Settlements occupy 2.38% of the area. Barren land covers 0.3% of the territory (Figure 9b).

3.3. Spatial Distribution of Erosion Coefficient (Z), Specific Annual Gross Erosion (W) and Specific Sediment Transport (G) before and after ECW

According to the erosion map created for the Zapadna Morava watershed in 1968, before the construction of the dam and the formation of the Ćelije reservoir, 9% of the watershed area was extremely endangered by erosion (

Table 7. Surface areas according to the intensity of soil erosion in 1968 and 2022.

Year	1968		2022
Erosion Category	Area km2	%	Area km2	%
Excessive Erosion	54.74	8.99	0.03	0.00
Intensive Erosion	173.96	28.56	3.26	0.54
Medium Erosion	191.4	31.42	145.09	23.82
Weak Erosion	189.05	31.04	418.27	68.66
Very Weak Erosion	-	-	42.50	6.98
Total	609.15	100	609.15	100

). All categories of erosion were represented except for very weak erosion. Intensive erosion covered 28.56% of the watershed area, while medium and weak erosion were present on 31.42% and 31.04% of the watershed area, respectively. With a mean erosion coefficient of 0.62, the watershed area was affected by intensive erosion processes, indicating that the erosion intensity in the watershed was in the third category (Figure 10a).

According to the erosion map for the year 2022, the area under excessive erosion is almost negligible, with intensive erosion covering 0.54% of the area. Medium erosion is present on 23.82% of the area. Weak erosion dominates, covering 38.66% of the area, while very weak erosion is present on 6.98% of the area. The erosion coefficient is 0.35, indicating that the erosion intensity in the watershed is in category 2 (Figure 10b).

Specific annual gross erosion was approx. 1189 m³/km⁻²/year⁻¹ in 1968 and approx. 540 m³/km⁻²/year⁻¹ in 2022, while specific sediment transport into the reservoir was approx. 554 m³/km⁻²/year⁻¹ in 1968 and approx. 254 m³/km⁻²/year⁻¹ in 2022 (

Table 8. Specific annual gross erosion and sediment transport.

Specific Annual Gross Erosion		Specific Sediment Transport
m3/km−2/year−1		m3/km−2/year−1
1968	2022	1968	2022
1189.12	554.20	540.18	253.55

).

Total annual gross erosion was approx. 724,352 m³/year in 1968 and approx. 337,590 m³/year in 2022, while total sediment transport into the reservoir was approx. 329,000 m³/year in 1968 and approx. 154,500 m³/year in 2022 (

Table 9. Total annual gross erosion and sediment transport.

Total Annual Gross Erosion		Total Sediment Transport
m3/year		m3/year
1968	2022	1968	2022
724,351.84	337,589.51	329,052.69	154,450.81

).

4. Discussion

The results of erosion control works performed in the Celije reservoir watershed have decreased the intensity of erosion, sediment production and transport during the observed period. Biological works have altered the land use pattern in the watershed, impacting the intensity of erosion processes. According to the status as of 2022, the distribution of land use has changed significantly: forests now cover nearly two-thirds of the watershed.

The change in land use structure has led to a reduction in the intensity of erosion processes and a decrease of the erosion coefficient. The specific annual gross erosion in the Celije reservoir and the specific sediment transport decreased by 46%.

The effects of biological and biotechnical structures were more noticeable in later periods. Reforestation works are considered the most effective in reducing erosion intensity. The dependence of erosion coefficient reduction on the percentage of areas treated with biological works is significant, as these works directly affect changes in land use, thereby reducing erosion intensity. The results of the similar research in Serbia show the significance of demographic and land use changes in the control of the intensity of erosion. In the period between 1971 and 2016 the value of erosion coefficient in the territory of the Jablanica river basin decreased from 0.432 to 0.36 [69], while on the territory of the Grdelica Gorge, the value of erosion coefficient decreased from 0.84 to 0.32 in the period 1953–2016 [23].

According to Spalevic et al., land use changes in the last 50 years influenced a decrease in the soil erosion intensity for 14% in the Miocki Potok River Basin, Montenegro [70]. On the other hand, changes in land use can lead to increased erosion intensity, primarily changes in the structure of arable land [71]. Reforestation and improvement of scrublands and grasslands quality could be the first step in reducing sediment generation rate and water velocity and thus in decreasing sediment load of reservoirs. Increasing forest cover by 6%, combined with improvement of quality of the scrublands, in North Greece has led to a decrease in the mean annual sedimentation by 15% [72]. Forest restoration significantly improved water infiltration capacity in the soil. Results confirm that forest restoration plays an important role in the ecosystem services that regulate the hydrological conditions of the basins, as well as in the properties of the soil [73]. Afforestation of less-productive agricultural land provides many benefits, including soil erosion control, water retention, carbon sequestration and increased biodiversity [74]. Numerous studies on a global scale, using various models and research techniques, indicate the significant impact of Land Use and Land Cover Change (LULCC) on soil erosion. The erosion models can also be practiced in accordance with the climate change scenarios [25], especially for an efficient erosion management practice [75,76,77,78,79].

Most studies conducted in Serbia and the wider region show that in addition to the implemented erosion control measures, demographic changes in the watershed also had a significant impact on reducing soil erosion intensity, primarily due to the decrease in population and unfavorable age structure [23,24,27,69]. All of this has led to a reduction in agricultural production, resulting in the formation of grassland and forest cover on abandoned fields and orchards through the spontaneous growth of vegetation. This vegetation cover positively contributes to mitigating erosion processes.

Check dams have been used throughout the world for a variety of purposes, including torrent control, water supply enhancement, agricultural land development and watershed restoration [80,81,82,83,84,85]. These structures have a bed stabilization function, and contribute to hillslope consolidation, decreasing slopes, retention and sediment transport regulation. Some check dams have several functions at the same time [86]. Check dams are considered to be an efficient short-term solution for sediment control [87]. The research results on the effects of check dams show that these structures can effectively reduce the sediment transport capacity in the watershed [88,89,90].

Analyzing the effects of check dams from the perspective of retained sediment quantity showed that the greatest effect is achieved by selecting the appropriate location for constructing the object. Dams of the same effective height yield different effects depending on the valley width and bed slope.

The formed slope of siltation significantly influences the effects of dams. If the ratio of slope of siltation to natural bed slope is close to one, the silting up area is longer, and therefore the amount of retained sediment is greater (Figure 11).

Based on the research results, a statistical analysis was conducted to determine the relationship between the natural bed slope and the measured slope of siltation. The natural bed slope was defined based on the project documentation, while the slope of siltation was defined in the field. The Pearson’s correlation coefficient is 0.965, indicating a high correlation between the parameters. The coefficient of determination value explains 93.1% of the total variation. The dependence of the slope of siltation on the natural bed slope is represented by a linear regression Equation (9) in which the parameters are significant at the given level of significance:

I_z = 0.759 I_t − 0.794

(9)

Equation (9) cannot be applied for values of the natural bed slope less than 1.046% because negative values of the slope of siltation are obtained. It has a regional character and can be applied in torrential streams with similar natural characteristics.

This method was applied in research on the territory of Serbia for the Toplica river drainage basin [34], the torrents in Southeast Serbia [35,36], the Drina river drainage basin [37] and the Nišava river drainage basin [38].

However, despite research results showing that the erosion intensity in the watershed is satisfactory, it is important to consider that the physical characteristics as well as the geological and pedological substrate of the watershed make it highly susceptible to erosion processes. According to research, climate change from 1961 until today has shown significant changes in temperature and the distribution and quantity of precipitation [91].

The entire territory of Serbia has suffered the problem of the lack of implementation of erosion control measures and maintenance of existing structures since the early 1990s, caused by the economic crisis.

Erosion control management was implemented 50 years ago, but the approach to erosion protection has changed in the meantime, primarily due to a lack of financial resources, leading to inadequate implementation of erosion control measures. On a smaller portion of arable land, there are areas where erosion control practices are applied, but there are also areas where this principle is compromised, leading to an intensification of erosion processes.

5. Conclusions

The Ćelije reservoir was created in 1979 by building a dam on the Rasina River in the Zlatar Gorge. Significant erosion control works were carried out to protect the reservoir from sedimentation in the watershed. Research shows that these works have significantly reduced erosion intensity in the watershed, sediment production and transport to the reservoir, thereby extending the use of the reservoir.

The relationship between slope of siltation and natural bed slope obtained through research can be utilized in the design of projects within the study area, as well as in watersheds with similar natural characteristics.

Based on the current state, it may appear that erosion issues in the watershed have been minimized. However, the natural characteristics of the watershed make it highly susceptible to erosion processes. While erosion control management was implemented over 50 years ago, the approach to erosion protection has changed since the 1990s due to a severe economic crisis in Serbia, leading to reduced financial resources and inadequate implementation of erosion control measures. Maintenance of existing structures and the repair of damaged ones have also been neglected for the same reasons.

Considering the significance of the Ćelije reservoir for water supply, it is essential to continue with erosion control works and measures in the watershed. This is necessary not only to reduce the amount of sediment entering the reservoir but also to prevent chemical and mechanical pollution that affects water quality.