Influence of reservoir shape upon the choice of Hydraulic vs. Hydrologic reservoir routing method

. The objective of the paper is to compare the hydrologic and hydraulic reservoir routing methods in term of assumptions, equations, numerical computation procedures, necessary data and advantages-disadvantages of their use. To test the results provided by the two methods, a set of two reservoirs from Romania was chosen: one long and narrow and the other one roundly shaped. Corresponding inflow hydrographs were chosen, and similar conditions were imposed for the outflow dam control structures, namely the initial water level in the reservoir to be at the Spillway Crest Level (SCL) and no outflow control. For the hydrologic method the Puls procedure was used and a program was written in Scilab to solve the continuity equation in finite differences. For the hydraulic method HEC-RAS software was used to solve the 1D Saint-Venant equations. Outflow and stage hydrographs at the dam were compared together with the stage hydrograph at the reservoir tail for the hydraulic method. Results show that the hydraulic method should be used for the long and narrow reservoirs, as it considers the backwater effect, whereas the hydrologic method can be efficiently used for all other reservoirs where this effect is negligible.


Introduction
Flood attenuation is one of the roles of storage, flood control or multipurpose reservoirs created by dam impoundment, as a key component of hydropower development schemes. Attenuation has a major protective impact upon downstream areas by decreasing the peak discharge and delaying the time of its arrival, therefore diminishing depth and extent of downstream flooding, protecting river banks of erosion and the floodplain of scour, at the same time limiting the transport and deposition of large quantities of sediments [1]. Reservoir (or storage) routing is a mathematical procedure for computing changes over time in magnitude and shape of a flood wave transition through a water retention facility. The wave travels from the reservoir's tail to the outflow structures of the dam. Computations are made considering the elevationstorage and elevation-discharge characteristics of the reservoir and of its outflow structures, respectively.
Applications are in the flood forecasting, regulation and protection, design of reservoir and its outflow structures and reservoir operation, since the flow downstream the reservoir is usually limited to a predetermined maximum value [2].
For computing the outflow hydrograph two reservoir routing methods are used: (i) hydrologic method and (ii) hydraulic method. The first conceptual method solves the continuity equation using a relationship between storage in the reservoir and the outflow discharge under the assumption of a horizontal water surface elevation (known as level pool reservoir routing). The second one, solves both the continuity equation and the momentum conservation equations (describing the physics of water movement), i.e. the nonlinear, hyperbolic, partial differential Saint Venant equations.
From mathematical point of view, the hydrologic methods involve simplified numerical techniques (usually finite difference) and steady-flow hydraulics [2,3,4]. For solving the hydraulic Saint-Venant equations, which take into account the backwater effect of the water surface profile along the reservoir, various numerical methods have been proposed, such as: direct finite difference (implicit and explicit) -considered as the standard for 1D approximation, method of characteristics, finite volume and, only in case of 2D and 3D models, finite element methods [5].
In terms of necessary data/input, both methods need the inflow hydrograph I(t). The hydrologic method however requires only a relationship between elevation (or stage) and storage (elevation-storage/capacity curve) and an elevation-discharge relationship for the outflow (usually a spillway), while the hydraulic method needs detailed topobathymetric data for the reservoir (continuous bottom surface elevations for the 2D approximation or cross-section profiles with distance between them for the 1D approximation), roughness coefficients, and the same elevation-discharge relationship for the outflow as a downstream boundary condition.
Therefore, if the hydraulic method is more data consuming and its equations are more difficult to integrate than the hydrologic one, a justified question pops up: for which reservoir is worth using the hydraulic routing? Fenton J.D., [6] and Haktanir, T. & Ozmen, H. [7] state that for long and narrow reservoirs where backwater effect of water surface profile is important the hydraulic routing should be used, while for the others there is no difference in applying the hydrologic routing. The objective of this paper is to prove this by analysing comparatively the two methods of reservoir routing for two differently-shaped reservoirs from Romania: one long and narrow and one roundly shaped.

