To What Extent Can a Sediment Yield Model Be Trusted? A Case Study from the Passaúna Catchment, Brazil

Abstract

Soil degradation and reservoir siltation are two of the major actual environmental, scientific, and engineering challenges. With the actual trend of world population increase, further pressure is expected on both water and soil systems around the world. Soil degradation and reservoir siltation are, however, strongly interlinked with the erosion processes that take place in the hydrological catchments, as both are consequences of these processes. Due to the spatial scale and duration of erosion events, the installation and operation of monitoring systems are rather cost- and time-consuming. Modeling is a feasible alternative for assessing the soil loss adequately. In this study, the possibility of adopting reservoir sediment stock as a validation measure for a monthly time-step sediment input model was investigated. For the assessment of sediment stock in the reservoir, the commercial free-fall penetrometer GraviProbe (GP) was used, while the calculation of sediment yield was calculated by combining a revised universal soil loss equation (RUSLE)-based model with a sediment delivery ratio model based on the connectivity approach. For the RUSLE factors, a combination of remote sensing, literature review, and conventional sampling was used. For calculation of the C Factor, satellite imagery from the Sentinel-2 platform was used. The C Factor was derived from an empirical approach by combining the normalized difference vegetation index (NDVI), the degree of soil sealing, and land-use/land-cover data. The key research objective of this study was to examine to what extent a reservoir can be used to validate a long-term erosion model, and to find out the limiting factors in this regard. Another focus was to assess the potential improvements in erosion modeling from the use of Sentinel-2 data. The use of such data showed good potential to improve the overall spatial and temporal performance of the model and also dictated further opportunities for using such types of model as reliable decision support systems for sustainable catchment management and reservoir protection measures.

1. Introduction

Soil is a dynamic system that is highly dependent on the variations of the surrounding environment. Erosion-induced changes are the dominant processes in terms of landscape and terrain shaping [1]. Erosion has multiple environmental and economic impacts. The first and most obvious impact is the degradation and productivity loss of fertile soils. Population growth goes hand in hand with a growth in food demand. The removal of the natural vegetation, deforestation, and the intensification of crop cultivation have increased the vulnerability of soil towards erosion [2,3]. Based on the results by [4], only during the last century, the per-capita removed soil has increased by around 400%. In comparison to 2000 years ago, the per-capita removed amount of soil today is around 2000% higher. In contrast, soil formation is extremely slow. Under tropical and temperate agricultural conditions, 200 to 1000 years are needed for the creation of 340 t ha⁻¹ of soil. The yearly renewal rate is around 0.2–2 t ha⁻¹ a⁻¹, while the soil loss in intense agricultural regions fluctuates from 10 to 100 t ha⁻¹ a⁻¹ [5]. A recent review study of [6] suggests a crop yield loss of up to 10% is to be expected by the year 2050 if the actual rates of soil loss continue. With such high differences between soil loss and renewal rates and also the high impact that soil loss has on the global food availability soil conservation practices become essential concerning the global food economy.

Water is the main natural erosive agent, as it is responsible for 80% of soil erosion worldwide [7]. Erosion has severe impacts on the aquatic ecosystems and water budget in reservoir systems. Sediment input and related nutrient flux due to erosion are the main factors deteriorating the water quality, threatening aquatic biodiversity, and reducing the lifetime of river impoundments. Therefore, soil loss is not just an issue concerning only food scarcity, but also water scarcity.

The cross-scale characteristics of the erosion phenomenon with its high spatial-temporal variation may cause high costs for the adequate quantification of soil loss by monitoring programs. Hence, alternatives like modeling are often considered for quantification of soil loss and localization of hotspots. A vast range of model types (physical, stochastic, or empirical) has been developed, but these models are normally specific to local or regional environmental conditions, and the performance varies based on the data availability and quality [8].

The aim of this study is to validate sediment input modeling by using validation measurements of sediment stock from the long-term siltation estimate in a reservoir. Large river impoundments represent the perfect opportunity as they often have trapping efficiencies >95% and consequently may serve as validation points for transported material [9,10,11]. For this study the sediment stock was assessed by high-resolution sediment magnitude measurements in the reservoir via a dynamic penetrometer [12]. For calculation of the revised universal soil loss equation (RUSLE) factors, satellite data were used for the land-use and land-cover (LULC) assessment, two soil sampling campaigns to define the soil properties of the area, and available datasets from the literature.

2. Materials and Methods

2.1. Study Area

The Passaúna Reservoir catchment (152.6 km²; 25°31′43″ S and 49°23′37″ W; 25°18′15″ S and 49°21′03″ W) is located near the Metropolitan Region of Curitiba and is part of a water supply system that provides water for more than three million people. About 30% of the population of the Metropolitan Region of Curitiba is supplied by this catchment. In 2001, the Environmental Protection Area of Passaúna was established, comprising 16,060 ha of territory Even so, anthropic pressure on the catchment has continued over the years (Figure 1). The Passaúna river composes 65.6% of the contribution area of the reservoir, followed by the contribution of the small sub-basins < 1 km² (8.4%), the Ferraria river (6.9%), the reservoir area (5.9%), the runoff lands around the reservoir (4.0%), the Eneas river (3.6%), and two other unnamed sub-basins with 3.2% and 2.6%, respectively [13].

