*Operational Problems in the Reservoir*

Operational problems have been present in the Stare Miasto reservoir since its construction. The first of them is caused by the internal dam splitting the reservoir into two parts. As mentioned above, the headwater level should be constant over the entire year (Figure 3a–c). The flow through the internal dam from the preliminary sedimentation zone to the main part of the reservoir is regulated by the sluice installed in the internal dam (Figure 4b,c). Due to the internal dam instability and the too small span of the sluice, equal to 3 m, the sluice was totally opened almost all the year. The problems with the internal dam increased in 2014, when huge flows in the Powa river occurred (Figure 4b) after heavy rainfalls in the upper watershed, which caused local flooding.

**Figure 3.** Stare Miasto reservoir in the Powa river: (**a**) The internal dam in 2011; (**b**) the internal dam destroyed after the flood wave propagation in 2014; (**c**) the uncovered reservoir bottom near the internal dam (author: J. Wicher-Dysarz).

High water stages and a fast-flowing flood wave reached the reservoir and met the internal dam. Due to the huge force acting on the dam during the flood wave propagation in 2014, the structure was damaged (Figure 3b).

**Figure 4.** Stare Miasto in the Powa river: (**a**) The bridge with highway A2; (**b**) the sluice in the internal dam in 2011; (**c**) the sluice in the internal dam in 2014 (author J. Wicher-Dysarz).

In 2014, it was decided that the internal dam would be reconstructed and the sediments would be removed from the upper part of the reservoir. The construction works were carried out in 2015. In the upper part the layer of sediments from the reservoir bottom was removed. The depth of the layer was 40 cm. The area of removal covered about 2 ha and the estimated volume of the removed sediments was 8000 m3. The soil removed was distributed more or less uniformly along the banks. The sluice gate was reconstructed in such a way that inverted elevation was decreased.

#### **3. Materials and Methods**

### *3.1. Available Data*

In the research four types of data were applied: (1) Hydrologic data for the Posoka gauge station; (2) the measurements of reservoir bathymetry made in 2006, 2013, and 2018; (3) a digital elevation model (DEM); and (4) sediment samples collected in 2011 and 2018.

The hydrological data for the Posoka gauge station were obtained from the Institute of Meteorology and Water Management (IMGW—Polish: Instytut Meteorologii i Gospodarki Wodnej). Although the available daily data was for the period 1971–2017 (Figure 2), the single flood wave observed between 10 April and 19 July of 2014 was used for simulations. The DEM applied in this research was obtained from the Head Office of Geodesy and Cartography (GUGiK—Polish Główny Urz ˛ad Geodezji i Kartografii). Its spatial resolution was 1 × 1 m. The vertical accuracy was 15 cm.