Hydrologic method
Level pool routing is the procedure for calculating the outflow hydrograph from a reservoir with a horizontal water surface (figure 1), given its inflow hydrograph and storage characteristics [8]. The characteristics of the inflow hydrograph, depicted in figure 1 are: I(t) and O(t) input and output volumetric flow rates respectively (inflow and outflow), S(t) = water storage volume in the reservoir, T is the total duration of the flood wave, Ti and Tr rising and recession limb duration respectively, tlag the difference in time between the peak inflow and peak outflow (known as translation), and a, attenuation, or the flow difference between the two. Neglecting precipitation, evaporation, infiltration and tributary inflow, the continuity (mass conservation) equation is given by: (1) The reservoir storage may therefore be written as: where t1 is the intersection time of the two inflow and outflow hydrographs (figure 1). For this time, I(t1) = O(t1) and therefore water storage in the reservoir, S(t) has reached its maximum, / 0 dS dt = , meaning that the period of accumulation in storage has ended and the release/depletion period from storage begins. If at the initial and final time, t = 0, water level is at the dam spillway crest the two highlighted surface areas in figure  1 representing storage and outflow volumes have to be equal.
For two successive moments 1 and 2 corresponding to the beginning and the end of a time step, equation (1) written in finite differences becomes ( ) According to the Puls method [9] equation (3) may be rewritten in terms of storage volumes as follows with all the terms in the left-hand side being known at the beginning of each time step and all terms in the right-hand side being unknown at the end of the time step. Alternatives of the Puls method are: Modified Puls method, which postulates that storage depends only on outflow rate through a nonlinear function, S=f(O) [10], Goodrich method (also semi-graphical) and the Standard Fourth-Order Runge-Kutta Method (SRK) [2,11]. Since equation (4) has two unknowns, (storage and outflow at the end of the time step, S2 and O2), another relationship to relate S=f(O) is needed. This may be found if two other relationships are known: S(h) or the elevationstorage relation and O(h), or the elevation-discharge equation of the outflow hydraulic structure, where h is the reservoir water depth (alternatively, water surface elevation in the reservoir, z, may be used).
A unique (invariable) relationship between storage and outflow holds for reservoirs having horizontal water surface. Such reservoirs have a pool that is wide and deep compared with their length in the flow direction [8]. However, the invariable storage relationship S=f(O) requires that a fixed discharge from the reservoir is released for a given water surface elevation at the dam, which means that the outlet works must be uncontrolled, or, if controlled, the gates should be held at a fixed position [8].
Therefore, the data needed for Puls method applied for level-pool routing of uncontrolled reservoirs are: (i) S(h) = depth -storage relationship; this is developed from topographic data prior to building the dam and may be expressed through a table or a curve/plot. Usually, the elevation-storage curve may be approximated for values around the Normal Pool Level (NPL) elevation value by a power-law equation as follows: where k and m are coefficients determined experimentally; (ii) O(h) = elevation-discharge relationship that depends on the geometry of outflow structure, such as an uncontrolled overflow spillway (or a morning glory spillway, a gate-controlled spillway with constant opening): where Cd and C are the non-dimensional and the dimensional discharge coefficients (the latter including the acceleration due to gravity) and hd -the nappe height. Alternatively, for an orifice (bottom or lateral outlet) the outflow is expressed as: where A is the area of the orifice and hO is the depth of the orifice axis; (iii) I(t) = Inflow hydrograph of the flood wave is given by the rainfall event with its peak flow, Ip, having a certain exceeding probability, p(%). For the present study the inflow hydrographs were taken as being the synthetic ones to produce the largest allowable outflow in the river channel downstream the dam. The total flood volume is represented by the area beneath the flow hydrograph and computed with: The computation time step, t must be small enough so that the hydrographs can be assumed to be straight lines between two successive time values. As a rule of thumb, the time step has to be smaller than 0.2 of the time to peak of the inflow hydrograph [12]. Equations (4), (5) and (6) may be solved simultaneously semi-graphically or by inverse interpolation of tabular values [6]. With the advance of computerization graphical methods have been replaced by functional methods so that computational procedure can be automated [8].
For the present study, the Puls procedure was used and a program was written in Scilab to provide a numerical solution of the equation (4) and thus obtaining the outflow and stage hydrographs.

Hydraulic method
The hydraulic method consists in routing the inflow hydrograph between the tail of the reservoir and the outflow structures at the dam (similarly to the flood routing along river reaches). This is accomplished by applying the Saint Venant equations for 1D unsteady flow conditions embedded into the HEC-RAS software program [13]: In these equations, x is the distance along river channel/reservoir, h the depth of the reservoir in each cross-section, A the cross-sectional area, J the bed slope,  the kinetic energy coefficient, Jf the friction (energy) slope due to roughness and ql the lateral inflow into the reservoir from other tributaries (usually equal to 0).
Under unsteady flow conditions the Chezy-Manning relationship becomes biunivocal, thus leading to looped rating curves: where n is the Manning friction coefficient and P is the wetted perimeter. If all terms are considered in the reservoir routing process, the method is also called dynamic routing [14]. Necessary data to build a numerical model and solve these equations through the finite difference method with the help of HEC-RAS software (implemented with the 4-point implicit finite-difference numerical scheme) are: (i) the geometry of the reservoir, dam and outflow structures. The topobathymetry of the reservoir must be described by multiple cross-sections with known distance between them; (ii) upstream boundary condition consists of the inflow hydrograph, I(t), at the tail of the reservoir whereas as the downstream boundary condition is considered the elevation-discharge equation of the dam outflow structures, O(h) (equation (6) or (7)); (iii) as initial conditions a water surface profile for the initial value of the inflow hydrograph along the reservoir must be developed for steady gradually varied flow starting from a water surface elevation at the dam's Spillway Crest Level (SCL).
Space steps may be different along reservoir and depend on the variability of the cross-section, whereas computational time steps have to be equal. To avoid numerical instabilities (even if the implicit scheme is more stable), the time step should be smaller than the time of travel of the flood wave between any two successive cross-sections (Courant condition: , where c is the wave celerity).
This method is supposed to be better for the long and narrow reservoirs (and also to routing in open channels and streams), where the water surface profiles are curved due to backwater effect and a variable relationship between Q and S applies [8].

