Long-Term Water Management Model for Preserving Sustainable Useful Capacity of Reservoirs

In this study, we present a long-term water management model which can be an important tool for the ordinary management of a reservoir. The model can be made simple and applicable in a graphic form. The key point of the model is that the sedimentation rate of solids can be modeled according to a power law. The basic equation of the model is tested using bathymetric surveys of different dams located in different parts of the world. Robustness and predictive power are evaluated both numerically and by comparison with satellite results. Then the model is transformed into the form of a management plot and applied to the real case of the Camastra dam.


Introduction
The dams are built for drinking, irrigation, industrial purposes as well as for electricity generation and flood protection (Khagram 2018). However, due to the sedimentation phenomenon, their useful capacity decreases over time (Covelli et al. 2020). So the minimum capacity necessary to meet the water demand, fails.
Dams modify sediment transport in the upstream and downstream riverbeds (Mirauda et al. 2018), resulting in loss of reservoir storage and reduced usable life, and depriving downstream reaches of sediment essential for channel form and aquatic habitats (Cimorelli et al. 2021). The reduction of useful capacity has a direct effect on water sustainability, i.e. on the ability of individual countries to meet the water demand and to cope with climate change. The correct development of a reservoir management plan ensures that water reserves are maintained for a long time (Bufalo et al. 2018). In order to draw up a management plan for a reservoir, which is capable of combining the technical term of useful capacity and the social term of water sustainability, herein we introduce the key-concept of sustainable useful capacity of the reservoir.
The basic parameter for monitoring the capacity of a reservoir at a given time is the volume of sediment deposited up to that precise time. Since this quantity can be estimated either from the amount of solids in the dam or from soil erosion, two predictive approaches have been developed. The first of these considers the sediment thickness to be trappingefficiency result (Brown 1943;Churchill 1948;Brune 1953;Chen 1975;Heinemarm 1981;Campos 2001;Espinosa-Villegas and Schnoor 2009;Xiaoqing and Yang 2003;Castro and Mantilla 2021). It is employed to plan extraordinary management operations.
The second approach estimates the sediment input to resevoir by Revised Universal Soil Loss Equation (RUSLE) (Renard et al. 1991;Kim and Julien 2006;Ranzi et al. 2012;Tanyaş et al. 2015;Gaubi et al. 2017). It is used is for estimating in soil losses under changes in climate or land use.
Both approaches require a large computational effort but, thanks to modern electronic computers, the effort does not by itself present insuperable difficulty. The major disvantages of these approached are not computational; rather they are due to the fact that we do not have a equation satisfactorily applicable to reservoirs over the full variability-range of parameters.
By describing a difficult, real problem it is sometimes possible to obtain a simple solution to the problem not in terms of immediate physical reality, but in terms of mathematical quantities which are suggested by an abstract description of the real problem . Our concern in this paper is concentrated exclusively on long-term water management of a reservoir in terms of its sustainable useful capacity. A model where the lifetime of the sedimentation process follows a power law is derived and discussed. The theoretical parameters of the model are related to different operating conditions of dams working in different locations in the world. Finally, the robustness and predictive capacity of the method were evaluated by comparing our numerical results with satellite surveys.
One of the main advantages of the method is that few bathymetric results, even approximate ones, can be converted into information for the management project. Based on bathymetric measurements, it is quickly usable and does not require additional information, therefore it can be systematically extended to dams operating under any conditions.

Model Details
The physical and chemical characteristics of solid sediments make the sedimentation process extremely complex. The spatial distribution of sediments as well as water management procedure can influence the settling process in the reservoir. Herein the sedimentation rate is modeled with a lifetime variable according to a power-law.

