imPlementAtion of the Adh hydrodynAmiC model on the WłocłaWek ReseRvoiR

The variation of water velocity in an artificial dam reservoir is influenced not only by the inflow discharge, but also by the bathymetry of the reservoir and the water level at the dam. The depiction of spatially complex variations in flow velocity through a reservoir would not be possible without the use of hydrodynamic models. A reliable hydrodynamic model of the reservoir is an effective tool for predicting and analyzing changes in the reservoir geoecosystem in an age of changing climate and risk of water stress. A depth-averaged two-dimensional AdH model was used to visualize the hydrodynamics of the Włocławek Reservoir. Running the model for eight different hydrological conditions delivered consistent results and allowed to calibrate the model parameters. Additionally, it provided a way to verify the data regarding the rating curve of the Vistula River upstream the reservoir.


Introduction
The close relationship between the natural processes occurring in artificial reservoirs and the hydrological regime of the river means that their functioning depends on the variability of the flow of the rivers recharging the reservoir and changes in water level on the dam (Kimmel & Groeger, 1984;Straškraba & Tundisi, 1999;Friedl & Wüest, 2002;Kennedy, 2005;Ye et al., 2007;Xu et al. 2009;Zhang et al., 2010;Hayes et al., 2017).These changes affect not only the flow velocity of the water, the rate of water renewal, and the water circulation, but also the amount of river sediment loads that has entered the reservoir, the physicochemical properties of the water, the structure, and abundance of aquatic organisms, and the rate of biological production (e.g., Straškraba & Tundisi, 1999;Wagner & Zalewski, 2000;Friedl & Wüest, 2002;Ahearn et al., 2005;Wei et al., 2009;Vörösmarty et al. 2010;Liermann et al. 2012;Zhao et al., 2013).Disruption of continuity of the channel processes by dams cause modifications, in terms of rate and direction, of the geomorphological and sedimentological processes that shape the bed and banks of the river channel below the dam, in the reservoir, and upstream its backwater (e.g., Babiński, 1994;Phillips, 2003;Walling & Fang, 2003;Vörösmarty et al., 2003).The main changes in the functioning of river ecosystems are visible above dams, where a greater or lesser transformation of riverine/lacustrine conditions occur due to a decrease in the water table longitudinal slope.The unit power of the river stream, determining the transport of river sediments and geomorphological activity, becomes less significant, making way for wave, wind, and density currents, among others.As a result, an internal water circulation system, independent of the river flow, is created in the lower parts of the reservoirs (e.g., Timčenko, 1989;Avakân et al., 1994;Litvinov, 2000;Matta et al., 2017).The importance of the streamflow current, however, is still significant in the backwater part of the reservoirs and, in the case of rheolimnic reservoirs, along their entire length.However, even in such reservoirs, in the parts covering areas outside the former channel, the water circulation is independent of the river flow (Širokov & Lopuch, 1986;Kajak, 1998, Dubnyak & Timchenko, 2000).
The authors decided to use the 2D depthaveraged AdH (Advanced Hydrodynamics) model described in detail later in this paper.It has been successfully applied to three other Polish dammed reservoirs: Tresna, Dobczyce and Goczałkowice (Hachaj & Tutro, 2014;Hachaj, 2018;2019).This model allows simulating estuarine and riverine flows, hydrodynamics in reservoirs, and flow-through lakes, flows due to dam breaches, flood flows and more.Additionally, it gives the user an option to include: • Wind impact including changes of the wind velocity in space and time and wind induced surface waves; • Subglacial flow in selected regions as well as the impact of ice floe; • Drying and wetting of the land along with the changing water surface level; • Suspended particles transport, with the overlay PTM model (Macdonald et al., 2006); • Impact of hydrotechnical structures and floating vessels on the time-dependent velocity field.
Although the Włocławek Reservoir has been in operation for more than 50 years, its hydrodynamics is still not sufficiently characterized.In studies on various aspects of the reservoir's functioning, only average values of flow velocities in cross-sectional profiles calculated with empirical formulas have been most frequently used (e.g., Banach, 1985;Grześ, 1983).Instrumental measurements of speed and direction of water flow were conducted incidentally in the limnic part of the reservoir (Wierzbicki & Ujda, 1986), and in the upper-river part of the reservoir during river icing (Majewski, 1987).More systematic hydrological measurements, in various hydrometeorological situations and sections of the reservoir, were conducted from 2000 to 2005 (Gierszewski, 2018).The simulations of water flow conditions through Włocławek Reservoir carried out so far with the help of twodimensional hydrodynamic models CCHE2D (Bogucka-Szymalska & Magnuszewski, 2007) and Aquadyn 3.1.(Gierszewski, 2006(Gierszewski, , 2018) ) for average and low flows of the Vistula River into the reservoir, have not been verified by real values of water flow velocities.Hydrodynamics of the reservoir were also indirectly characterized by lithodynamic analyses of the sedimentary environment (Bogucka & Magnuszewski, 2006;Gierszewski et al., 2006;Gierszewski, 2018).
The studies presented in the article have completed and detailed the existing knowledge of the hydrodynamics of the Włocławek Reservoir.Their main objective was to visualize the reservoir hydrodynamics for any real conditions.It was essential to achieve the best possible compatibility of the modeling results with the results of velocity measurements during low flows of the river Vistula.Under such hydrological conditions, the reconstruction of the river ecosystem into a lake ecosystem is the most complete (Straškraba, 1999;Błędzki & Ellison, 2000;Soares et al., 2008;Cunha-Santino et al., 2016;Shivers et al., 2018).A reliable hydrodynamic model of the reservoir that also considers its operation under such conditions will be an important tool for predicting and analyzing changes in the reservoir's geoecosystem in an age of a changing climate.