Examples of the data applied to create a computer model of bathymetries are presented in Figure 5. Figure 5a shows the reservoir shape and the location of the area for presentation of details. This area was chosen as a part of the reservoir near the internal dam. The bathymetry reflecting the year 2006 was reconstructed on the basis of topographic maps used for the reservoir design purposes (Figure 5b). The digital elevation model was built on the basis of elevations and isolines presented in the maps. In 2013, the bathymetry of the reservoir was measured in cross-sections with two devices: (1) GPS; and (2) an Echotrac CVM depth sounder. The total number of cross-sections was 98 and the distances between them varied from 50 to 100 m. Examples of these data is shown in Figure 5c. The last bathymetry was measured in 2018 as a part of the ISOK project (Polish: ISOK = Informatyczny System Osłony Kraju przed nadzwyczajnymi zagrozeniami, English: IT system of the Country's Protection ˙ Against Extreme Hazards). The second part of the project started in 2017 and the data from this part were applied for the presented research. The total number of cross-sections was 14 and the distances between them were in the range 200–300 m. Measurement points locations are shown in Figure 5d.

The last set of data applied in the research included sediment samples measured in 2011 and 2018. In 2011, 36 sediment samples were collected. In 2018 the measurements were repeated and 30 samples were taken from the reservoir bottom. The granulation based on 66 samples measured in 2011 and 2018 was determined as a combination of sieve and aerometer analysis, according to Polish norm PN-R-04032:1998 "Soils and minerals—samples and analysis of granulation." The sieve analysis was performed in wet conditions with a normalized set of sieves. The aerometric analysis was performed with a set of aerometric devices produced by Eijkelkamp. According to the above-mentioned Polish norm, the soil fractions taken into account were as follows: sand (2–0.05 mm), silt (0.05–0.002 mm), and clay (0.002 mm).

**Figure 5.** Examples of data used for the reconstruction of the reservoir bathymetries: (**a**) Location of the selected area; (**b**) design map from 2006; (**c**) measurement points from 2013; (**d**) measurement points from 2018.

#### *3.2. Applied Methods*

The ArcGIS 10.5.1 software developed by Esri Inc. was applied which is described in e.g., Docan [35]. It enables quite easy processing of GIS data, such as vector and raster layers. An integral part of ArcGIS is the ArcToolbox. It is a module including the main external tools and methods. Some of them are concurrent to the methods available in the basic ArcGIS interface, while others extend the capabilities of the standard interface. Extension of ArcGIS is also possible with specific plug-ins installed as ArcGIS toolbars. One such plug-in is HEC-GeoRAS [36], which is a toolbar with methods designed to support preparation of the river flow model. The tools available in HEC-GeoRAS may also be useful in preparation of a 2D model for simulation of flow in a reservoir.

The second plug-in applied is RiverBox, which is a set of geo-processing tools developed in Python by Dysarz [37]. Although the software was designed primarily for reconstruction of a river bed, after preliminary data processing it has been also successfully applied for this reservoir. The basic idea of interpolation in the channel-oriented coordinates is shown in Figure 6a. The bottom elevations are reconstructed from measurement cross-sections (green lines) through linear interpolation in two directions: the longitudinal and the transversal. The plug-in includes a number of algorithms for measurement data processing and three algorithms for reconstruction of the bed. These are described in more detail in the quoted publication [37]. The scheme of the third algorithm applied here is presented in Figure 6b. The main idea was to create a number of random points over the interpolated area, channel-oriented interpolation of the elevations in the points, and finally application of "typical" spatial interpolation. In the third algorithm of the RiverBox, the linear interpolation in the two dimensions was applied for this purpose. It wa implemented by adoption of two tools from ArcToollbox: creation of TIN (Triangular Irregular Network) and transformation of the TIN into standard raster. The choice of the third algorithm was determined by its effectiveness analyzed in Dysarz [37].

**Figure 6.** Main ideas of the bathymetry reconstruction with RiverBox [37]: (**a**) Scheme of interpolation; (**b**) applied algorithm.

Three bathymetries of the analyzed reservoir were reconstructed from available data, they represent the reservoir bed from 2006, 2013, and 2018. The first was composed from detailed design maps of the reservoir, with application of geo-referencing and digitization. The next two were reconstructed with the help of the RiverBox plug-in for ArcGIS Desktop software.

The HEC-RAS is a well-known hydrodynamic model for rivers and water reservoirs. This program was designed at the Hydrologic Engineering Center (HEC). The acronym RAS means River Analysis System, which well defines the application area. The concepts applied in the package are described by [38]. The main features of HEC-RAS are 1D simulations of flow and transport processes in river networks, including floodplains and reservoirs, as well as 2D simulations of pure flow. The HEC-RAS package also includes several useful tools for data preparation and results processing. These tools include the module for GIS data processing called RAS Mapper [39].

The 2D flow module was applied in this research. The basic model for such simulations implemented in HEC-RAS is the full dynamic wave [38] presented below

$$\frac{\partial H}{\partial t} + \frac{\partial (h\nu)}{\partial x} + \frac{\partial (h\nu)}{\partial y} = 0 \tag{1}$$

$$\frac{\partial \underline{u}}{\partial t} + u \frac{\partial \underline{u}}{\partial \underline{x}} + v \frac{\partial \underline{u}}{\partial y} = -g \frac{\partial H}{\partial \underline{x}} + \theta\_l \left( \frac{\partial^2 \underline{u}}{\partial \underline{x}^2} + \frac{\partial^2 \underline{u}}{\partial y^2} \right) - c\_f \underline{u} \tag{2}$$

$$\frac{\partial v}{\partial t} + \mu \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = -g \frac{\partial H}{\partial v} + \theta\_t \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right) - c\_f v \tag{3}$$