Passaúna reservoir initiated operation in 1989. The total water surface area is 895 ha and the reservoir has an actual volume of 69.3 hm³, considering the spillway level at 887.2 masl. The intake is located approximately 3 km upstream of the dam. The impoundment structure is a 1200 m long and 17 m high rock-fill dam with a clayey core.

2.2. Sediment Yield Model

The sediment input (or sediment yield) is calculated as a product of soil loss from the hillslopes and a sediment delivery ratio:

SI = A·SDR

(1)

where SI stands for sediment input (or sediment yield), A for soil loss, and SDR for sediment delivery ratio.

As mentioned, the soil loss is calculated based on the RUSLE model. The universal soil loss equation (USLE) originated from [14] to assess the soil erosion in US agricultural land. Research for quantifying the soil loss started in 1940 in the Corn Belt and ended with the final publication by [14], where figures and relations were added for calculating each of the parameters. The next development in USLE happened in 1997, when [15] published the revised form of the Universal Soil Loss Equation (RUSLE). In the new version of RUSLE, the core philosophy of USLE was retained, even though significant changes in the calculation of the single parameters were included. The idea of USLE/RUSLE consists in the parametrization of the factors that affect erosion (terrain geometry, soil physical properties, rain characteristics, land use/land cover, and conservation practices).

In this study, due to the adequate data availability, a model in a monthly time resolution was used. Mathematically, RUSLE is presented in the following form:

A = L·S·R·C·K·P

(2)

where

• A is the soil loss at the investigated area
• L is the slope length factor
• S is the slope steepness factor
• R is the rainfall-runoff erosivity factor
• C is the cover management factor
• K is the soil erodibility factor
• P is the support practice factor

2.2.1. Topographic Factor LS

LS expresses the expected ratio of soil loss per unit area from a field slope to that from a 72.6 ft (22.13 m) length of uniform 9% slope under otherwise identical conditions [14]. The relation was adapted by [16], especially for the L Factor. The basis for calculation of the pixel-based topographic factor was a digital elevation model (DEM) of an accuracy of 10 m available from TanDEM-X service (Figure 2). For calculation of the LS Factor, the open source platform inVEST was used [17]. The LS Factor was calculated as follows:

LSi = (((A_i-in + D²)^m + 1 − A_i-in ^m + 1)/(D ^m + 2·x_i ^m·22.13 ^m))·S_i

(3)

where:

• S_i is the slope factor calculated from terrain slope $\theta$ in radians as showed below
  S = 10.8 sin θ + 0.03 when θ < 9%
  S = 16.8 sin θ − 0.50 when θ > 9%
• D is the gridcell dimension
• A_i-in is the contributing area (m²) at the inlet of a grid cell which is computed from the d-infinity flow direction method
• x_i = |sin a_i| + |cos a_i| when θ > 9%
  and a_i is the aspect direction for grid cell i
• m is the length exponent factor (

  Table 1. Values of dimensionless factor m.

  Slope % [s]	m
  s < 1	0.2
  1 < s < 3.5	0.3
  3.5 < s < 5	0.4
  5 < 9	0.5
  s > 9	m = β/((1 + β)) 1

  ¹ β = ((sinθ/0.0986))/((3∙sinθ ^0.8 + 0.56)).

  )

2.2.2. Soil Erodibility Factor K

The K Factor corresponds to the soil erodibility or the soil’s intrinsic susceptibility to erosion, which reflects the spatial variability of possible soil erosion depending on its structural and compositional characteristics [18]. This factor can be determined through experiments, and carried out in field plots using a specific measurement setup [19]. Alternatively, it may be obtained from predefined estimates based on the soil classes documented in the published literature reporting soil erodibility values for soil classes observed in different regions of Brazil (

Table 2. K Factor values from literature data base.

Soil Class	K Factor Value (t h MJ−1 mm−1)	Soil Class
Haplic Inceptisol	0.03	[20]
Humic Inceptisol	0.0175	[21]
Oxisol	0.018	[22]

).

In order to determine the K Factor, two soil sampling campaigns were organized in the Passaúna catchment with a total of 22 soil samples (Figure 3Left). The texture (silt, clay, and sand fractions) and loss on ignition at 550 °C (LOI550) were defined for each sample. For each point, three subsamples were taken as replicates within a radius of 5 m. Disturbed material was dried and sieved in 2 mm mesh, and the texture analysis was undertaken by the Bouyoucos hydrometer method [23] based on the classification by the United States Department of Agriculture (USDA), which addresses that the particle sizes between 0.05–2 mm are sand, between 0.002–0.05 mm are silt, and smaller than 0.002 mm are clay. For the samples of the first campaign, also some physical parameters of the soil were measured. All soil samples were used to calculate the K Factor at each location. For this study, Equation (4), proposed by [24] for the sample points collected covering Ultisol, Red Oxisol, and Typic Eutraquox classes, was used.

K = ((SAN + SIL)/CLA)/100

(4)

where SAN, SIL, and CLA are sand, silt, and clay fraction in percentage, respectively.