Study area
The Włocławek dam was put into operation in 1970.In its assumptions, it was to be the first of the eight dams forming the cascade of the lower Vistula.This 1,200-meter-long structure consists of five elements: a ground dam, weir, 160 MW hydroelectric power plant, fish passage, and a navigation lock, as well as a side dam on the left bank of the Vistula.As a result of damming up the Vistula by 14 m, the largest dam reservoir in Poland in terms of area (75 km 2 ) was created.The length of the reservoir is 57 km, the average width 1.2 km, the average depth 5.5 m, and the maximum depth is about 15 m.Bigger depths occur in the line of the former thalweg of the river, running mostly along the right bank (Fig. 1).Shallows are associated with inundated sandbars and Vistula islands, as well as with the left floodplain (Gierszewski, 2018).The current capacity of the reservoir is estimated at 370 million m 3 .The Włocławek Reservoir is a typical valley one in its nature, it is narrow and long, with little diversified banks.The area of the reservoir covers the former river channel of the high water along with the flood terraces of the left bank (Fig. 1).Only 14% of the reservoir's area coincides with the part unrelated to the former Vistula River channel.
The Włocławek dam is currently used mainly for hydropower purposes, while the reservoir is used for recreational purposes, as well as, to a small extent, for flood protection.The average volume of water inflow to the reservoir is 931.8 m 3 •s -1 (in the years 1971 to 2015), and is slightly higher than the outflow, which reached 905.7 m 3 •s -1 .The bigger flow differences between the inflow to, and the outflow from the reservoir are associated with the functioning of the reservoir during flood periods, or larger discharges of water through the dam in order to improve navigation.In the Włocławek Reservoir, as in the case of other channel-type reservoirs with a hydroenergetic function, water levels are usually balanced.Their annual fluctuations in the lower part of the reservoir during its functioning ranged from 52 cm to 187 cm, averaging 101 cm.The high stability of water levels in the reservoir shows that on average, over 203 days a year, their daily amplitudes are less than 5 cm (Gierszewski, 2018).

Bathymetric map
Bathymetric measurements of the Włocławek Reservoir were made on from September 29 -October 3, 2012.The work was carried out from a motor pontoon equipped with a Teledyne ODOM E2 single beam echosounder and a Trimble 5800 GPS GNSS receiver.Data was collected during a stable low flow (from 353 to 367 m 3 •s -1 ) in cross sections separated by approx.100 m.Receiver.GPS GNSS recorded the positions and ordinates of the water table elevation on an ongoing basis.To obtain a digital elevation model (DEM) with bathymetry of the river (Fig. 1), the measurement data was imported to SAGA GIS v.2.3.2.where spatial interpolation of the distributed data was performed using the multilevel B-spline algorithm method proposed by (Lee et al., 1997).The developed DEM with bathymetry has a resolution of 3.0 x 3.0 m.

Flow velocity measurements
An Acoustic Doppler Current Profiler by Aanderaa was used to measure water velocity and flow directions.The measurement was made in verticals every 1 m in depth for 30 seconds.The velocity values that were to be used to calibrate the 2D hydrodynamic model are the average velocities from each measurement vertical.The results from 208 hydrometric verticals were used in this study.The measurement results represent conditions corresponding to flows ranging from 240 to 2320 m 3 •s -1 and damming level of 56.89 to 57.38 meters above sea level.