where the Equation (1) describes mass balance and the next two, Equations (2) and (3), are momentum balance equations. The independent variables are *t*—time and *x*, *y*—spatial dimensions. There are three dependent variables, namely *h*—the water depth, and *u*, *v*—the depth-averaged components of flow velocity. The water surface elevation *H* depends on the depth and *g* is acceleration of gravity. There are also two coefficients of more complex nature. The first is kinematic eddy viscosity coefficient *ϑ<sup>t</sup>* defined as follows

$$
\vartheta\_! = D \text{ltu}\_\* \tag{4}
$$

where *D* is non-dimensional empirical constant describing mixing intensity and *u*<sup>∗</sup> is well known shear velocity. The second coefficient is bottom friction coefficient *c <sup>f</sup>* derived from Chezy or Manning formulae as follows

$$\mathbf{c}\_f = \frac{\mathbf{g}|V|}{\mathbf{C}^2 R} \quad \mathbf{c}\_f = \frac{n^2 \mathbf{g}|V|}{R^{4/3}} \tag{5}$$

where *R* is hydraulic radius and |*V*| is magnitude of velocity vector. There are also velocity coefficient for Chezy equation *C* and roughness coefficient for Manning equation *n*.

The interpretation of mass balance in Equation (1) is relatively simple. In arbitrary control volume, the inflow and outflow fluxes described by the second and third terms of the equation should be in balance with the increase or decrease in the water volume stored, which is expressed in terms of water surface elevation changes in the first term. The left sides of Equations (2) and (3) describe the change in the momentum stored in the control volume, including the momentum stored locally and the momentum fluxes. The first terms of the right sides represent the pressure and gravity forces, while the second is responsible for internal eddy stresses and mixing. The last terms of the right sides describe the friction force generated on the contact surface of the liquid and bottom.

In HEC-RAS, the simulation may be performed with a simplified version of the model (1)–(3). If the left side of momentum equations and the eddy viscosity terms are neglected, the model (1)–(3) becomes the so-called diffusive wave model [38].

$$
\frac{\partial H}{\partial t} - \frac{\partial}{\partial x} \left\{ \beta \frac{\partial H}{\partial x} \right\} - \frac{\partial}{\partial y} \left\{ \beta \frac{\partial H}{\partial y} \right\} = 0 \tag{6}
$$

In the above equation *β* plays a role of diffusion coefficient and it is defined as

$$\beta = \frac{R^{5/3}}{n|\nabla H|} \tag{7}$$

where |∇*H*| is the magnitude of water surface gradient or, in other words, the water surface slope. The model (6) with (7) is the default choice due to stability and efficiency. The simulations presented here were performed with the diffusive wave model.

The equations are approximated with a hybrid scheme based on a combination of finite difference and finite volume methods. Such an approach is applied due to the fact that the computations mesh is composed of rectangular elements inside the flow domain and multi-edge irregular elements near the boundaries [38]. The applied boundary conditions are a flow hydrograph in the reservoir inlet and a constant stage hydrograph in the reservoir outlet.

The obtained results were transformed into the cell-oriented shear stress *τ* and stream power *SP* values. These two variables are defined as follows

$$
\pi = \rho g R S\_f \quad SP = \pi |V| \tag{8}
$$

where *Sf* is energy grade line slope (friction slope).

For the purposes of the present research, three models were developed. All of them were 2D models of the entire reservoir. Elements such as the bridge, internal dam, and main dam were reconstructed by modification of elevations in a digital terrain model representing the reservoir bottom. The basis for the 2D simulations was the numerical mesh covering the entire reservoir area, approximately equaling 1 km2. The imposed cell size of rectangular cells was 2 × 2 m, but the cells near the boundaries may have different shapes. Hence, the minimum cell area was 1.74 m<sup>2</sup> and maximum was 7.87 m2, when the total average was 4.01 m2. The total number of cells covering the reservoir was 241,699. Three bottoms were tested, which were the reconstructions of the real reservoir for 2006, 2013, and 2018. The reconstruction process is described above. In all cases the flood wave observed in 2014 was simulated with constant head water kept in the reservoir outlet. The elevation of the head water simulated was 93.50 m a.s.l. This is the so-called normal head water level in the Stare Miasto reservoir. The simulation time step was chosen due to the stability requirements and it equaled 30 s. Because the semi-implicit scheme was used for time discretization, the weighting factor was set to 0.75.