Afterwards, the values were interpolated using the inverse distance weighting (IDW) approach in order to obtain the information for the full coverage of the watershed.

2.2.3. Rain Erosivity Factor R

Based on the availability of data, two approaches to calculate the R Factor were investigated.

• Based on literature findings

For lack of 10–20 min frequency precipitation data for the Passaúna catchment, initially literature findings were used to determine the rainfall erosivity in the catchment [25] studied extensively the relations between the rain erosivity calculated from pluviographic and pluviometric data. Optimally, the rain erosivity is calculated by using long-term pluviographic (disdrometric) data, even though this type of data is mostly unavailable. The pluviometric data is often more easy to access but has a major disadvantage as it gives no information about the duration of the rain [25]; derived three different equations for three different locations in Parana to relate the erosivity calculated from the pluviometric data (R_Pm) with the erosivity calculated from the pluviographic data (R_Pg). Based on the aforementioned research [26], calculated the erosivity factor for the whole state of Paraná in a monthly resolution (Figure 4). In their research [26], integrated data from 114 pluviometric and pluviographic stations with more than 20 years of data (1986–2008).

The values used for this study were extracted from the monthly erosivity maps for the area of Curitiba. For the whole catchment with its 150 km², a constant value of R was used for each month.

2.

Based on Pluviometric Data of Daily Frequency

For calculation of the R Factor using the second approach, the data of two pluviometric stations in the catchment was used. The stations are part of the hydrological information system of Instituto das Águas do Paraná. The station of Colonia Dom Pedro was located in the central part of the catchment while the other station Barragem Sanepar (Dam), is located in the Southern part of the catchment near the dam (Figure 5). For both stations, precipitation data from 2000 until 2018 were available with a daily resolution.

For calculation of $EI$, thus the R Factor, the approach by [27] (Equation (5)) was applied. The precipitation patterns at both locations are similar; therefore, only one value of erosivity factor was used for the whole catchment (Figure 6).

EI = 68.730∙(C_c )^0.841

(5)

C_c = (p^2)/P

(6)

where p is the average monthly precipitation in mm and P is the yearly average precipitation in mm.

2.2.4. Cover and Management Factor C

The land-cover factor C is one of the most important factors, when it comes to what causes the highest inconsistencies in the outputs of a RUSLE-based model [28,29,30]. Optimally, the C Factor is determined from experimental soil erosion plots under natural rainfall conditions [31,32]. This type of data is often expensive to produce and in most cases, C Factors are derived from the literature. One of the most important drawbacks for the use of constant C Factors is the high variability of values for the same land-cover class among different literature sources. A literature review by [33] showed that the C Factors among the same class could differ by up to a factor of 100 (

Table 3. C-factor values for five land-use/land-cover (LULC) classes in Brazil from a literature review by [33].

Land Use	Cmax	Cmin	Cmin
Bare soil	1.000	0.696	0.100
Impervious areas	1.000	0.257	0.000
High vegetation	0.090	0.008	0.00004
Low vegetation	0.630	0.099	0.008
Water	0.000	0.000	0.000

). Another major disadvantage of constant C Factor values is the inability to capture the spatial and temporal variability of the factor values among the same LULC class. With the developments in satellite-based earth observation systems and the increase in data availability during the last decade, more scientists base their approaches on remote-sensing data [34,35,36].

For calculation of the C Factor in this study, the Sentinel-2 data was processed, and spatial information about LULC, urban soil sealing, and NDVI was derived.

For generation of the LULC maps, the Random Forest algorithm was used for pixel-wise labeling of a Sentinel-2 time series raster stack [37]. The scenes were selected based on image quality criteria and with the aim of representing different phenological phases. Train and test sample data were collected through visual interpretation of aerial images and field work. The estimate of overall accuracy based on a hold-out test set is 84%.

NDVI was the core parameter derived from the Sentinel-2 dataset. The use of NDVI values for calculation of the C Factors enabled a model setup with a monthly temporal resolution. NDVI is related to the vegetation density, biomass, and productivity [38]. It was calculated based on the 10 m red (Band 4) and near-infrared (Band 8) bands of Sentinel-2. An automated processing chain was established comprising the download, preprocessing (atmospheric correction), optimized cloud masking, scene selection, and processing of land surface variables. The automated processing was focused not only on NDVI but also on other variables like the degree of soil sealing or LULC.

Degree of soil sealing or imperviousness is defined as the fractional coverage of artificially sealed ground which impedes water from infiltrating into the ground. The calculation of imperviousness is based on a strong inverse relationship between vegetation cover and impervious surface as well as on the idea that an urban landscape can be linearly decomposed into vegetation, impervious layer, and soil [39,40]. The imperviousness layer was calculated based on a min-max rescaling of the NDVI derived from satellite acquisitions between the maturity and senescence onsets. The rescaling was guided by a visual comparison of results with submeter resolution aerial images as well as findings by [40], who studied the linear relationship between NDVI and imperviousness across several European cities.

Two NDVI-based approaches were considered for calculation of the C Factor in this study: [34,41]. As shown in [32], for the conditions prevailing in Brazil, the methodology derived from [34] (Equation (7)) produces more reliable results. Therefore, this approach was used for calculation of the C Factor.