Sedimentation Volume
Sedimentation is the process under which a mixture of fluid and solid particles is separated into in solid and liquid phases under the gravity influence. In this separation process, the fluid flows upwards while the solids move downwards and within the suspension the volume fraction of solids, , is a function of the position z and time t.
The continuity equation for the solid phase, considering only uni-dimensional movement Tiller and Khatib (1984), is where (−u s ) is the positive value of sedimentation velocity of solids ( u s is negative because the motion direction is contrary to the positive direction of axis z).
Let [H(t)] be the boundary separating solid and liquids phases, then instantaneous average of Eq. 1 gives where At the solid-liquid interface, the falling speed of the solids is zero so that we can rewrite Eq. 2 as Moreover, the concentration of solids in the suspension is related to the useful volume by a simple mass balance where V is the falling volume of the solid particles and C s their mass-concentration (kg ⋅ m −3 ). If the partial specific volume of solids, v s , is considered to be time-independent, Eq. 5 is transformed into hence combining Eq. 6 with Eq. 4 one obtains where, for the sake of simplicity, we have set (0, t) = 0 (t) and u s (0, t) = u s0 (t) . The falling volume can be expressed as the difference between the total capacity, V ∞ , and the volume of solids, V s , sedimented at time t, i.e.
We specify that V ∞ is the volume that can be filled with sediment in an infinite time, namely the original capacity of the reservoir. Therefore, Eq. (8) becomes ( Direct integration of Eq. (9) yields the key equation of the model Obviously, such an equation has a operational validity only if the analytic expression of the (t) function is known. We assume a power law of the type where 0 is the characteristic time scale of the reservoir filling and a parameter which measures the level of complexity of the settling process.
It has been known for some time that power law models are good (albeit empirical) descriptors of the kinetics of complex systems. Power-law kinetic confers greater flexibility to the form of the rate functions than classical kinetics. It can be argued that power-law kinetics provide approximations to kinetics in general (Stroud et al. 2006).
Under these conditions, Eq. (10) can be easily integrated to yield At a glance one might be tempted to conclude that the model (Eq. (12)) requires three adjustable parameters. A second look shows that this conclusion is not justified because V ∞ is a constructive characteristic of a dam, hence it is a known parameter. Accordingly, only the parameters ( , 0 ) have to be determined to monitor the sedimentation process in a reservoir. Nevertheless, the calculation procedure depends on the number of bathymetric surveys available to the dam manager. Moderate number of surveys Due to the difficulties in bathymetric measurements, the number of data available to the dam manager is often insufficient to implement a nonlinear least squares procedure. In these situations a simple and efficient estimation procedure can be implemented to evaluate the model parameters. For this purpose it is convenient to introduce a new variable by so that Eq. 12 is transformed into Thus, the straightness in the plot of log Z versus ln t , provides an estimate of model reliability.

Large number of surveys
Approximate estimate from only two surveys For recently built dams only two reliefs are available so that the direct application of the regression procedures loses significance. In these cases the dam manager, may estimate the parameters ( , 0 ) by taking advantage of the linearity of Eq. (14). Indeed, let t 1 and t 2 be the times at which the bathymetric measurements V s1 and V s2 were carried out. Hence, results

Sustainable Useful Capacity of a Water Reservoir
In the absence of extraordinary events, the useful capacity of the reservoir, after t years, is given by While the concept of available useful capacity is intuitive, the definition of sustainable useful capacity is less obvious due to the difficulty of identifying the future functional period of the reservoir. The consumption of water represents the needs for drinking, civil and industrial use. Generally, the period of maximum demand for water is one semester per year. Consequently, we define semiannual consumption as where 0.5 indicates that the semester begins six months after the start of water storage (half a year).
It is obvious that the useful capacity will satisfy the users if and only if on the contrary, the reservoir will be in serious risk if To define sustainability, we assume as a basic value the maximum of the S function for all the years of operation of the reservoir, i.e.

Results and Discussion
In order to test the model in calculating the long-term reduction in the useful capacity of a reservoir, the volume of accumulated sediment was calculated for several dams distributed in different parts of the world. The model was applied to calculate the useful capacity and sustainability of different dams distributed in different countries of the world. Since the model reconstructs the sedimentation profile on the basis of bathymetric surveys, the choice of dams was made in terms of available bathymetric data. Bathymetric data is the only information needed to run the model so for each dam investigated we only present the geographic location, initial capacity and year of commissioning. Therefore, the dams have been divided on the basis of available bathymetric surveys, i.e. according to the procedure used to calculate parameters ( , 0 ).

Nonlinear Procedure
For these reservoirs a relatively large number of bathymetric surveys are available and a non-linear fit procedure can be employed. The quality of the fittings was assessed using a non-linear fitting correlation coefficient R 2 . The dams that meet this requirement are Welbedacht and Pong.
Welbedacht is a gravity dam located in the South Africa. It was completed in 1973 with original capacity V ∞ = 1.14 ⋅ Therefore, if one has a sufficient number of bathymetric data, the model provides a tool capable of establishing the amount of sediments and evaluating how to manage the water in the reservoir.
Camastra dam is a barrier of the Camastra river, a tributary of the Basento river, in Potenza (Italy). Its construction was completed in 1967 with an original capacity V ∞ = 35.3 ⋅ 10 6 m 3 . This dam deserves a particular comment, in as much as, although the number of bathymetric reliefs is not high, the available data must be processed with the non-linear procedure.  (India). For these dams a suitable number of bathymetric reliefs is available so that the model was tested with a non-linear fit. The dots are the results of the bathymetric surveys and the red curves the fitting to the model. c Camastra dam (Italy). The red curve is the fitting to the available bathymetric reliefs. For the blue curve in the fitting a value of estimated under "anomalous"conditions (circle) is also used. The closeness of the curves indicates the "robustness"of the method As Fig. 1c reveals, for this dam the first bathymetric survey was carried out after 26 years of the dam's operation. On the other hand, to apply the linearization procedure, used in Eq. (14), the point (0.0) has to be excluded from the analyzes. Here, this would involve the calculation of the slope of the line from points very far from the origin of the axes, with large errors in calculating and 0 . On the contrary, in the nonlinear fitting of bathymetric data, the points (0,0) and (0, V ∞ ) are fixed and the curve V s (t) has fewer degrees of freedom. This robustness of the non-linear method allows it to produce good results although the number of experimental data is not high enough. In this case, the value provided are = 0.4 ± 0.1 and 0 = (1.3 ± 0.6) ⋅ 10 2 y and t 1 2 = (5 ± 3) ⋅ 10 y and the calculated sedimentation profile covers the bathymetric results very well (see Fig. 1c).
Results are in agreement with the technical finding that after 50 years of operation the useful capacity of the reservoir has been reduced by 50% (Cimorelli et al. 2021). The erosion phenomena, which developed on the slopes overlooking the lake, contributed to speeding up the sedimentation in the first 20 years of the reservoir's life. In the following years, the reduction of the deposit of sediments was due both to the raising of the bottom of the river section upstream of the reservoir inlet, and to the effectiveness of the river network works carried out before and after the dam construction. In order to further investigate the method robustness, we decided to repeat the non-linear fitting by introducing another bathymetric survey, carried out after 21 years of operation in anomalous conditions and therefore with reduced reliability (the empty circle in Fig. 1c)). As one can see, the two curves differ very little, indicating that the points (0, 0) and (0, V ∞ ) strongly condition the calculated profile.
Gobind Sagar dam, of an original capacity V ∞ = 9.87 ⋅ 10 9 m 3 , dams the Satluj river (India) and went into operation in 1963 (Shukla et al. 2017). For this dam a relatively large number of bathymetric surveys are available to use the non-linear procedure, however, it cannot be applied as it is. Indeed, if the model equation is fitted to experimental data it produces = 2.23 ± 0.03 and 0 = 34.3 ± 0.2 . It is important to point out that > 1 is physical nonsense as it would mean that the lifetime (t) decreases over time. That is to say small particles settle first and then the larger ones. Thus, the result > 1 is a strong signal of some anomalous condition. On the other hand, as one can see in Fig. 2a the first points deviate markedly from the model. This led us to hypothesize a delay in the recording of bathymetric measurements so that the model was generalized as follows where t 0 is the delay time.
Therefore, the robustness of the non-linear method is well suited to different real situations and, moreover, it is able to signal any delays in bathymetric measurements. Table 1 collects all the parameters determined with non-linear fitting under different real conditions.

Linear Procedure
For several dams the number of bathymetric surveys available is not so numerous as to make the non-linear procedure reliable. In these cases one can take advantage of the linearization of the model. In other words, Eq. (14) is fitted with a linear regression procedure and the model parameters estimated given regression coefficients.
Panchet dam was built in Dhanbad district's Panchet (India) and completed in 1959, with an orginal capacity V ∞ = 1.67 ⋅ 10 9 m 3 . The plot of log Z vs. log t displayed in Fig. 3a exhibits a good straightness. From this straight line the parameters = 0.42 ± 0.04 and 0 (4.1 ± 0.1) ⋅ 10 3 are extracted with correlation coefficient R 2 = 0.98 . The fitting parameters were then used to recalculate the sedimentation profiles. Results displayed in Fig. 3a' indicate a good predictivity of the experimental curve.
Rocky Gorge dam is located on the Patuxent River in Howard County, Montgomery County, Maryland (USA). It was completed in 1952 by a original capacity of 24.2 ⋅ 10 6 m 3 (Ortt et al. 2007). For this reservoir, the linear plot shown in Fig. 3b provides = 1.1 ± 0.3 and 0 = (2.7 ± 0.7) ⋅ 10 2 y with a correlation coefficient R 2 = 0.99 . The sedimentation profile recalculated by these parameters carefully reproduces the bathymetric reliefs (Fig. 3b').
Triadelphia dam is located on the Patuxent River, whitin the Humpstead and Glenwood district (Maryland, USA). The Triadelphia dam was completed in 1942 with original capacity V ∞ = 2.727 ⋅ 10 7 m 3 . The linear regression gives = 1.2 ± 0.5 and 0 = (3.1 ± 0.9) ⋅ 10 2 y with a correlation coefficient R 2 = 0.97 (see Fig. 3c). The sedimentation profile calculated with the linear regression parameters, shown in Fig. 3c', indicates that the relationship accurately predicts the experimental profile. Fig. 2 Gobind Sagar dam a Non-linear fitting of bathymetric data leads to a concavity-up, devoid of engineering significance; b Assuming a delay t 0 , the model fits experimental data with the right concavity and a correlation coefficient R 2 = 0.999

Estimation of Parameters from Two Bathymetric Measurements
In some real situations, there may be management problems such that it is not possible to carry out bathymetric surveys with regularity. In these cases we have only two experimental values such as, for example, for Chimahanda and Altinapa dams.
Chimhanda dam is located on the confluence of Runwa and Mwera rivers (Africa). It was completed in 1988 with original capacity 5.2 ⋅ 10 6 m 3 (Dalu et al. 2013). Bathymetric data for this catchment are limited and the non-linear fitting is unreliable. Based on Eq. 15, one finds = 0.32 and 0 = 260 y. The sedimentation profile recalculate with these parameters is shown in Fig. 4a. These estimated parameters predict t 1∕2 = 82 y, that is to say that in 2050-2070 half of the reservoir is filled with sediments. However, in 2015, or after 27 years of dam operation, 35% of the reservoir was silted up. Although the forecast may seem distant, one must take into account that the dramatic phenomenon of deforestation currently gripping Zimbabwe to obtain agricultural land and mining activities, has eroded the basin upstream of the Chimhanda reservoir. Sedimentation processes have been accelerate abnormally (Tundu et al. 2018).
Altinapa dam is located in Konya (Turkey) and it was completed in 1967 with an original capacity V ∞ = 3.76 ⋅ 10 7 m 3 (Ceylan et al. 2011). For this dam parameter Eq. (15) provides = 0.77 and 0 = 262 y. The recalculated sedimentation profile shown in Fig. 4b exhibits t 1 2 = 162 y. If the parameters are used to calculate the sediment volume after 42 years of dam operation, i.e. in 2009, one finds that only 21% of the initial capacity is lost. Fig. 3 a- This value can be explained keeping in mind that Altinapa is a flood control reservoir requiring frequent bottom outlets openings with partial loss of sediments.

Robustness and Predictive Capacity of the Method
In order to evaluate the robustness of the method and the predictive ability, both numerical analyzes and comparison with satellite surveys were performed.

Numerical Analysis
The numerical analysis on bathymetric data consisted in fitting the model to an ever smaller number of data points. In other words, the experimental points were eliminated one at a time and the remaining measurements were always fitted with an exponential stretch. As an example, Fig. 5a illustrates the case of the Welbedacht dam, where one can note that while eliminating 89% of the points the model still predicts the measures a good accuracy. Indeed, if the fitting is carried out on experimental data collected up to 1978 (three experimental data points) and the volume of sediment projected in 2018 is calculated, one finds V s = 11 ± 1 against the measured value of 11.1 (De Villiers and Basson 2007). A difference of about 1% indicates an effective robustness of the model. This is because, regardless of the number of experiments performed, the fitting equation fixes the starting point (start of operations) and the point of maximum capacity of the reservoir (end of operations) which condition the shape and curvature of the sedimentation. As a further example in Fig. 5b the forecast on the Camastra dam is displayed. The experimental data available up to 2015 were fitted with the nonlinear model. The 2017 experimental value of 17.7 ⋅ 10 6 m 3 was compared with the predicted one of 17 ± 0.5. Results of the numerical analysis for all investigated dams are collected in Table 2.

Comparison with Satellite Surveys
Remote Sensing has become widespread as a tool for investigating the status and quality of waters. This method exploits the interaction of water bodies and sunlight along the electromagnetic spectrum to reconstruct the characteristics and content of the optically active substances present in the investigated water column. Thanks to its ability to collect data frequently from instruments aboard satellites, Remote Sensing allows monitoring of the concentrations of suspended solid sediments (Tsolakidis and Vafiadis 2019). In order to compare results of the two techniques, the sedimentation profiles were calculated from the bathymetric data and, then, the predictions compared with values of satellite surveys. In order to test the predictive power of the model, the calculated sediment volume was compared with that measured by satellite surveys for the Vaigai and Ramganga dams (Merina et al. 2016).
Vaigai dam collects the waters of the river Vaigai and it was completed in 1958 with an original capacity V ∞ = 1.95 ⋅ 10 8 m 3 (Khan and Nazar 2015). Sedimentation profiles, shown in Fig. 5c, were obtained from linear regression on the bathymetric surveys for years 1958, 1976, 1981, 1983 and 2000. The extracted parameters = 0.58 and 0 = 1003 allow to determine the sediment volume in 2012 V s = 3.3 ⋅ 10 7 m 3 , while the value established by Merina et al. with satellite measurements it turns out to be 3.2 ⋅ 10 7 m 3 (Merina et al. 2016) and is shown in the plot of Fig. 5c.
Ramganga dam was built on the Ramganga river in 1975 with an original capacity V ∞ = 2.5 ⋅ 10 9 m 3 (Jain et al. 2005). The sedimentation profile was calculated by bathymetric data available for the years 1974, 1988, 1997. The parameters = 0.73 , and 0 = 1522 y were obtained from linear fitting, hence the sedimentazion profile displayed in Fig. 5d was calculated. The sediment volume estimated in 2000 by a satellite survey . This point is in agreement with the trend of the sedimentation profile, indicating a good predictive ability of the model.

Application to Ordinary Management of a Dam
An important question is how these results are related to ordinary management of a reservoir. This involves controlling the evolution of the dam and determining the time at which it enters the risk zone. That is to say the time when the reserve is no longer able to meet  the demands of civil society. Although the discussion can be kept general and extended to any type of dam, for the sake of simplicity we will refer to the case of Camastra dam whose model parameters are collected in Table 1. At this point it should be noted that the semester of maximum consumption in Italy runs from May 1st to October 31st. Then, having established the parameter ( 0 and ), Eq. (16) may be used to deduce the reservoir capacity W(t − t 0 + 0.5) , however results are more effectively illustrated graphically rather than with fitting equations (management plot).
As shown in Fig. 6, this capacity is a function that decreases over time (black curve). The maximum possible consumption is a horizontal line in the same plot (red line). The two curves intersect at point P at time t B . It is evident from that for for t < t B is W(t − t 0 + 0.5 > S max so that the reservoir is able to satisfy all the required needs. On the contrary, for t > t B is W(t − t 0 + 0.5 < S max and the reservoir does not have enough water to meet global needs. At the time t = t B , the dam manager has to make important decisions in order to get out of the risk of closure. The manager can make the decision D 1 to reduce consumption by eliminating some users. In this case the S max curve is lowered (red dashed line) and the intersection point moves to t ′ B . The life of the reservoir has increased by t � B − t B years, although with reduced consumption. However, the D 2 decision can be made. The segment P Q in Fig. 6 represents the amount of sediment deposited at time t B , therefore the decision D 2 is equivalent to the removal of the volume P D 2 from the reservoir. This is equivalent to a volume of available water corresponding to point U, at time t R . In other words, the removal of P D 2 amount of sediment returns the reservoir to the condition of time t R , thus its life has increased by t B − t R years.

Conclusions
A mathematical model, describing the sedimentation rate as a power law, is proposed. It uses only two external parameters to be calculated with a fitting procedure. The model calculates the sedimentation profiles by fixing the starting and ending points of the curve, this makes it very robust and effectively predictive. The equations were tested using bathymetric surveys from several dams located in different parts of the world. To validate the model, both a numerical approach and a comparison with satellite surveys were used. Finally, it is shown how to render the model in graphic form to make it accessible to dam managers and to keep water sustainability under constant control. The large variability of the geographical contexts of the dams investigated assure us of the model reliabiliby in a wide range of real situations. This wide applicability is the manifestation of the fact that the parameters and 0 eccompass the specific dam-reservoir features. The model provides a framework which enables us to correlate data in sensible manner, it tells us what to plot against what, the coordinates we must use to get a smooth line. For enginering work such a framework is extremely useful because it enables to interpolate and extrapolate limited bathimetric reliefs and to manage reservoir's water in a sustainable way.
Author Contributions All authors (BM, AD, AM and LA) contributed to the conception and design of the study, the analysis of data and the writing of the paper. The model was designed and managed by BM and LA. All authors read and approved the final manuscript.
Funding No funding was received for conducting this study.

Availability of Data and Materials
The data used in this research are taken from the literature. They are the official data of dam managers in different countries of the world.
Code Availability Custom code based on MATLAB R2021a has been used.

Declarations
Ethics Approval All work is in compliance with Ethical Standards.

Consent to Participate
Authors consent to their participation in the entire review process.

Consent for Publication
The authors give their permission to publish.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.