Data sources
Average daily values of Vistula water discharge from water gauging cross sections in Kępa Polska and Włocławek were collected from publicly available resources of the Institute of Meteorology and Water Management -National Research Institute (IMGW PIB).Average daily values of water levels at the dam were taken from State Water Holding Polish Waters.

The AdH model
The AdH (Adaptive Hydraulics model) modeling tool (Berger et al., 2013) has been created by The US Army Corps of Engineers and it is available in the SMS (Surface-water Modeling System) package distributed by the AQUAVEO company (AQUAVEO, 2018).AdH belongs to the group of two-dimensional depth averaged (Whipple, 2004) finite element modeling tool.This hydrodynamic model is based on adaptive constrained mesh (Branets et al., 2009).AdH is a natural successor to the previous models from the SMS package.
As the name of Adaptive Hydraulics model implies, this model adapts the computational mesh in the regions where it is necessary to obtain convergence.This is done by refining (i.s.splitting) the elements into smaller triangles.Such a procedure generates additional vertices that are used in the given time step.The time step itself can also be automatically reduced if needed, then released to its original value.
The basic input data for this model include: bed shape (bathymetry); bed friction parameters (interpreted as Manning roughness coefficient); initial water surface level; physical parameters of the water (density, viscosity); inflows and outflows discharges (steady, unsteady, surface level dependent).Additional input may include wind impact, ice cover and more (Berger et al., 2013).
Finite elements in AdH are triangles with discretisation nodes located in their vertices.The main result of the model run is a twodimensional planar field of the mean average horizontal velocity of water.Vertical movement and other vertically dependent phenomena are neglected.This field may be static or time-dependent according to the conditions imposed upon the domain.
The AdH model solves numerically conservation equations of mass (continuity equation) and momentum (1).The continuity equations use the cross-section area and the discharge in any given cross-section as the variables of interest (2): where: A -cross-section area (m 2 ), Q -discharge (m 3 •s -1 ), U -vertically averaged velocity of the flow (m•s -1 ), H -depth (m), B -section width (m).
The mass conservation equation (3) for any given planar element is: (3) where: v x , v y -velocity components in x and y directions, respectively (m• s -1 ).
The AdH model uses the finite element method to solve two-dimensional momentum conservation equations for water flow (4)in Eulerian frame for both x and y directions (Berger et al., 2013;Witek, 2013).These equations are as follows: If necessary, these basic equations may be both simplified or enhanced.For standard calculations of water flow in retention lakes, the external pressure, treated as atmospheric, is considered to be constant over the reservoir's water surface and the Coriolis effect is neglected (is small for most water reservoirs).Consequently, it allows to ignore the mentioned factors and simplify the equation.On the other hand, other forces caused by wind or ice may be added if needed.The program solves the equation along the streamline upwind (Petrov-Galerkin) scheme.It uses the Runge-Kutta method with Cash-Karp error control routine which produces a fifth-order accurate solution (Savant & Berger, 2012).
As its input the model takes the bathymetric map of the flooded area, water surface elevation at its main outflow (Dirichlet boundary condition), and discharge values at all the inflows and secondary outflows, if there are any (Neumann boundary conditions).All the boundary conditions may be constant, time-dependent or mixed.As standard starting condition the model usually takes flat water surface at the elevation not lower than where: v x , v y -velocity components in x and y directions, respectively (m•s -1 ), f x -unitary force in the x direction, represented by wind-induced surface friction, wave impact (m•s -2 ), f y -unitary force (wind, waves…) in the y direction (m•s -kinematic viscosity coefficient (m 2 •s -1 ), p -external pressure (Pa), ρ w -water density (kg•m -3 ).
the downstream boundary condition.It is also possible to use a previously obtained solution as the set of starting conditions for a new calculation.The results of the model yield: scalar water surface elevation field and vector planar water velocity field (averaged along the vertical dimension).
In principle every constant present in the governing equations of the model may be modified if necessary for calibration and validation purposes.But there are no published studies involving changes of gravitational acceleration, or freshwater density.If Coriolis effect, wind impact, and differences of the atmospheric pressure are neglected (as they are in this study), the parameters remaining for calibration are water viscosity and bed roughness.The latter is accessible to the user as the Manning roughness coefficient n.The Chezy coefficient is then derived from the Manning one via the formula (5): (5) where: n -Manning roughness coefficient (m For the purpose of this study a fragment of the Vistula River from the railway bridge in Płock (km 629+500) to the dam in Włocławek (km 674+850) was mapped in the model as a free mesh of 19,828 triangular elements.The solution stability has been checked for the water surface level in range from 56.89 to 57.38 m a.s.l. at the dam, and for discharges up to 2,320 (m 3 •s -1 ).Two secondary inflows were introduced (as indicated in Fig. 2), each of them assumed to carry 1% of the total discharge.

Model calibration
The calibration was performed on 8 velocity measurement series summarized in Tab. 1 and having 10 to 37 measurement points each.Figure 2 shows the distribution of the measurement points for 4 series exposed on the calculated water velocity maps for appropriate conditions.The calibration target variable was the vertically averaged velocity of the water column in given points.This value is a direct output from the model, and it has been measured in each of the points as described above.The target function to minimize was Root Mean Square Error (RMSE) of computed versus observed velocity.To avoid misinterpretation, as RMSE may be overinfluenced by a few highly divergent values, Mean Absolute Error (MAE) was also calculated for each calibration run.Plots  were also created to provide a visual safeguard against misreading the calculated results (Anscombe, 1973).
The input variables of calibration were two global model parameters: turbulent viscosity and bed roughness, and one boundary condition parameter: water discharge (Q) as total inflow.The former two have their default values in the model, which were taken as the starting value for the calibration procedure.For the last one the starting value was taken from the reservoir management log for specific days.
As for the turbulent viscosity coefficient, the calculations proved that the model is to a great extent insusceptible to changes of this parameter.Such insensitivity allowed to exclude the viscosity coefficient from further analysis.In all the further calculations the default value of the turbulent viscosity coefficient, 1.13•10 -6 (m•s -2 ), was used.This value is then also recommended for future calculations.
For the purpose of calibration, the roughness coefficient was changed with the interval of 0.005 starting from its default value of 0.025.The total inflow adjustments were taken with the 5% increment starting from the measured value.The best fits (i.e., the results of the lowest RMSE) are shown in Tab. 1.The plot of fitted variables versus total discharge is shown in Fig. 3. Comparison of MAE for all the results and visual checking appropriate calculated versus observed velocity plots did not reveal anomalies.
The average adjustment of the inflow discharge ranges from -10% to +40%, with the average of +11.8%.Following the 5% step of inflow adjustment, the authors recommend a +10% correction for further calculations.
The bed roughness coefficient of the best fits varied in range from 0.030 to 0.045 with the average of 0.036.Although it is possible that a minor declining trend is present there, the number of data points is too low to fully back such a claim.The authors thus recommend using the value of 0.035 for all future analyses.
In order to compare errors between different series, relative RMSE was calculated by dividing RMSE by total inflow (last row of Tab.1).This value remains in range from 6.57•10 -5 to 16.12•10 -5 with a slight decreasing tendency with the total discharge.

Interpretation of the results
Out of two global model parameters: the turbulent viscosity coefficient and the bed roughness the model appeared to be invariant with respect to the first one.
As for the bed roughness coefficient, the calibration runs resulted in fitting the value of 0.035, which is higher than the default (0.025).
To the authors it appears to be a numerical case, and not a physical one.While in the governing equations of the model this coefficient is treated as the Manning one, it should not be interpreted simply as such.For low velocity flows (millimetres or centimetres per second), bed interaction formulas provide the main damping not only for physical effects but also for numerical instabilities.This is why the fitted bed roughness coefficient may be higher than one expected from Manning coefficient tables.Additional fact that supports this interpretation is, that the fitted roughness coefficient slightly decreases with the discharge: For higher discharge values the best fit for bed roughness was 0.030, and lower ones it was in the range from 0.035 to 0.045.The most probable cause is that for lower velocities the numerical artefacts in the calculated velocity field are relatively higher, thus require a stronger damping not to propagate into subsequent iterations.The R 2 value for a linear approximation of this decrease of the roughness coefficient vs total discharge is 0.36.
When taking into account the upper and lower boundary conditions, one can assume the lower one, i.e., the water surface level measured at the Włocławek dam, to be accurate enough not to require additional tuning.Indeed, this value is provided with a 1 cm accuracy while model sensitivity checks indicate that changing it even by several centimeters does not induce significant changes in the modeled velocity field within the reservoir.However, the upper boundary condition, i.e., the main inflow discharge from the Vistula River should be treated with caution.It is estimated via the rating curve at Kępa Polska water surface level measurement post.Such measurements have the accuracy good enough for hydrologic purposes, but too low to serve as firm data for model calibration runs.Instead of treating the inflow as fixed data, it should be treated as a variable to fit during the calibration procedure.As an added value, calibrating a numerical model with flexible upper boundary condition provides a verification of the rating curve being currently in use at the given measurement post.
For each of the 8 calibration series the inflow values were tested with a 5% increment/decrement, relative to their original values.The correction, for which the best fit was obtained, ranges from -5% to +40% with the average of +11.875%.It leads to a supposition that the discharge values obtained from the rating curve are on average about 10% lower than the real ones.
This result should be verified during the next stages of modeling when wind effects are taken into account.Indeed, wind may cause local currents that increase velocity without the need of rising the inflow value.It is noteworthy however, that the highest percentage of discharge adjustment was obtained for measurement series 7, for which the velocity measurement points were located close to each other in the middle part of the reservoir (see Fig. 2D).Additionally for series 7 the wind impact was relatively small: the average wind velocity was 1.05 m•s -1 , which is lower than the wind velocity average for all the series (2.27 m•s -1 ).
Within such a setup local currents would be visible if present, but it was not the case.
The value of velocity RMSE relates directly to the absolute value of velocity, which is -on average -proportional to the total discharge.In order to compare this value between the series and to assess the model accuracy for higher and lower discharge values, the "relative RMSE" was calculated for each series by dividing the minimal RMSE value by the total inflow for each series.Results are shown in Tab. 1 and Fig. 3.All the errors are within the same order of magnitude; the difference between the highest and lowest values is less than threefold.The relative RMSE seems to slightly decrease with the total discharge, however the R 2 value for the linear fit of this dependence is only 0.26.It means that the model can be used for a wide range of the water discharge within the reservoir.

Model limitations
The authors consider the model being able to reproduce the reservoir's hydrodynamics properly with the proposed values of parameters and adjustments.However, it should be stated that the prediction of the models has their limitations; some of them are intrinsic, while others can be worked out during future research: • No predictions can be better than the input data.This general rule should be taken into account when comparing numerical models of field objects with models of laboratory ones.For the former the level of data uncertainty is usually much higher, thus the tolerance for modeling accuracy should be set higher as well.
• The AdH model is a two-dimensional, planar one.As such all the vertical flow components are neglected; this is an inherent feature of such models.Some 3D-related effects may still be simulated to some extent (Hachaj & Szlapa, 2017) and the vertical velocity component may be artificially adjusted, if necessary, particularly for the purposes of sediment transport (Wilk et al., 2022).• It is not certain whether the calibrated parameter values can be used for discharges significantly out of the calibration interval (from 240 to 2320 m 3 •s -1 ).
To assess this issue, additional measurement sessions of real velocity would have to be conducted.• The current approach does not include any wind impact which is one of the driving forces of reservoir hydrodynamics, especially for low discharge values.Calibrating the model to properly consider the wind impact is a natural next step of this research.

Conclusions and perspectives
Two-dimensional depth-averaged Adaptive Hydrodynamics (AdH) model was implemented at the Włocławek Reservoir.The stability and sensitivity of the implementation were both checked.Model predictions regarding depth-average planar water velocity field were compared with real measurement data.The optimal values of model global parameters (bed roughness and turbulent viscosity) were found.Additionally, a correction towards the inflow discharge measured at Kępa Polska was proposed.
The model allows simulations of the flow within the Włocławek Reservoir working under different hydrologic conditions.The range of possible values of water surface elevation and inflow discharge parameters corresponds with their typical values changing from low flow to flood circumstances.The output of the model is an unsteady planar velocity field of the water contained in the reservoir calculated for the assigned conditions.
The research described in this paper can be treated as the first step, or rather a base of a longer program.Further modeling work is intended to perform, including: • Assessment of the wind impact on the hydrodynamics of the reservoir, calibration of wind-related parameters; • Preparation of the "winter" model of the reservoir with ice cover or ice floe on the water surface, calibration of ice-related parameters; • Checking the impact of different local bed roughness (especially caused by aquatic vegetation along the banks) on the reservoir hydrodynamics, introducing local corrections of the model parameters; • Research regarding water exchange and retention time in the reservoir working under different hydrometeorological condition (identifying potential hydrodynamic settling basins inside the reservoir); • Application of the overlay sediment tracking model (PTM), modeling of sediment accumulation and passing through; • Scenario-based approach of simulating extraordinary events within the reservoir, like propagation of flood waves, spreading of pollutants, significant lowering of water surface level, dam removal, etc.;

Figure 1 .
Figure 1.Digital Elevation Model of the Włocławek Reservoir

Figure 2 .
Figure 2. Sample calibration series: I.The distribution of the measurement points for 4 series exposed on the calculated water velocity maps for appropriate conditions (A-series 1, B-series 3, C-series 4, D-series 7; see Tab.1).II.Plots of calculated versus observed velocity for these series

Figure 3 .
Figure 3. Best fits of calibration variables and relative RMSE for all the series plotted versus the total discharge

Table 1 .
Summary of all the measurement series.