C = (−NDVI + 1)/2

(7)

The previously mentioned satellite-derived data was used to calculate the C Factor also in non-sealed urban areas. As can be seen from Figure 7, in the urban areas, the NDVI is in the range of 0.25, which would result in a C Factor of 0.35–0.40, which corresponds to the C Factor values of arable land. Therefore, a filter was applied to the data with the simple logical condition that if a pixel in the urban areas had more than 60% soil sealing, the NDVI at the same location was set to 0.999, as it was assumed that no or very little sediment can occur from sealed areas.

2.2.5. Conservation Practices Factor P

During several field trips in the Passaúna catchment, many agricultural properties were visited. Support practice was observed at almost none of them (Figure 8). Therefore the P Factor was set to a constant value of $P=1$.

2.2.6. Sediment Delivery Model

The sediment delivery ratio plays a crucial role in the outcome of the final sediment amount reaching the river, as it is directly related to a large number of factors (amount of soil displacement, geometry of the transporting paths, land cover of the surrounding area, or amount of surface runoff) [42]. For this study, the SDR was calculated based on the flow connectivity approach by [43]. Hydrological connectivity is a term often used to describe the linkages between runoff and sediment generation in the upper parts of catchments and the receiving waters [44]. The use of the connectivity index as an input parameter for SDR has shown satisfying results globally [45,46,47,48]. For this study, the calculation of the connectivity index, and subsequently SDR, was implemented in ArcMap 10.5 as described in [49] based on the following formula:

SDR_k = SDR_max/(1 + exp((IC_0,k − IC_k)/K_IC,k))

(8)

where SDR_max is is the maximum attainable SDR coefficient at kth cell, set to 1, as soil in the Passaúna catchment has a high percentage of clay (0.002 mm), silt (0.002–0.05 mm), and fine sand (0.05–0.25 mm). IC_k is the index of connectivity as explained in [49], IC_0,k is a calibration parameter with a value of 0.5 [43,50] and K_IC,k is a calibration parameter with a value of 2.0 [43,50].

2.3. Sediment in the Reservoir

RUSLE results represent averaged long-term (mostly over decades) soil loss and sediment input (when RUSLE is combined with the SDR) from the catchment. When dealing with the results, the challenges are encountered mainly in having a reliable assessment of the sediment input in cases where no monitoring station is available. Due to their high trapping efficiency, large reservoirs with long residence times act as sinks for the incoming sediment. Therefore, they can be used as suitable long-term validation points for sediment input modeling. Several studies were conducted in this regard [9,51,52,53]. In order to acquire fast and accurate sediment information in areas where no previous data are available, a so-called portable free-fall penetrometer (PFFP) was used to assess the sediment distribution in the reservoir. PFFPs are not new in marine research (mostly sediment management in harbors), while their application in freshwater is still limited [54,55,56,57,58,59,60,61,62,63].

For this study, the commercial system GraviProbe (GP) produced by the Belgian company dotOcean (Figure 9) was used. Often, the spatial distribution of sediment thickness in a reservoir can be determined by bathymetric surveys, if precise pre-impoundment bathymetry is available [9,64,65]. Most of the time, this information is missing, and even when it is accessible, often the accuracy of the data does not allow for proper sediment stock estimation. One of the main advantages of the GP is its independence from previous data. The GP can deliver rapid results on penetration depth, cone penetration resistance, and shear strength of the sediment for each deployment. The GP was deployed at 134 points in the Passaúna reservoir (Figure 10). To determine the sediment magnitude, the information from the dynamic cone penetration resistance (DCPR) at some locations in the reservoir was related to the sediment magnitude data from core samples [12,66]. Characteristic changes in the sediment density/composition were identified and then related to the changes in the DCPR from the GP. The sediment thickness could be derived for all other points, where GP data was available because the relation between DCPR and sediment-pre-impoundment soil interface was established. This was possible due to the fact that the share of sand and coarse particles in the Passaúna sediment is smaller than 5% at most locations and a full penetration could be achieved. From the relation between core samples and GP information, it was observed that a DCPR of 200 kPa is the threshold between the sediment and the pre-impoundment soil. More information about the method can be found in [12].

3. Results

3.1. C Factor

For each of the available NDVI maps, the C factor was computed. The highest seasonal changes in the C Factor values were observed for the cropland (Figure 11). Between January and February, which is harvesting time, and November, which is seeding time, there is a change of almost 100% in the average C Factor of the catchment (from 0.15 to 0.28). A high interannual change in the C Factor was also observed in the scrubland/grassland areas. Winter and spring are characterized by a low vegetation coverage, whereas summer and partially autumn reveal a high vegetation coverage. Forests showed moderate changes mainly because a small percentage of the trees in humid subtropical regions lose their leaves in winter. The seasonal change in the forest average C Factor (change of <0.05 between maximum and minimum average C Factor) can also be related to the misclassification of certain areas with other LULC into forest class. Pasture and meadow follow also a similar land cover pattern. In summer and autumn, the vegetation cover is high, and in winter and spring, it diminishes. Bare soil has the smallest changes of all classes. There is a seasonal change of a maximum of 0.05 among the months and this can be attributed to the errors of the LULC classification process.

