Assessment of uncertainties in a complex modeling chain for predicting reservoir sedimentation under changing climate

Long-term predictions of reservoir sedimentation require an objective consideration of the preceding catchment processes. In this study, we apply a complex modeling chain to predict sedimentation processes in the Banja reservoir (Albania). The modeling chain consists of the water balance model WaSiM, the soil erosion and sediment transport model combination RUSLE-SEDD, and the 3d hydro-morphodynamic reservoir model SSIIM2 to accurately represent all relevant physical processes. Furthermore, an ensemble of climate models is used to analyze future scenarios. Although the capabilities of each model enable us to obtain satisfying results, the propagation of uncertainties in the modeling chain cannot be neglected. Hence, approximate model parameter uncertainties are quantified with the First-Order Second-Moment (FOSM) method. Another source of uncertainty for long-term predictions is the spread of climate projections. Thus, we compared both sources of uncertainties and found that the uncertainties generated by climate projections are 408% (for runoff), 539% (for sediment yield), and 272% (for bed elevation in the reservoir) larger than the model parameter uncertainties. We conclude that (i) FOSM is a suitable method for quantifying approximate parameter uncertainties in a complex modeling chain, (ii) the model parameter uncertainties are smaller than the spread of climate projections, and (iii) these uncertainties are of the same order of magnitude as the change signal for the investigated low-emission scenario. Thus, the proposed method might support modelers to communicate different sources of uncertainty in complex modeling chains, including climate impact models.