**Table 1.** The dependence between absolute roughness and characteristics of sediment grain sizes.


Besides the geometry reconstruction and flow boundary conditions, the important elements determining the flow process in the analyzed area were roughness coefficients. A simplified approach was used for this purpose. The sediment samples collected from the reservoir bottom were the basis for determination of the equivalent roughness coefficient for the entire reservoir bottom. Two step calculations were applied. In the first step, the dimensionless resistance coefficient *λ* from the Darcy–Weisbach formula was determined for full turbulent flow conditions:

$$\lambda = \left[ -2\log\left(\frac{k\_s}{14.84\ R\_{\rm li}}\right) \right] \tag{9}$$

where *Rh* is the hydraulic radius and *ks* is the absolute roughness. The second coefficient was defined theoretically as the average height of bottom irregularities. In practice its value was calculated on the basis of sediment characteristics, but a number of formulae have been proposed in literature for this purpose. The chosen and applied approaches are presented in Table 1 [40]. The symbol *dx* in this table means specific diameter of the grain sample determined on the basis of the sieve curve. The next step was calculation of Manning's roughness coefficient *n* from the formula:

$$m = R\_h^\ddagger \sqrt{\frac{\lambda}{8g}} \left[ \text{sm}^{-1/3} \right] \tag{10}$$

In Figure 7a, the box-and-whiskers plot, representing the relationship between the five methods listed in Table 1 and the roughness coefficients, is presented. The sediment characteristics were determined on the basis of the 66 sieve curves obtained from sediment samples, presented in Figure 7b, in which dots represent the mean values. The boxes show the values of the mean ± standard deviation. The whiskers denote the values between mean ± 1.96 standard deviations, which define the 95% confidence level. The plot presented in Figure 7a summarizes 330 calculations of roughness coefficients. The obtained average value was 0.038 sm−1/3, while the minimum and maximum values were 0.026 and 0.059 sm−1/3, respectively.

**Figure 7.** Estimation of roughness coefficients: (**a**) Statistics of roughness coefficients calculated on the basis of approaches presented in Table 1; (**b**) sieve curves applied for the calculation of roughness coefficients.

#### **4. Results and Discussion**

On the basis of 2D simulations, the depths, flow velocities, shear stresses, and stream power were determined for the entire reservoir with each of the three different bathymetries. Due to the huge amount of data, the results were kept in three groups: maps, graphs of selected cross-sections, and graphs with statistical analyses. In Figure 8, the locations of the main elements applied to analyze the results is shown. The reservoir was divided into five parts as presented in Figure 8a. They were, (1) the main reservoir 1, (2) the main reservoir 2, (3) pre-reservoir, (4) area under bridge, and (5) the internal dam. The four cross-sections were located in characteristic places of the greatest parts. There was 241,699 computational cells in the numerical grid. Because this number of cells was difficult for processing, the sample points (Figure 8b) were used for monitoring the obtained results. The total number of sample points was 1670, which was about 0.67% of the available locations of computed values. Although this number was much smaller than the number of computational cells, the sample points covered the entire reservoir with sufficient density. Because the reservoir parts differed in size, their distribution over the parts was not uniform. In the main reservoir 1, the main reservoir 2, and pre-reservoir there were 500 points, while 20 points were generated in the area under the bridge and 150 points generated over the internal dam.

**Figure 8.** Elements used in the analyses: (**a**) Parts of the reservoir and location of cross-sections; (**b**) location of sample points.

**Figure 9.** Analysis of depth distributions over the reservoir bottom in the Stare Miasto reservoir: (**a**) Bathymetry 2006; (**b**) bathymetry 2013; (**c**) bathymetry 2018.

The maps presented in Figure 9, Figure 11, and Figure 12 show the spatial distribution of depth, velocity, and shear stress over the three reservoir bathymetries, namely 2006, 2013, and 2018. As may be observed, the reservoir bottom changed during 12 years. Figure 9 presents the maps of maximum depths obtained as the results of simulations. The following maps (Figure 9a–c) show the results for bathymetries of 2006, 2013, and 2018, respectively. The comparison of bathymetries is made in Figure 10—part A (Figure 10(a1–a4)), where the geometry of cross-sections marked in Figure 8 is shown. In part B of Figure 10 (Figure 10(b1–b4)), the differences of the bottom elevations related to the initial bathymetry 2006 were analyzed. The blue line denotes increase in the elevations, from the initial bathymetry to the bathymetry measurement in 2013. The red lines represent differences between measurements of 2018 and the initial data.