The difference among the seasons can be clearly observed also in the spatial distribution of the C Factor (Figure 12). The western area of the catchment, where most of the agriculture activity is located, shows higher values in July than in January. In July (winter period), the soil is mostly uncovered and has an average C Factor greater than 0.3. While in January (summer and wet season), vegetation covers most of the catchment area. Only sporadic parts of the agricultural areas, which were not seeded, had high C Factor values also in January.

3.2. R Factor

The R factor computed based on precipitation data showed results different from those calculated by [26]. The largest differences are observed in January, April, and October. The data by [26] shows high differences among the months and overestimates substantially for the month of January. The precipitation data from both pluviometric stations show a more uniform distribution than what is suggested by [26] (Figure 13). The calculations by [26] also include a margin of error due to the low density of weather stations. In certain regions, these coarser maps cannot represent accurately the rain erosivity when brought in a mesoscale plot. Therefore, for the final calculation of erosion and sediment input, the R Factor calculated from the pluviometric data in the Passaúna catchment was used.

3.3. K Factor

In general, the soil in the Passaúna catchment shows a low erodibility factor (<0.02 t ha h ha⁻¹ MJ⁻¹ mm⁻¹). The most erodible soils, dominated by Distrofic Latossol (Oxisol), are located in the northern part of the catchment, according to the soil map provided by the Brazilian Agricultural Corporation (EMBRAPA) [67]. The results are also aligned with further literature values in that geographic area, which also assessed that the K Factor for Latossol (Oxisol) is in the range of 0.019–0.026 t ha h ha⁻¹ MJ⁻¹ mm⁻¹ [68,69,70,71]. The western part of the catchment, which is also dominated by Oxisol, showed low soil erodibility with values reaching up to 0.013 t ha h ha⁻¹ MJ⁻¹ mm⁻¹.

In general, as shown in Figure 14, the soil has a similar texture pattern throughout the catchment area. The silt-clay content of the samples was always larger than 50%. The sand content in the soil is also relatively high (reaching up to 50% at some locations). Most of the catchment is covered in sandy clay, which has a low to average erodibility.

3.4. Sediment Delivery Ratio

Based on the physiographic characteristics of the Passaúna catchment, the SDR was calculated for each of the investigated months by applying the approach developed by [43] (Figure 15). In general, the calculated SDR could reach values of up to 0.15 in the dry months and rarely in some locations of above 0.15. The interannual vegetation cover which characterizes the region contributes to having low SDR values throughout the catchment. The highest SDR was observed in unprotected soil areas near the river stretches and at high slopes. The largest part of the catchment has SDR values lower than 7.5% in both dry and wet seasons. As explained in [49], connectivity, thus SDR, varies in both time and space. To define the change in the spatial patterns, the mean SDR was computed for each of the months. The results show low differences between the months. The mean SDR for the dry month of July was calculated to be 6%, while for the wet month of January, it was 5%.

3.5. Sediment Input—Initial Model Run

The results show high sediment input in all of the wet months. The highest sediment input happens in January with 14,000 t, even though the vegetation coverage of the catchment is rather high. In the month of August, despite the low vegetation cover, the overall sediment input from the catchment is the lowest (Figure 16). By comparing the soil loss distribution to the LULC map (Figure 17), it can be observed that high sediment input occurs from the forested areas. Even in the months of winter, when precipitation is low, there is significant sediment input from the forested areas. By comparing the mean value of the calculated C Factor for the forest areas in the Passaúna catchment with literature values, it was found that the assessed C Factor is significantly overestimated (Figure 18) from the method used in this paper. The average C Factor found by [33] (Table 3) is between 10 to 20 times lower than the calculated C Factor values during the months of July and October (Figure 18) from this study. The values of Max. and Average in Figure 18 refer to the maximum C Factor found in the literature review (Min = 0). The approach developed by [34] seems to overestimate the C Factor in forested areas, even though it is not sure from their research whether the empirical approach they developed can be used in forested areas. Therefore, an arbitrary correction factor of 0.05 was applied to the C Factor by multiplication. This value of 0.05 was chosen, as the calculated values of C Factor were 10 to 20 times higher than the average C Factor of forest areas in that region. Furthermore, the new R Factor, calculated from the pluviometric data, was also included in the new equation for giving the final results shown in Figure 19 and Figure 20.

The C Factor correction decreased the overall amount of sediment input into the Passaúna reservoir by 30% from an initial 94,300 t a⁻¹ to 57,300 t a⁻¹. After the inclusion of the new R Factor calculated from the daily precipitation data, the sediment input decreased by a further 5% to 54,800 t a⁻¹ (Figure 20b). The use of the new R Factor shifted also the seasonal dynamics of the sediment input. The final model indicates that the most important month in terms of sediment yield is not January but October (Figure 20a). The month with the lowest input is April and not August, as suggested by the initial model results. The spatial distribution of sediment input changes significantly between the final model run (C and R Factor correction) and the initial model run (Figure 19). The sediment input from forested areas is reduced substantially in the final model run to less than 0.1 t ha⁻¹. For the overall operational time of the Passaúna reservoir (30 years), according to the modeling results, the accumulated sediment stock should be approximately 1.6 × 10⁶ t.