Introduction
Albania possesses adequate conditions for hydroelectricity production because of its location in the Balkans and its mountainous topography.Considering that the country's power supply is, with a value of 99%, almost exclusively based on hydroelectric power (IEA 2022; Lehner et al. 2005;Statkraft 2019), the Devoll River offers great potential for hydropower development due to its large streamflow.However, as a result of the active erosion processes taking place in the catchment and the consequent transport of sediments, it is also considered the most turbid river that drains into the Mediterranean Sea (Ardıçlıoğlu et al. 2011).Constructed reservoirs along the river interrupt the sediment continuum, resulting in the deposition of sediments and progressive reservoir sedimentation.
In addition, climate change may not only directly (e.g., change in temperature and precipitation) but also indirectly (e.g., change in land use) influence erosion processes in the catchment, resulting in higher sediment loads and amplified sedimentation processes in these reservoirs (e.g., Plate 1993;Walling and Fang 2003).Therefore, the lifetime of planned and constructed reservoirs and the efficiency of hydropower production may decrease (Mahmood 1987).These trends are also predicted for other areas and catchments within Europe.According to Wagner et al. (2017), the average annual electricity generation for the Alpine region will slightly decrease by the year 2050, due to the effects of climate change.Additionally, Panagos et al. (2021) stated that climate change is the main cause of the increase in mean soil erosion rates on agricultural land in Europe by the year 2050.Hence, climate change will have negative impacts on soil erosion, especially 1 3 in regions where this process can already be considered critical, such as the Devoll River basin (Li and Fang 2016).
While several studies show an expected annual increase in sediment yield due to future climate change scenarios (Azari et al. 2016;Chen et al. 2020;Li et al. 2022), others indicate that sediment yield could even decrease in some regions (Bussi et al. 2014;Hirschberg et al. 2021).These contrasting findings suggest that the behavior of hydrological and morphological processes under the impacts of climate change depends not only on regional characteristics but also on the models involved, especially the climate models that provide the forcing data.As part of the DIRT-X project (https:// dirtx-reser voirs 4futu re.eu/), the dynamics of the different hydrological and morphological processes in the Devoll catchment are identified and their response to climate change is quantified.Hence, a special focus is given to the reservoir inflow (runoff and transported sediments) and sedimentation processes within the reservoir.Since the latter relies on the behavior of the inflow boundaries, a modeling chain is required to predict the future development of all involved variables and finally to predict the bed level changes within the reservoir (target variable).
Even though several studies were carried out for different regions worldwide to predict the erosion processes in the catchment and the resulting sediment yield, they all focused on a single model.In addition, they focused mainly on the catchment and no hydro-morphological processes in the river or reservoir were considered.For example, Shrestha et al. (2013), Azari et al. (2016), Zettam et al. (2017), Santos et al. (2021) and Li et al. (2022) employed the Soil and Water Assessment Tool (SWAT) to predict the sediment yield coming from catchments in Laos, Iran, Algeria, Brazil, and China, respectively.In the study carried out by Bronstert et al. (2014), the Water Availability in Semi-arid Areas with SEdiment Dynamics (WASA-SED) model was implemented to predict water and sediment fluxes in semi-arid environments.Other examples can be found in Nerantzaki et al. (2015) and Nunes et al. (2013), who also applied the SWAT model in Mediterranean catchments for predicting suspended sediment transport and erosion dynamics, respectively.Although the latter applied a chain of models combining SWAT with a physically based distributed erosion model, reservoir sedimentation processes were not included.Wagner et al. (2017) also focused on a modeling chain applied to an Alpine catchment.However, the modeling chain consists only of a hydrological and a hydropower model, thus the focus is not on intermediate erosion and sediment transport processes as in this study.More recently, Wild et al. (2021) developed a Python-based framework to simulate runoff, sediment, and hydropower production.Although several processes are considered in the model, the focus is on decision-making and the evaluation of possible reservoir configurations.Furthermore, the model has some limitations regarding the representation of physical processes (e.g.runoff generation), where it still relies on the output of other external models.
The novelty of our study is the development and application of a process-based modeling chain composed of three different models that aim to predict the sedimentation processes in the Banja reservoir under changing hydro-climatic conditions.With this modeling chain, we ensure a detailed representation of the physical processes leading to reservoir sedimentation, while exploiting the capabilities of state-ofthe-art models tailored to particular processes.Since each model works independently and has different input and output variables (e.g., runoff or sediment yield), the subsequent models rely on accurate output variables to ensure the applicability of the modeling chain for predicting bed elevation as the final target variable.
However, when more than one model is involved in predicting a target variable, superposition effects of uncertainties from different sources may occur, resulting in a propagation and an increase in the uncertainty of the final target variable.Since the selection of model parameters and their associated values might be challenging, e.g., due to a lack of measurements, it is of interest to know not only their impact on the final simulation results but also the confidence of these results (Moges et al. 2021).Other types of uncertainties are related to the Global Climate Models (GCMs), to the downscaling techniques used for linking the large scale of GCMs to the regional models (Regional Climate Models, RCM), and finally to the model scale (Prudhomme and Davies 2009).
The significance of the aforementioned uncertainties is investigated and presented in this study.The question that arises at this point is: Are simulation results more affected by perturbations in the model parameters (parameter uncertainty) than by the spread of climate projections (climate model uncertainty)?To answer this question, approximate uncertainties related to model parameters are calculated using a simplified method and compared to the spread of climate projections by analyzing three GCM/RCMs and the Representative Concentration Pathway RCP2.6.The selection of RCP2.6 is motivated by the fact that this (mitigation) scenario is subject to the smallest change signals amongst all available emission scenarios.Consequently, results from parameterizations with only small variations are not compared to higher change signals that would result from high emission climate scenarios.

Study area
The study area is located in the Devoll River catchment, upstream of the Banja reservoir in Albania, approx.70 km south of the capital city Tirana.The catchment covers a surface area of 3140 km 2 with varying topography and altitudes, ranging from 100 to 2000 m a.s.l. Figure 1a shows the location of the study area within Albania and the catchment area.In this figure, the catchment is subdivided into two subcatchments (Fig. 1c), which were delineated according to the topography and location of the gauging stations (Banja and Kokel).The main city located in the study area is Korça and is also shown in the figure .Since Albania belongs to the Mediterranean climatic belt, the climate in the study area is mainly characterized by dry and hot summers and mild, rainy winters (Eftimi 2010).The mean annual temperature in Korça is 10.3 °C, with a mean value of 19.9 °C during summer (July) and 0.9 °C during winter (January) (climate-data.org2019).The mean annual precipitation in the highland plain near Korça is around 800 mm yr − 1 while the Western part receives up to 1600 mm yr − 1 (Almestad 2015;climate-data.org 2019;Eftimi 2010).However, in higher altitudes, snowfall is common during the winter months.Snow cover depths and days with snow cover vary strongly, depending on the location within the catchment (Mouris et al. 2022).Although the majority of the study area is forested (30%) and covered by scrubs and herbaceous vegetation (25%), agriculture predominates in the Korça plain.Since parts of the catchment are characterized by barren and steep slopes, loose soil results in high erosion and subsequently high sediment loads entering the Devoll River from the catchment.Hence, reservoir sedimentation will be a severe challenge for existing and planned reservoirs along the Devoll River.The river is dammed approximately 160 km from its source, forming the Banja reservoir located near the town Gramsh.The embankment dam with a clay core has a maximum height of 80 m (between 95 and 175 m a.s.l.) and was impounded in 2016, with a maximum storage capacity of the reservoir of approx.400 million m 3 .

Meteorological data
Within the study area, four meteorological stations are in operation and record daily values of precipitation, temperature, and wind speed.Measurements are available for the period from 09/2015 to 08/2020.Other variables, such as radiation and relative humidity, are not recorded by these stations.Thus, the prediction of the hydrological response of the catchment becomes quite challenging, since a longterm time series of all meteorological variables would be required.As a viable alternative, the ERA5 reanalysis dataset is used for the hydrological simulations.This dataset is available with an hourly resolution on a grid of approx.31 × 31 km since 1959 onwards (Copernicus Climate Change Service 2017).

Measured data
Hourly measurements of runoff are available at the Kokel gauging station, covering the period from 03/2016 to 04/2018.Furthermore, suspended sediment concentration measurements are available at the same gauging station.The runoff and the suspended sediment concentrations are obtained by using the acoustic backscatter signal from two side-looking H-ADCPs (Horizontal-Acoustic Doppler Current Profiler; 0.6 and 1.2 MHz).The approach used to calculate suspended sediment concentrations based on acoustic backscatter data for this study site is described in Aleixo et al. (2020).Nevertheless, these measurements are only available for water depths at the gauging station exceeding 1 m.In the final step, the suspended sediment load is calculated by using the measured suspended sediment concentration and the associated runoff.
In addition, there are Digital Elevation Models (DEM) available from two bathymetric surveys of the Banja reservoir.The first survey was carried out in 2016, shortly before the impoundment of the reservoir.It was a drone survey of the terrain and a subsequent structure-from-motion postprocessing.The second survey was conducted in 2019, after 3 years of operation, by moving ADCP measurements and was used to calibrate the reservoir model.A digital elevation model of differences (DoD) of the two surveys shows a general sedimentation trend, with an average deposition height of 2.7 m in the upstream part (> 5 km distance to the dam) of the reservoir.

Modeling chain
The prediction of sedimentation processes within a reservoir involves several preceding catchment processes, which need to be considered in the simulations.To tackle these processes in a reliable manner, hydrological, soil erosion, sediment transport, and hydro-morphodynamic models are necessary.In some cases, some of these processes can be simulated in a simplified way by a single model (e.g. with SWAT).However, limitations often arise regarding the representation of physical processes, and this must be considered when analyzing the results.Due to the progressive development of modeling tools and their specialization in certain processes, the use of multiple state-of-the-art models in a modeling chain seems to be a promising approach to increase the quality of simulation results.For these reasons, a modeling chain composed of different and independent models is applied in this study to benefit from several models with a specialization in particular processes.
The schematic modeling chain is depicted in Fig. 2.Besides the models used, the target variables (output variables) of each model are shown, which then serve as input for the subsequent model.For example, the target variable of the soil erosion and sediment transport model (RUSLE-SEDD) is the suspended sediment load, which serves as input for the hydro-morphodynamic reservoir model (SSIIM 2).The final target variable of the modeling chain is the bed elevation along the thalweg of the Banja reservoir.A description of each model is presented in the following section.

Model setups
This section summarizes the main processes included in the three models used in the modeling chain for the Devoll catchment and the Banja reservoir.

Water balance model
The hydrological processes are simulated by the Water Flow and Balance Simulation Model (WaSiM, Schulla 1997, 2021).It is a physically based distributed model capable of representing the water cycle above and below the land surface.WaSiM uses physically based modeling approaches for the simulation of the different hydrological components (Schulla 2021).In this study, the Richards version 10.04.07 is used, including the most recent snow canopy interception sub-model (Förster et al. 2018).The model domain has a spatial and temporal resolution of 1 km 2 and 3 h, respectively.The calibration period spans from 05/2016 to 04/2018, for which measured runoff is available.A first year (05/2015-04/2016) is considered as a warm-up period.
Table 1 summarizes the main processes involved in WaSiM and the selected methods for obtaining the values for each of them.As the ERA5 dataset is available with a spatial resolution coarser than the model grid, the values are interpolated using the methods described in Table 1.

Soil erosion and sediment transport model
Soil loss and sediment transport are calculated at the catchment scale with the Revised Universal Soil Loss Equation (RUSLE) model (Renard 1997) in combination with the Sediment Delivery Distributed (SEDD) model (Ferro and Porto 2000).The RUSLE calculates the gross soil erosion in the catchment, while the SEDD model estimates sediment transport and delivery.The model is spatially discretized into cells of 25 m x 25 m.Since this model combination enables the estimation of a monthly or annual suspended sediment load for any point in the river network, the sediment input to the reservoir is computed.Monthly suspended sediment loads from 05/2016 to 04/2018 are used for the calibration of the presented RUSLE-SEDD model.
The soil loss A [t ha − 1 yr − 1 ] is determined as a multiplication of six erosion risk factors (Eq.( 1)), which are summarized in Table 2.More detailed information on the input datasets, applied methods, and codes can be found in Mouris et al. (2022).

Reservoir model
The fully 3d numerical model Sediment Simulation In Intakes with Multiblock Option (SSIIM 2) is used to simulate flow characteristics, suspended sediment transport, (1)

Uncertainty quantification of model parameters
Several types of uncertainties are expected in modeling, arising from the complex behavior of environmental systems, simplifications in models, unknown boundaries, and missing or inaccurate input data (Shoarinezhad et al. 2020).Some of these uncertainties are related to parameters used to simulate different processes in each model.In this study, the variation in the simulation results (target variables) due to perturbations in selected model parameters is analyzed.The main objective is to achieve an approximation of the parameter's uncertainties by using a simplified method, which is more economic in terms of computing time compared to other stochastic methods (e.g., Monte-Carlo simulations).Hence, the analysis of the variations in the target variable is performed with a First-Order Second-Moment Method (FOSM) (Gelleszun et al. 2017).The FOSM method, which was successfully validated by Gelleszun et al. (2017), is based on the variance-covariance propagation and, according to Kunstmann et al. (2002), the results are comparable to the ones obtained by applying more sophisticated methods (such as Monte-Carlo methods).
The covariance matrix of the selected target variable y is expressed by the following Eq.( 2): where C yy is the covariance matrix of the calculated tar- get variable y, with size m × m ; C xx is the empirical covari- ance matrix of the selected parameters, with size n × n ; A is the Jacobian, sensitivity or functional matrix, with size m × n and contains the partial derivations of the model with respect to its parameters; m is the number of time steps and n is the number of parameters.
The variance of the target variable y can be obtained from the diagonal of the covariance matrix C yy , according to Eq. ( 3): where a ij are the elements of the Jacobian matrix A and c jk are the elements of the empirical covariance matrix of the parameters C xx .
The variance-covariance propagation (Eqs.( 2) and ( 3)) gives the confidence intervals of the model with respect to the perturbations of the selected parameters.The empirical covariance matrix of the parameters, C xx , can be determined with the following Eq.( 4): where S e 2 is the empirical residual variance (scalar value) that can be obtained for the entire simulation period according to Eq. ( 5): y obs is the observed data (of the target variable y); y sim is the simulated data (of the target variable y); u is the length of the ( 2) available observed and simulated data and n is the number of selected parameters.The Jacobian matrix A is calculated by numerical deriva- tion in the optimum (central differences as an approximation of the derivatives that cannot be determined analytically in the case of numerical models).Then, each parameter is changed ± 1% from its optimum value (Hill 1998), according to Eq. ( 6): where i is the selected parameter; y _id is the target variable y obtained as an approximation of the derivatives, which composes the Jacobian Matrix A; opt_value_i is the optimum value of parameter i, obtained after calibration of the model; y _il is the target variable y obtained with the lower value of parameter i ( opt_value_i -1%) and y _iu is the target variable y obtained with the upper value of parameter i ( opt_value_i + 1%).
The empirical variance S e 2 (Eq.( 5)) gives an idea of the parameter perturbations in relation to the chosen optimization algorithm, which is used during the calibration of the model, to obtain the set of parameters that best simulate the target (output) variable y in each model.The empirical standard deviation, which has the same units as the target variable y, can be obtained as S e = √ S e 2 .
Finally, the approximate uncertainties of the model parameters can be represented with the standard deviation of the target variable, expressed by the root square of the variance (Eq.( 7)):

Target variables and selected model parameters
A target variable (output) is selected for each model.Furthermore, we choose a maximum of five parameters per model for the analysis of uncertainties to constraint computing times.Table 3 shows the target variable (output) for each of the models, whereas Table 4 summarizes the selected model parameters and their corresponding optimum values, which were obtained from the calibration processes for the single models.In addition, the lower value refers to the perturbation when the parameters have been decreased by -1%, ( 6) while the upper value refers to the perturbation when the value has been increased by + 1% from the optimum value.
For spatially distributed parameters, such as the C factor, or seasonal factors, such as the R factor, the respective mean values are given in Table 4.In all cases, the selected model parameters are the most sensitive ones and have the greatest impact on the simulation results in each model.A perfect agreement between observed and simulated bed levels in the reservoir is not to be expected since the WaSiM and RUSLE-SEDD models were calibrated for the Kokel gauging station and not for the reservoir (see Fig. 1).Consequently, the deviations in reservoir bed elevation may be closely related to under-or overestimation of the runoff and sediment load entering the reservoir, since they were not measured directly at the reservoir inflow.

Workflow
In total, 11 runs were carried out with WaSiM, 17 runs with RUSLE-SEDD, and finally 23 runs with SSIIM 2 to capture the changes in the parameters.Figure 3 shows the selected workflow applied to the modeling chain.

Model simulations under different climate projections
The modeling chain is used to predict the catchment's response under different future climate projections.Table 5 summarizes the GCM and RCM model combinations used under different Representative Concentration Pathways (RCP).These datasets are provided with a spatial resolution of 0.11 degrees (EUR-11 grid, WCRP 2009) and with a 3-hourly time step.As reference data, the ERA5 reanalysis dataset is used, considering a reference period from 01/1981 to 12/2010.The bias adjustment is performed according to the MultI-scale bias AdjuStment (MidAS) tool, v0.2.1, which provides cascade adjustments in time and space, using a day-of-year scaling step (Berg et al. 2022).
RCPs represent climate projections under different greenhouse gas concentrations that might lead to an increase in radiate forcing by the end of the century.For example, RCP2.6 is the lowest of all RCPs and expects a radiative forcing of 2.6 W m − 2 by 2100.Furthermore, each RCP is related to a global mean temperature increase compared to a reference period considered from 1986 to 2005.In the Bed elevation m a.s.l.Simulated bed elevation along the thalweg of the Devoll river within the reservoir case of RCP2.6, the global mean temperature increase is 1.0 °C, for RCP4.5 1.8 °C, and for RCP8.5 3.7 °C (Collins et al. 2013).
The main objective of our study is to scrutinize whether the spread of the climate projections for future simulations is broader than the variations in simulation results as a consequence of perturbations in the model parameters.In other words, we aim to analyze whether the uncertainties inherent to model parameterization are higher than uncertainties from climate projections.The perturbations in the model parameters are represented by a ± 1% change of their optimum value, according to the method previously described.For the comparison with the spread of the climate projections, RCP2.6 is selected.In this way, a tangible comparison (1% parameter variation vs. 1.0 °C increase in temperature) is carried out.Furthermore, we would like to highlight that in this case, the smallest climate change signals are generated and thus the results of this RCP are used for comparison to the parameter perturbation approach.By doing so, we avoid overestimating climate change signals.If the spread of the climate projections exceeds the approximate parameter uncertainties in the low emission scenario, this is also to be expected for the high and medium emission scenarios.

Results
The bed level changes along the upper part of the thalweg (> 5 km distance to the dam) of the Banja reservoir are presented, considering uncertainties in selected model parameters, but also different climate projections.The results of the climate impact simulations refer not only to the final target variable of the modeling chain (bed elevation along the thalweg of the Banja reservoir) but also to intermediate results (climate variables, monthly runoff, and sediment yield).Finally, a comparison between the uncertainties of the model parameters and the climate projections is performed.Figure 4.

Uncertainty quantification of model parameters
To assess the uncertainties associated with model parameters and their impact on the simulation results, the final target variable of the modeling chain is analyzed.Figure 3 shows the measured bed elevation along the thalweg of the Banja reservoir in 2016 and 2019, together with the simulated bed elevation in 2019.In addition, the standard deviation (dark gray shaded area) is shown, which considers the spreading of the simulation results from a total of 23 model runs, and may hence be related to the approximated uncertainties of the model parameters.The figure also includes the standard deviation obtained only for the reservoir model and without considering the variations of parameters from the previous models (yellow shaded area), thus focusing only on the parameters of the reservoir model.In both cases, the values refer to the standard deviation, which was obtained after applying Eq. ( 7).
The average value of the standard deviation for the reservoir model only (after Eq. ( 7)) is 0.28 m.The average standard deviation for the entire modeling chain, considering the uncertainties of all 11 parameters, is 0.64 m.Hence, it becomes obvious, that the largest uncertainties in the modeling chain result from the reservoir model.Figure 3 indicates in addition that higher variations of the target variable are located near the head of the reservoir, where a delta formation is visible.

Model simulations under different climate projections
Since precipitation and temperature are important forcing variables for the generation of runoff, soil erosion, and the consequent transport of sediments into the reservoir, a special focus is set on future changes in these variables.Figure 5a shows the decadal changes in the mean monthly precipitation regarding the reference period .The values are taken as an average value of the entire catchment and refer to the ensemble mean of the three GCM/RCMs and for RCP2.6.A positive change indicates an increase in the mean monthly precipitation values (green color), whereas a negative change indicates a decrease in the values (red color).Although there is no clear trend in the changes between the decades (y-axis), a reduction in precipitation during the summer months is visible (x-axis), whereas at the same time an increase during the winter months will occur.
Figure 5b shows the mean monthly changes in temperature, where a positive change suggests that temperature will increase in the future (red color) and a negative change suggests a decrease (blue color).The figure makes it visible that the mean monthly temperature will face an increase in the future, reaching higher values, especially during the spring months (April-May).
Finally, changes in mean monthly runoff at Kokel are analyzed (Fig. 5c).In this figure, a positive change suggests that the runoff will increase in the future (purple color) and a negative change suggests a decrease (brown color).In this case, there is a clear trend in the decreasing mean monthly runoff during the spring months (April-May), becoming larger by the end of 2050 and 2090.The reduction of the mean monthly runoff during spring is related to the rise in mean monthly temperatures, which leads to an increment in the evapotranspiration values and therefore less water will be available as runoff.Furthermore, the early melting of snow (shifted to late winter months) and the decrease in snow storage contribute to the reduction of runoff during spring.On the contrary, an increase in the mean monthly runoff is predicted mostly for the winter months.This is also due to the rise in mean temperatures, leading to early snowmelt, less snowfall and more rainfall.

Uncertainty assessment: parameter vs. climate model uncertainties
In this section, the uncertainties in model parameters and their impacts on simulation results are compared with the spread of climate projections for RCP2.6.In addition to the final target variable (i.e.bed elevation along the thalweg of the reservoir), results are also shown for the output variables of each individual model in the chain.

Mean monthly runoff at Kokel
The diagrams on the left in Fig. 6a show the mean monthly runoff at Kokel obtained from the water balance model, for three different periods as an ensemble mean of the three GCM/RCMs and RCP2.6 (rows, from bottom to top): 2011-2040, 2041-2070, and 2071-2100 and its corresponding standard deviations.In addition, the mean monthly runoff for the reference period (1981-2010, black dashed line) together with the mean monthly standard deviation due to uncertainties in model parameters (gray shaded area) are presented.
The hydrographs plotted on the bottom row for the period 2011-2040 indicate that the mean monthly runoff will not change significantly in the near future.However, a slight increase is expected during winter (an increase of 9.4 mm for January).For the rest of the year, the values will be on a similar level.
When looking at the second half of the century (middle row, period 2041-2070), there is a clear shift in the maximum value from spring (April) to late winter (March).Furthermore, the peak is below the values of the reference period (a decrease of 7.2 mm is expected for April).Similar results are obtained for the last period (upper row, period: 2071-2100), where the shift in peak flow from April to March manifests itself and a decrease of almost 9.0 mm is expected for the peak runoff.
Observing Fig. 6a, it is possible to conclude that the approximate uncertainties arising from model parameters in the water balance model are by far smaller (almost 5 times) than the ones coming from the climate impact models (spread of climate projections, measured as the standard deviation of the ensemble mean).On average, these values rise from 3.8 mm month − 1 to 15.2 mm month − 1 (300% rise), 18.6 mm month − 1 (389% rise), and 19.3 mm month − 1 (407% rise) for the first, second, and third period, respectively.

Mean monthly sediment yield at Kokel
The diagrams on the right in Fig. 6b show the mean monthly sediment yield at Kokel, obtained from the soil erosion and sediment transport model, for three different periods as an ensemble mean for the three GCM/RCMs and RCP2.6 (rows, from bottom to top): 2011-2040, 2041-2070 and 2071-2100 and its corresponding standard deviation.Similarly to the mean monthly runoff, the black dashed line indicates the mean monthly sediment yield for the reference period .In addition, the mean monthly standard deviation regarding uncertainties from model parameters (gray shaded area) is shown in the figure.
The sediment yield behaves similarly to runoff.Consequently, the maximum values are expected from February to April (between the end of the winter season and the beginning of the spring months).In general, the mean maximum sediment yield will not experience great changes, except for the near future, where an increase of around 20,000 tons is expected for March, which corresponds to the predicted increase in runoff.
Similar to the mean monthly runoff, the values are not expected to change significantly during summer (low-flow season) because erosion is strongly correlated with precipitation and thus with runoff.Furthermore, the standard deviation of the sediment yield regarding perturbations in the model parameters (gray shaded area) is also smaller (almost 5 times) than the spread of climate projections, given by the standard deviation of the ensemble mean.In this case, the values increase on average from approx.13,600 tons month − 1 to 44,300 tons month − 1 (225% rise), 55,700 tons month − 1 (310% rise), and 86,900 tons month − 1 (539% rise) for the first, second, and third period, respectively.

Bed elevation along the thalweg of the Banja reservoir
Figure 7 shows the bed elevation along the thalweg of the upstream part of the Banja reservoir, simulated with the The simulated bed levels in the near future ( 2036) show an increase in bed elevation within the reservoir, especially in the upper part.Here a clear delta formation is visible (compare Morris and Fan (1998)).For the mid-term period, in the year 2066, on one hand, the bed elevation will increase, but also a delta progression becomes visible, which is in accordance with literature.The highest depositions are observed at approx.10,500 m distance from the dam.Although in Fig. 7 it is not possible to gain insight into the impact of seasonality on the evolution of the bed elevation, in general, a higher accumulation of sediments occurs during months with high runoff and higher sediment transport (e.g., Fan and Morris (1992)).
Finally, by the year 2100, the maximum deposition height will increase from approx.10 m in 2036 up to approx.30 m.There is no clear difference between the maximum bed levels, but the largest increase occurs at a distance of around 9000 m upstream of the dam, which also indicates that the delta migrates further into the reservoir, when comparing the location of the delta at the end of the mid-term period (10,500 m distance to the dam).Hence, the deposition regime moves further downstream, whereas a sediment balance between erosion and deposition is established in the upstream part (9000-14,000 m distance to the dam).
In the case of the bed elevation, the spread of the climate projections is determined by the maximum difference between the ensemble members, represented finally as an average value for all the x-locations (distance from the dam).Similar to the previous simulation results (mean monthly runoff and mean monthly sediment yield at Kokel), the impact on the bed elevation along the thalweg of the Banja reservoir (target variable) is smaller than the spread of climate projections due to perturbations in the model parameters.The average change in bed elevation (along the thalweg) considering uncertainties in the model parameters is 0.64 m, whereas the average changes due to uncertainties in the climate projections are 0.96 m, 2.33 m, and 2.38 m for the years 2013, 2066, and 2100, respectively.These values represent an increase of 50%, 264%, and 272% compared to the average change of the bed elevation due to uncertainties in the model parameters.

Discussion
A complex modeling chain composed of different individual models is adopted to predict the sedimentation processes in the Banja reservoir in this study.The main advantage of applying such a chain lies in the detailed and process-based representation of each intermediate process.The results of the modeling chain are satisfactory and can be used for predicting the sedimentation processes under future climate conditions.
Nevertheless, uncertainties cannot be neglected.In our study, we focus mainly on parameter uncertainties and compare them to those inherent in climate models.The approximate uncertainties related to model parameters are determined by using the method developed by Gelleszun et al. (2017).The modeling chain includes three representative models that are used to study the impact of 11 sensitive parameters on the bed elevation changes of the Banja reservoir.Within this study, these 11 parameters were changed in the range of ± 1%, resulting in 23 model runs of the final model.According to the mentioned study, this approximate method proves to be robust and efficient, thus reducing dramatically the computing times that other sophisticated methods (such as Monte-Carlo simulations) might require.
Even though the ± 1% uncertainty of selected model parameters proved to work well, it needs to be considered that the possible change is strongly parameter-dependent, which means that for some parameters this change may be a major change, whereas for others, it might be considered as only a minor change.Hence, future studies should focus on determining how much the model parameters can deviate from their optimal value, to ensure that the approximate uncertainties are smaller than the spread of climate projections.
Additionally, uncertainties that might arise from climate model predictions are analyzed and compared to the approximate uncertainties in model parameters.In this study, the climate predictions are presented for RCP2.6, corresponding to three different GCM/RCMs realizations.The analysis of the results reveals that the approximate uncertainties in model parameters of the water balance model is significantly smaller than the uncertainties coming from the different climate projections (Fig. 6a).The same can be concluded when analyzing the simulated mean monthly sediment yield (Fig. 6b) and the final bed elevation along the thalweg of the Banja reservoir (Fig. 7).The results agree with the ones found by Kingston et al. (2011) and Wagner et al. (2017).They studied different uncertainties influencing their model results.In both studies, the authors conclude that the uncertainties arising from model parameterization are remarkably smaller than the ones generated by the climate projections for the case that the models were calibrated and validated with existing data in a first step.Among all RCPs, RCP2.6 is seen as the lowest mitigation scenario.Although reaching RCP2.6 emission values by the end of the century may be technically feasible, urgent actions are required to achieve this.For example, reducing emissions rapidly during the first decades of the century and increasing the use of renewable energy sources, for which countries beyond the Organization for Economic and Co-operation Development (OECD) are also required to participate (van Vuuren et al. 2011).Thus, climate projections under RCP2.6 may be too optimistic and the need of contemplating other However, RCP2.6 was chosen here from a methodological perspective: we intend to scrutinize whether even the spread of results from a set of climate projections with low change signals exceeds typical variations in results achieved through perturbation in the model parameters.The climate change signal is defined as the absolute difference between the ensemble mean values obtained for the future period and the reference period (historical climate), respectively.For example, for monthly runoff in the last period of the twenty-first century (2071-2100), the climate signals are (for the ensemble mean) 3.25 mm month − 1 for RCP2.6 (i.e., the difference between the colored and the black dashed lines), whereas a value of 6.20 mm month − 1 is obtained when considering RCP8.5.This suggests that the climate signals in RCP2.6 are on average in the same order of magnitude as the approximate uncertainties related to model parameterization, where a mean value of 3.8 mm month − 1 was obtained.Indeed, on a seasonal level, climate change signals can still  2066, and 2100 (20, 50, and 84 years after impoundment of the dam).The ensemble means of the three GCM/RCMs for RCP2.6 are shown (dark colored solid line), together with the mean value of spread between ensemble members (colored shaded areas).The black dashed line indicates the bed elevation in the past (for the year of finalization of the dam construction, 2016), together with its standard deviation, obtained from the impacts on simulation results due to perturbations in model parameters be higher.This is especially evident in the spring months, where the climate signal for runoff is almost 2.5 times higher.These findings show that low climate signals might be masked by model parameter uncertainties.On a minor note, the + 1 °C increase in the global mean temperature in RCP2.6 corresponds to ± 1% changes in model parameters as a mere and thus tangible comparison of numbers.
When predicting the response of a variable in the future (such as the bed elevation along the thalweg of the Banja reservoir), the use of multiple scenarios or ensembles is recommended, in order to cover further possibilities on how the climate is predicted (Collins et al. 2013).According to Teutschbein and Seibert (2010), complex Ensembled Regional Climate Models (E-RCMs) considering more than one RCM and a range of GCMs and RCPs are useful for hydrological simulations.Although our study can be classified into the mentioned group of E-RCM, considering further GCMs and RCMs might be interesting to understand how the climate spread changes and influences the prediction of the hydrological variables.
It is also worth mentioning that in our study, we focus only on the approximate uncertainties related to model parameters.However, other uncertainties might arise when applying such a modeling chain, such as the selected calibration approach (e.g., manual or automatic approach), the errors in measured data (e.g., in runoff or sediment yield), or the meteorological forcing data used as input.Additionally, changes in land use may contribute to alterations in the hydrological response of the catchment and should also be considered when predicting a catchment's response in the future (Li and Fang 2016).

Conclusion
A complex modeling chain was set up to predict the bed elevation along the thalweg of the Banja reservoir in the Devoll River (Albania), by considering hydro-climatic changes, monthly runoff, and sediment yield coming from the catchment.Despite the challenge of using three different models for predicting the final target variable, we benefit from the main features and accuracy of three process-based state-of-the-art models.As each model predicts a target variable, which serves as input for the subsequent model in the chain, the final target variable of the modeling chain depends strongly on the reliability of the antecedent results.To see how well this reliability can be ensured, model parameter uncertainties are studied for the entire modeling chain by using a simplified approach based on the FOSM Method.
These approximate parameter uncertainties (measured as a standard deviation) for predicting the bed elevation changes along the reservoir increased from 0.28 m (reservoir model parameters only) to 0.64 m when considering the uncertainties of all 11 parameters of the three models.Despite this increase, the values are still comparatively small (only 0.19% and 0.44% of the mean measured elevation bed in the year 2019, reprectively), and it can be concluded that the perturbations in the model parameters are not a significant source of uncertainty for the final simulation results.
Furthermore, three combinations of GCM/RCMs for RCP2.6 were selected to study the behavior of the involved variables in the future (until the year 2100).The spread of the climate projections is compared to the approximate uncertainties resulting from the perturbations in the model parameters.The results show that the spread of climate projections by the end of the century is on average larger than the approximate parameter uncertainties, being 408%, 539%, and 272% higher for the prediction of runoff, sediment yield, and bed elevation, respectively.However, as demonstrated in the case of runoff, they are in the same order of magnitude as the climate change signal inherent in the mitigation scenario RCP2.6.Nevertheless, and given that change signals are higher in other RCPs, the use of such a complex modeling chain is a valuable tool for predicting sedimentation processes in a reservoir for different climate change scenarios.The method described in this paper highlights how the parameter uncertainty for each model is quantified approximately, whilst demonstrating their robustness when comparing the larger spread imposed by climate projections.This is in particular helpful to guide modelers and practitioners to communicate different sources of uncertainties in complex modeling chains including climate models, and to highlight how uncertainties compare to climate change signals.

Fig. 1 a
Fig. 1 a Location of the study area within Albania; b important tributaries and outflow boundaries of the Banja reservoir; c topography, subcatchments, gauging stations and location of the Banja reservoir.(adapted from the European Environment Agency (2016))

Fig. 2
Fig. 2 Selected modeling chain and target variables for the study case of the Devoll catchment and the Banja reservoir

Fig. 3
Fig. 3 Selected workflow applied for the modeling chain of the Banja reservoir located in the Devoll catchment (Albania).The three models, their parameters, perturbations and the number of model runs are shown as well

Fig. 4
Fig. 4 Measured (solid red line for the year 2016, solid blue line for the year 2019) and simulated (dotted blue line for the year 2019) bed elevation along the thalweg of the upstream part of the Banja reservoir as a result of executing the entire modeling chain.The standard

Fig. 5
Fig. 5 Decadal changes of mean monthly values relative to the reference period (1981-2010): a change of the mean monthly precipitation as an average for the entire catchment [mm month −1 ] ; b change of the mean monthly temperature as an average for the entire catchment [°C month −1 ]; c change of the mean monthly runoff at Kokel [mm month −1 ].All values refer to the ensemble mean of the three GCM/RCMs and RCP2.6

Fig. 6
Fig. 6 Mean monthly runoff a and mean monthly sediment yield b at Kokel station for three different periods (rows, from bottom to top: 2011-2040, 2041-2070, 2071-2100).The values represent the ensemble mean of the three GCM/RCMs and RCP2.6 (dark colored solid lines ± standard deviation of the ensemble mean).The black

Fig. 7
Fig. 7 Bed elevation along the thalweg of the Banja reservoir [m a.s.l.] for the years 2036,2066, and 2100 (20, 50, and  84  years after impoundment of the dam).The ensemble means of the three GCM/RCMs for RCP2.6 are shown (dark colored solid line), together with the mean value of spread between ensemble members (colored shaded areas).The black dashed line indicates the bed elevation in the past (for the year of finalization of the dam construction, 2016), together with its standard deviation, obtained from the impacts on simulation results due to perturbations in model parameters

Table 2
Summary of the soil loss erosion risk factors

Table 4
Summary of selected model parameters for each model, including definition and units *All sub-processes belong to the unsaturated zone model.The optimum value refers to the value obtained from the calibration and lower and upper values refer to the ± 1% variation of the optimum value