**Figure 10.** Selected cross-sections: (**a**) Bottom and water surface elevations for three bathymetries; (**b**) difference in bottom elevations between 2006 and 2013, as well as between 2006 and 2018.

As it can be seen, during the life of the reservoir, since 2006, the reservoir bottom was considerably transformed. Figure 10, graph b1 shows the cross-section located in the inlet of the Stare Miasto reservoir. The differences in the bottom elevations varied from 2 m of accumulated sediments to 2 m of removed sediments. The comparison with the bathymetry of 2006 shows that the greatest differences occurred in 2018 as a result of natural processes, as well as artificial removal of sediments in the end of 2014 and the beginning of 2015. At this time, a 0.4-m layer of sediments was removed from the bottom of the upper part. Figure 10(b2,b3) show similar bottom increases for the bathymetries of 2018 and 2006. In the last cross-section located near the main dam, the sediments were accumulated during 12 years of reservoir life. The only place where erosion occurred was located near the valve tower and spillway.

The graphs presented in Figure 11 illustrate the spatial distribution of maximum magnitudes of velocity in the reservoir. According to the maps, the maximum velocities occurred near structures such as the internal dam, the sluice in the internal dam, the flow below the bridge, and the reservoir outflow. There were also greater velocities in the reservoir inlet, though the area is so small that it was not visible well in the scale of Figure 11. The internal dam separates the preliminary part of the reservoir from its main part. The single and very small sluice caused local increase in the velocity magnitude and erosion of the alluvial bed. The maximum velocities over 3 m/s occurred also below the bridge. The reason was the too narrow bridge span.

**Figure 11.** Spatial distributions of velocities in the Stare Miasto reservoir for three bathymetries: (**a**) Bathymetry 2006; (**b**) bathymetry 2013; (**c**) bathymetry 2018.

The results presented ealier were confirmed by the results shown in Figure 12a–c, illustrating that the distribution of the shear stresses in the reservoir was similar to that of velocity. The maximum values occurred at the reservoir inlet, near the internal dam, and in the bridge span. Additionally, the reservoir inlet was more prone to shear stress increase. The maximum shear stresses occurred for the bathymetry of 2013 (Figure 12b). In the bathymetry of 2018, the corresponding values were lower, especially, in the pre-reservoir. Supposedly, the artificial change in the bottom configuration as a result of dredging in 2014/2015 was the reason for such a change. Additionally, the very specific cross-sections marked in Figure 12 are cross-sections near the inlet, in the internal dam, and under the bridge. They were applied for analyses presented in Figure 13.

**Figure 12.** Shear stress calculated for the Stare Miasto reservoir and three bathymetries: (**a**) Bathymetry 2006; (**b**) bathymetry 2013; (**c**) bathymetry 2018.

**Figure 13.** Shear stress and stream power in chosen cross-section: (**a**) inlet; (**b**) internal dam; (**c**) bridge.

The next variable analyzed was the stream power [41–43], but it is shown only in the graphs created with one of two methods: (1) as cross-sections with results extracted from results in raster format (Figures 8a and 12) (2) as statistics of values in the random sample points (Figure 8b). In Figure 13a–c, the graphs of shear stresses (SS) and stream power (SP) in selected cross-sections are shown for comparison. The cross-sections applied were those shown in Figure 12. In Figure 13a, the inlet cross-section is presented. The shear stresses (SS) simulated there were over 20 N/m2, while the stream power (SP) values reached about 100 W/m. The highest values of shear stress as well as stream power were observed for the reservoir in 2013. Supposedly, the accumulation of sediments in the period 2006–2013 changed the configuration of the bottom in such a way that there was an increase in these two factors. The values of shear stress and stream power were lower for the conditions in 2018. The removal of sediments in the pre-reservoir in 2014/2015 induced another kind of change. The next specific cross-section is shown in Figure 13b, at the internal dam with its sluice gate. The shear stresses as well as stream power in the sluice were increasing rapidly, independently of the bottom configuration tested. However, the sluice gate is made of concrete and it should be resistant to such forces. The huge values of the analyzed variable are more dangerous for the rest of the internal dam. The maximum values of shear stress were about 40 N/m<sup>2</sup> there, while the stream power reached about 200 W/m. The values of these two factors were relatively small for the conditions in 2006, but in other configurations of the bottom, 2013 and 2018, high shear stresses and stream power were noticed. These values were quite high, especially in the lowland reservoirs where fine sand is deposited with grain sizes smaller than 0.5 mm. It means that, irrespective of the dredging and rebuild of the sluice gate, the internal dam is still exposed to huge forces and is prone to break. In the cross-section of the bridge, the shear stresses were smaller, about 7 N/m2 and seem to be stable. The same was noticed about the