3.6. Reservoir Sediment Stock

As explained in [12], a DCPR of 200 kPa is defined as the vertical consolidation threshold between the sediment overlay and pre-impoundment soil. The sediment in the reservoir showed a high spatial heterogeneity of the siltation patterns. Even points at a horizontal distance of 10 m showed different sediment thickness values, mainly because of the bottom topography. This underlines the need for large numbers of measurement points, to obtain representative estimates of the accumulated sediment. Near the deepest part of the reservoir, a sediment thickness of up to 1.8 m could be observed. The areas with the highest sediment accumulation are located near the dam and near the inflow (Figure 21). Also, the sidearm located in the southwestern part of the reservoir showed high sedimentation rates compared to the northern areas of the reservoir. By applying the inverse distance weighting interpolation technique, the spatial distribution of the sediment magnitude was obtained. In total, a stock of 3.4 × 10⁶ m³ of sediment could be measured in the reservoir. According to [66], the sediment has an average density of 1.12 g/cm³. Therefore, the total mass of sediment in the reservoir is approximately 3.8 × 10⁶ t.

4. Discussion

4.1. Comparison of the Approach to Literature Findings

Several studies were conducted in the Alto Iguacu area in regard to soil erosion [71,72]. Reference [72] conducted a similar study in the Passaúna catchment even though the methodology followed to calculate erosion was different. The soil loss and sediment input were calculated in a yearly time step and the calculation of the C Factor was based only on a LULC map. Despite the similarities in the spatial distribution patterns, the findings from our study indicate that the soil loss is lower than the amount calculated by [72] (

Table 4. Comparison of results by this study with results by [72].

Soil Erosion Classes (Annual Mean)	Present Study (%)	[72] (%)
Very Slight (< 2 t ha−1 a−1)	55	52.0
Slight (2–5 t ha−1 a−1)	3.5
Moderate (5–10 t ha−1 a−1)	3.7
High (10–50 t ha−1 a−1)	15.8	10.0
Severe (50–100 t ha−1 a−1)	9.0	5.0
Very Severe (100–500 t ha−1 a−1)	11.3	33.0
Catastrophic (>500 t ha−1 a−1)	1.4

). Our results show that almost 63% of the catchment had very slight, slight, or moderate soil loss, against 52% found by [72]. Major differences were also observed in areas with very severe and catastrophic soil loss. [72] calculated that 33% of the catchment had more than 100 t ha⁻¹ a⁻¹ of soil loss, while our results showed that only 12.7% of the catchment had more than 100 t ha⁻¹ a⁻¹.

Reference [73] conducted another study in the Passaúna watershed, but focused mostly on the continuous monitoring of suspended solids in the Passaúna river before entering the reservoir. In his study, [73] collected 33 large-volume river samples between February 2018 and July 2019. In his study, also measurements from one intensively measured high-flow event of October 2018 were included. The point where the measurements were conducted collects water from 55% of the overall Passaúna Reservoir catchment. For this case, [73] calculated an annual average flux of 10,800 t a⁻¹. This value is approximately 300% lower than the value calculated for sediment input from 55% of the catchment from this study. Reference [73] explains the relatively low flux values with the importance of episodic high-flow events, whose dynamics are not properly described and, in this case, are strongly underestimated by the derived rating curves of suspended solids.

Other regional studies such as that by [71] or the more holistic study by [36] show similar patterns of soil loss in the area of Parana and Alto Iguacu. However, the information presented in these studies is too coarse and cannot be directly compared with our findings.

4.2. Sediment Input from Catchment vs. Reservoir Sediment Stock

By comparing the results from the two approaches, it could be observed that the sediment stock is 229% higher than the overall sediment input from the catchment, as calculated from the model (Figure 22). The discrepancies in the results of the modeling are rather high. However, when we refer to the sediment stock, all the material entering the reservoir is included, including here the organic and mineral material that was inside the reservoir before impoundment or which was deposited during the construction phase of the reservoir. In addition, based on the definition by [14], USLE (and subsequently RUSLE) accounts only for the sheet and rill fractions of the soil loss. Therefore, further factors have to be considered to reach a complete estimate of sediment input from the catchment (

Table 5. Overview of factors creating inconsistencies.

	Factors Creating Errors
In reservoir	Internal production
	Existing biological stock
	Errors of the measuring concept
	Trapping efficiency of reservoir
In catchment	Errors associated with RUSLE calculations
	Errors associated with SDR calculations
	Non-inclusion of gully erosion in RUSLE
	Non-inclusion of channel erosion in RUSLE

).

One important factor that can create bias in the sediment budget is the error created from the measuring and processing technique in the assessment of reservoir sedimentation. After the interpolation, the frequency distribution of sediment magnitude values changes significantly as shown in the pie charts of Figure 21a,b. This indicates that the interpolation technique has a significant effect on the overall results. The average sediment magnitude of the raster is 40 cm, which is 30% smaller than the average of all measurements (57 cm). An underestimation of the average value from the interpolation technique shows an underestimation of the calculated sediment volume.