Sites and data
The first considered long and narrow reservoir (   Inflow hydrographs were chosen to produce the maximum allowable flood wave into the river channel downstream of the dam. For the Izvorul Muntelui reservoir this is the hydrograph with the peak inflow value having a 0.01% exceedance probability whereas for the Morii reservoir this is the hydrograph with the peak inflow value having a 1% exceedance probability. These hydrographs have also the same order of ratio between their inflow volume and the reservoir storage at NPL: 0.373 and 0.588 (Table 1).
For both routing methods the same discharge coefficients and geometry were considered in computations for the outlet structures.   The physical and technical characteristics of the two reservoirs and inflow hydrographs which were chosen for the analysis are presented comparatively in Table 1 [ [15][16][17][18].
One may see that the last ratio of reservoir length over maximum width L/Bmax factor differs by two orders of magnitude between the two reservoirs, due to their very different shapes.

Results and discussions
For the Izvorul Muntelui reservoir, 36 cross-section topobathymetric profiles surveyed along the reservoir were linearly interpolated at a maximum distance of 100 m in order to create the hydraulic model geometry ( figure 6). The dam has four bottom outlets and an intake for the hydropower plant (considered not to operate) and a 4 gated (Tainter gates) ogee spillway.  As mentioned at Method, the scenario considered in the present study consists of no control of the evacuated flow (spillway gates completely opened) and the initial water level at the dam being at the Spillway Crest Level.
For the Izvorul Muntelui Reservoir, in figure 8 are shown the outflow hydrographs obtained with the two routing methods, whereas in figure 9 are presented the stage hydrographs at the dam together with the difference between stage hydrographs at the reservoir tail and dam obtained with the hydraulic method.  As one may see from the above mentioned two figures and from Table 2, there is a small difference between the outflow hydrographs' peaks obtained with the two reservoir routing methods (120 m 3 /s in absolute values, which corresponds to 4.8% in relative value with respect to peak inflow value). However, there is a larger  difference between the two WSE hydrographs peaks at dam (0.96m in absolute values, i.e. 13.5% in relative value with respect to maximum difference).     difference with an average reservoir surface area (at NPL, for example). A total value of 750 000 m 3 is obtained.
The differences between the two methods in terms of peak outflow values are of 5 m 3 /s in absolute values, which translates into 6.4% relative values with respect to peak inflow, whereas the differences in terms of WSE at dam are negligible (of only 3.5 cm in absolute value, i.e. 2.6% relative to maximum difference).
On the other hand, for the Morii Reservoir, whose routing results are shown in figures 10 and 11 and Table  2, the hydraulic method gives practically no difference between the WSE at the dam and at the reservoir tail (since there is no backwater effect).

Conclusions
Reservoir routing is a mathematical procedure that routes an inflow flood hydrograph between the tail of the reservoir and the dam outflow structures. The shape of the hydrograph changes by having attenuated peaks and enlarged time bases.
For the present paper the hydrologic (Puls) and hydraulic (1D Saint-Venanat) routing methods were applied to test the efficiency of their use for two reservoirs from Romania under similar hydrologic input and same reservoir engineering data and outlet structures conditions: one long and narrow -Izvorul Muntelui reservoir on Bistrita River and one roundly shaped -Morii reservoir on Dambovita River, Bucharest.
A computer program was written in Scilab to perform the simulations of the Puls routing procedure. For the hydraulic routing Hec-RAS software was used to create the geometry of the model and make the simulations under boundary and initial conditions. Application of the two methods confirm that the hydrologic one is simpler and quicker in terms of mathematical procedure and necessary input data, whereas the hydraulic one is more costly and timeconsuming to prepare the geometry of the reservoir and dam outflow structures, based on cross-section profiles and engineering data.
As expected, results show very little difference between the two methods for the Morii reservoir, whereas for the Izvorul Muntelui, the significant backwater effect leads to important differences in computed water surface elevations at dam.
If the backwater effect is not very significant (depending on the time rate of change of flow through the system), the level pool routing method can also be applied with very similar results [8]. All the above findings confirm the conclusions of other researchers [19,6] that for long and narrow reservoirs the best routing method is hydraulic but is not worthwhile using it for others than these.
For the hydrologic reservoir routing method to accurately predict outflow and stage hydrographs (primarily relying on elevation-storage relationship), care must be taken that the storage-capacity curve must be periodically revised, since it is very sensitive to reservoir siltation. Also, for longer reservoirs, the elevation-storage curves differ between static and dynamic flow conditions.
To find quantitative criteria to better classify routing methods based on the reservoir shape (such as a critical ratio between the reservoir maximum width over length), studies like this should be extended to include more sites.