stream power values obtained there. The hydraulic conditions under the bridge seemed to be less dependent on the bottom configuration.

The next stage of the research was statistical analysis performed with Statistica 13. The values of maximum magnitudes of velocity, maximum shear stresses, and maximum stream power were read in the locations of the randomly generated points shown in Figure 8b, and then processed over the set of points and subsets included in particular reservoir parts (Figure 8a). Such results are presented as standard box-and-whisker plots in Figure 14a–c. The graph in Figure 14a shows the maximum flow velocities for different parts of the reservoir indicated in Figure 8a. The total medians were 0.32, 0.41, and 0.44 for bathymetries 2006, 2013, and 2018, respectively, while the mean values were 0.53, 0.67, 0.73 m/s. The significant differences between median and mean values for all bathymetries were caused by greater velocities calculated for the bridge and the internal dam. In the first case, the medians were 4.69, 4.81, and 6.38 m/s for the same bathymetries as previously. In the internal dam these factors reach values 2.33, 3.65, and 4.11 m/s. Such great irregularities indicate the importance of these two locations. The greatest values for all bathymetries were noticed under the bridge and over the internal dam. The velocities simulated with the initial bathymetry were slightly smaller than the others. The similarity of the internal dam median velocities for 2013 and 2018 shows that the dredging in 2014/2015 did not improve the safety of this dam much. Fortunately, the maximum velocities decreased from about 10 m/s in 2013 to 6–7 m/s in 2018. In the other parts of the reservoir, the velocities obtained were much smaller irrespective of the bathymetry tested.

The trends observed in the spatial distributions of shear stresses (Figure 14b) and stream power (Figure 14c) were different. Definitely the highest values of these two variables were observed in the internal dam. The discrepancy of the velocity magnitudes with shear stress and stream power observed under the bridge may be explained by relatively great narrowing of the flow area but small hydraulic slope governed by the head in the main dam. The internal dam influenced from the upstream and downstream was in different conditions. Hence, the shear stress and stream power confirm the observations of velocities.

**Figure 14.** Characteristic values of analyzed hydraulic parameters simulated for bathymetries 2006, 2013, and 2018, and recorded at random points: (**a**) Velocity; (**b**) shear stress; (**c**) stream power.

It is also important that there is a difference in trend between the median and maximum values. The median values of the shear stress as well as stream power more or less decreased between 2013 and 2018. The simulation over the bottom for 2013 provided median shear stress equal to 7.95 N/m2 and median stream power of 27.97 W/m over the internal dam. The median results obtained for bathymetry 2018 were 7.77 N/m2 and 19.27 W/m, respectively. Hence, the decreases were 2.2% and 31.1% for shear stress and stream power, respectively. The same results indicate significant increases of the maximum values in the case of these two factors. The maximum shear stress and stream power values for bottom for 2013 were 28.63 N/m2 and 653.78 W/m, respectively. The results for the bathymetry 2018 gave maximum values of the same variables equaling 36.75 N/m2 and 815.12 W/m. Hence, the increase of shear stress was 28.4% and the increase of the stream power was 24.7%. It confirms the previous conclusion that the rebuild of the internal dam and sluice gate did not improve the safety of this object.

The analysis of bed changes and shear stress or stream power in reservoirs is relatively difficult for modeling. In many papers there is a description of the shear stresses acting on the river bed in mountainous as well as lowland conditions [30,44,45]. The majority of such concepts is based on the simulation in 1D or calculated with empirical formulae [30,45]. In results provided by Glock et al. [30], who simulated flow in Austrian rivers and a laboratory flume, shear stresses varied from 0.4 to 93 N/m2. The course of the channel, meandering or straight, influenced the obtained values greatly. The results were also dependent on the type of model applied. They compared results of 1D, 2D, and 3D models. It is obvious that shear stresses are impacted by the roughness of the channel bed and