In order to properly compare the interpolated map with the measured values, the spatial component should also be taken into consideration. This means that if most of the measurements are located in the thalweg (disproportionally with its surface compared also to the bank slope areas), the average value for the measurements will be higher than the average from the interpolated values, as most of the accumulation is expected to be in the thalweg. Therefore, the reservoir was divided into two parts, thalweg and reservoir bank slope, as shown in Figure 23. For each of the compartments, the average sediment thickness from the GP measurements was calculated. Finally, an overall average value for the whole reservoir was calculated as shown in Equation (9).

M = (M_t∙A_t + M_b∙A_b)/(A_t + A_b),

(9)

where M_t and M_b are the averages of the measurements in the thalweg and reservoir bank respectively, while A_t and A_b are the areas of the aforementioned compartments.

Based on Equation (9), the average sediment magnitude measured in the reservoir is 62 cm, thus 36% higher than the mean raster average (40 cm). Hence, the interpolation can lead to an underestimation of the sediment stock of up to 36%.

Another important factor, which can affect the sediment balance in the Passaúna reservoir, is the contribution of internal production to the sediment stock. Apart from acting as a sink, the reservoir acts also as a source of particles. Due to the climatic conditions and the relatively high nutrient availability (mesotrophic state), the reservoir is productive for plankton communities. Therefore, the autochthonous material created in the reservoir can play an important role in the sediment balance of the system. In other studies, it was observed that the autochthonous material can account for up to 75% of the sediment stock [74]. Even if the exact share of internally produced sediment cannot be defined, the high LOI (>20%) in the main basin of the reservoir in comparison to average values of ~10% at the inflow underlines the importance of autochthonous production. Moreover, before flooding, the reservoir area was not cleaned from the existing biomass. Several trees and former vegetation areas are still visible at the reservoir bottom (Figure 24). This organic material (sometimes degraded) also plays its role in the bias created when comparing both approaches.

One of the most discussed limitations of RUSLE is its inability to represent also gully and stream bank erosion [75,76,77]. Even with the calculation of the SDR based on the connectivity index, the uncertainties about the prediction of gully and streambank erosion are still present. In comparison to sheet and rill erosion, gully erosion is generally less investigated. However, various studies [78,79,80] showed that gullies substantially contribute to the sediment budget at a catchment scale. They do not only contribute as a sediment source, but also increase the efficiency of sediment transport from uplands to the valley bottom and river channels, as most of the sediments generated from rill and inter-rill erosion that are not connected to gully structures are deposited at the foot of the hillslopes [81].

Reference [82] estimated that 47–83% of the sediment occurred from gully erosion. In addition, [81] indicated in a review study, that worldwide, gullies can represent 10–94% of the total sediment yield from water erosion. When referring to the sediment input from the catchment, it is still unknown to what extent the gully structures contribute to this budget for the case of Passaúna.

4.3. Limitations of Normalized Difference Vegetation Index (NDVI)-Based Approaches for the Estimation of the C Factor in Forested Areas

When water is the erosive agent, there are three main phases that characterize erosion. The first phase is the detachment of soil particles. In this phase, the potential energy of the raindrops due to its absolute elevation is transformed into kinetic energy. The free fall of the raindrops due to gravity causes remobilization of soil particles when the drops hits the soil surface. The second phase is the transport of the detached material by the accumulated flow, and the final phase of erosion is deposition, which occurs when the transport forces are depleted [14,83,84]. In the C Factor results before correction, a similarity in the values of plant-covered arable land and forest areas was observed. Despite the similarities in the cover canopy between planted arable land and forest (according to the NDVI values), the topsoil’s physical properties between these two classes are completely different. While in the erosion component associated with the rain splash, both LULC classes behave similarly due to the similar protection by the plant canopy, in the component of erosion associated with runoff, forest and arable land behave differently. The soil surface below the plant coverage in arable lands is basically bare and facilitates the detachment of soil particles from surface runoff. In contrast to this, the soil in intact forests is normally covered by low vegetation (grass or meadows) and leaf litter, which hinders the creation of runoff and reduces the soil particle detachment. In addition, the soil is more compact in forested areas than in arable land where tillage almost always takes place. In its original form, the C Factor has a direct relation to the soil loss ratio (SLR) [15]. The SLR is a product of five sub-factors, which are prior land use, canopy cover, surface cover, surface roughness, and soil moisture. All of the former factors, except the canopy cover, are associated with the conditions of the soil surface, indicating the importance of the top soil conditions for soil movement initiation. Therefore, the C Factor cannot be quantified only by taking into consideration the vegetation index (canopy cover) but should also include the properties of the soil surface, especially in non-agricultural areas [85,86,87,88].

4.4. Management Implications

In terms of management, an analysis was performed on how the soil loss and sediment input from the catchment could be reduced by afforestation in the most problematic areas characterized by high soil loss rates. For this purpose, three scenarios were investigated; more specifically, afforestation of areas with more than 100, 200, and 250 t ha⁻¹ a⁻¹ of soil loss, which account respectively for 12%, 5%, and 3% of the catchment area (Figure 25). From our calculations, it was found that with full afforestation of these areas, a reduction of 50%, 27%, and 26% of the annual sediment input could be achieved. Such a measure, apart from tackling the soil degradation in the catchment, can also contribute to significantly increasing the reservoir lifetime.

October and September are the most important months in regard to sediment input and soil loss (Figure 20). For October in particular, the combination of the RUSLE factors is the most effective for producing the highest amount of soil loss. Figure 26 shows the combination of C and R Factors for the three most characteristic months of the year. In the case of the C Factor, October has the same values as July, which is one of the driest and coldest months of the year and has the lowest vegetation cover. As far as the R Factor is concerned, the erosivity is as high as the erosivity in the month of January, which is the month with the highest rainfall. In the case of October, the worst possible combination is present as the rainfall erosivity is maximal, while the vegetation cover is minimal. This combination of factors produces the highest soil loss from a system. In the case of proper land management implementation, like crop rotation, applying crop residues and cover crops in the unprotected soil during the winter and spring months (April–October), a significant reduction of the sediment input could be achieved [89,90,91].

4.5. Uncertainties of the Sediment Yield Model for the Passaúna Catchment

RUSLE was developed as a tool for long-term soil loss calculation. By calculating the C Factor from a certain scene in 2017 or 2018, it is assumed that the LULC of that specific month has not changed during the last 20 years (rain data available for approx. 20 years). This is to a certain extent not correct. In Parana state, from 1990 until 2019, there has been an increase of almost 45% in arable land and a 5% annual increase in urban areas [92]. Most of this area that was transformed into agricultural land used to be forest, which suggests a gradual increase in erosion in the last 20 years.

This is one of the major drawbacks of the method. However, this drawback can also represent an opportunity. In the case of available precipitation and NDVI data for single months for the entire investigated period, RUSLE could be adapted also for calculation of the actual sediment input and soil loss from that certain month of that specific year. In this way, a calibrated model could be used to derive an accurate balance of sediment input for each month and not only a long-term monthly average of sediment input as in most applications so far.

4.6. Benefits from the Integration of Sentinel-2 Data in Erosion Modeling

The use of vegetation index for calculation of the land cover factor from freely available data is not new. Several studies were conducted based on this principle. However, the spatial and temporal resolution of the images (Landsat or MODIS) in most of the existing literature is relatively low (around 30–250 m) compared to the also freely available Sentinel-2 data [32,93,94,95,96]. Improved spatial accuracy and flyover frequency of satellite imagery leads to better erosion modeling results [97,98]. In the tropics and sub-tropics, it is likely that a sequence of satellite scenes show high cloud cover, leading to large data gaps. This emphasizes the importance of short flyover intervals in order to represent fast-changing conditions in the catchment. Furthermore, by the application of more advanced processing steps, more specific information can be derived about the investigated area (for example, the degree of soil sealing). Certain information can be used, as in the case of this study, for a better mapping of erosion and sediment input from urban or semi-urban areas.

Despite the prevailing discrepancy between the sediment yield and accumulated mass inside the reservoir, the model is fully capable of representing the spatial and temporal patterns of soil loss and sediment yield. The application of RUSLE-based models in monthly resolution as decision support systems, due to the increased performance of the model in both spatial and temporal dimension, can lead to an improved river basin management. Specifically, as shown in the previous section, the inclusion of NDVI data can enable the planning of management activities not only in high spatial accuracy but also aimed at the temporal dimension by highlighting the most problematic months of the year. The latter is of high importance as it may lead to reduced catchment management costs.

5. Conclusions

In this study, the modeled sediment yield from a catchment with the sediment stock in a reservoir at the outlet of the catchment is compared. For assessment of the sediment yield from the catchment, a RUSLE-based model was used. The error margins of the RUSLE model results can fluctuate considerably. Therefore, the models have to be coupled with validation measures. This study shows that reservoirs can be used as validation points, despite some limitations. The use of reservoirs as validation points represents a good opportunity, as they collect almost entirely the incoming sediment. Reservoir sediment stock measurements are often easier to achieve than conventional continuous sediment flux monitoring, which produce a high sampling effort and need to deal also with large errors due to the high variability in the river stretches. The assessment of the siltation status of the reservoir creates added value for every operator. In the case of complex systems, however, as shown in this study, several other factors can affect the reservoir sediment balance and, therefore, be misleading regarding the aim of the research. Reservoirs of lower process complexity (e.g., in mountainous areas, low organic material input, or low temperatures) can be more easily used as validation points.

This study showed that the most important factors that create discrepancies for the case of the Passaúna sediment budget are associated mostly with the sediment yield model. On the other hand, when including the errors because of the interpolation technique, the underestimation of sediment yield from the model may become even greater. Even though we fully agree that a RUSLE-based model can reproduce the spatial and temporal patterns of sediment yield from a catchment, the comparison of the approaches in this study shows that there are clear limitations of using modeling approaches for reservoir sediment stock or reservoir lifetime assessment. The two components of the model (RUSLE and SDR model) do not allow the usage of only one calibration point as it is impossible to track the source of error. In order to increase the accuracy of the results, models including the effects of channel, gully, and artificial ditches are needed, and these models need to be calibrated with complementary measures (river-suspended solids and bedload monitoring, calibration of models with erosion plots, or quantification of gully erosion).