flow velocities. In the Stare Miasto reservoir maximum velocities up to 10 m/s (Figure 14a) generated huge shear stresses, reaching values of 90 N/m2 (Figure 13b). Such magnitudes were noted in the internal dam and near the sluice installed there. Huge velocities and shear stresses caused huge values of stream power. In the sluice it was 900 W/m (Figure 13b). During the flood wave propagation in the reservoir, the most critical location was the reservoir inlet. The velocities simulated there equaled 5 m/s (Figure 11), which was related to huge shear stresses reaching values of 25 N/m2 (Figure 13a).

### **5. Conclusions**

The impact of hydraulic parameters, such as velocity, shear stress, and stream power, on the operational conditions in a two-stage reservoir was analyzed. The basis for the investigations was simulations assuming the 2D hydrodynamic model. Such an approach enabled better illustration of the processes and changes in the reservoir for different moments of time. The method applied was different from the 1D models usually applied, as it enabled better recognition of the advantages and drawbacks of two-stage constructions. The shear stresses were calculated for the entire reservoir basin and especially for the structures located in it, and visualization revealed significant differences in response to impacts such as flood wave propagation. The simulations were made for the bathymetries of 2006, 2013, and 2018. The conclusions drawn from the results indicated really great stresses exerted on the internal dam, which may lead to breakage of the dam structure. In general, the analyses proved that structures in the reservoir are prone to profound stresses. The magnitudes of shear stresses in specific flood conditions are not observed frequently in lowland rivers and reservoirs. The contraction caused by the bridge span is also very dangerous. In flooding conditions, it may cause great acceleration of water flow. The highway bridge is a massive construction and no direct threat related to great stress has been detected in this case. However, it is not possible to exclude the probability of erosion below the bridge and collapse of the construction due to bridge scour in future. One thing is obvious—the velocities, the shear stress, and stream power values are not typical of lowland conditions if the flood propagates through the reservoir. The presented analyses confirm that two-stage reservoirs are a good solution of some reservoir problems, if the inflows are moderate. However, the increase in inflows up to flood magnitudes may generate new problems related to the stability of the structures located in the reservoir basin.

The applied 2D hydrodynamic modeling enabled more detailed analysis focused on spatial distributions of important factors: the velocity, the shear stress, and the stream power. The obtained results are satisfactory and they may be applied in several ways. The approach presented enabled identification of potential threats, which is well shown on the example of the internal dam. However, such an approach required greater computational costs than in other cases—more simulation time, faster processors, more effective models, more memory on disks, etc. Although, the use of three different bathymetries and comparisons between them increased the complexity of the problem, it also enabled simplified verification of the object vulnerability to flood hazard. As indicated, the rebuild of the internal dam did not improve the safety of this structure and future breakage could be expected. Some application areas of the results presented seem to be obvious—improvement of reservoir design, optimization of structure shapes, invention of protective elements reducing detected threats. Others are not so obvious (e.g., operational control of reservoir to safely propagate the flood wave). However, such an approach requires greater computational power than that applied in this research.

As can be seen, the two-stage reservoirs are quite an interesting alternative approach to prevention of reservoir sedimentation. However, their construction should be carefully analyzed and computationally tested before the main dam is built. The results of the present research could help in the design of crucial elements, such as the sluice in the internal dam, the span of the bridge, etc. All such elements could cause unexpected effects dangerous for the stability of the reservoir construction and safety of the water management in the reservoir.

**Author Contributions:** Conceptualization; methodology; software; validation; formal analysis; investigation; resources; data curation; writing—original draft preparation; writing—review and editing; visualization; supervision; project administration; and funding acquisition, J.W.-D.

**Funding:** This research received no external funding.

**Acknowledgments:** This research would not be possible without the access to the data provided by such institutions as IMGW, CODGiK, and GUGiK. The access to the free data available in the Internet, such as Geoportal 2, was also very helpful.

**Conflicts of Interest:** The author declares no conflict of interest.
