Mathematical modelling of flow and sediment transport – the Tisa dam case study

High flow velocities downstream of the Tisa Dam at Novi Bečej (Serbia), produced by the large difference between water levels in the upstream and downstream reaches, cause bank erosion by the right bank immediately downstream of the emergency spillway, threatening stability of the dam itself. Due to complex interaction between the parameters of flow and sediment transport at the given location, an approach involving a spatial (3D) mathematical model of flow and sediment transport, aimed at solving such complex processes, is investigated in this paper. An optimum solution meeting all requirements for the stability of structures in question has been reached by the model, through comparison of the results generated by the existing conditions with the result provided by the measures proposed to fix the problems.


Introduction
The Danube -Tisa -Danube (DTD) Canal is a system of canals created to serve several purposes such as the flood control, drainage, irrigation, navigation, tourism, hunting, and fishing.The canal system spreads over a surface area of 12700 square kilometers between the Danube and Tisa rivers in Vojvodina, Republic of Serbia.The total length of the system is 929.0 km, and it includes 51 water management structures and 180 bridges.The largest and one of the most important structures in the Danube -Tisa -Danube system is a dam on the Tisa River built in 1977 near the town of Novi Bečej (Figure 1).The dam is located at the river kilometer 62+988 from the mouth of the Tisa discharging into the Danube, 2.4 km downstream of Novi Bečej.Its main purpose is to store the water in the upstream reaches to allow for the gravity feed of the DTD system.The dam serves other purposes as well like maintenance of navigation and waterway transport.The dam is made up of a concrete section, earth fill section, and a lock that can take 1000.0 t vessels.The concrete part of the dam consists of seven spillways equipped with segment gates each 24.5 m in width.The segment gates are equipped with flaps for fine regulation of water discharge.The crest of the spillway is at 68.0 m a.s.l.The earth-fill part of the dam includes three concrete-lined auxiliary spillways each 24.0 m in width.Their spillway crest level is at 76.5 m a.s.l.The segment gates are designed for the 100-year flood, with the tailwater level at 80.35 m a.s.l.The discharge of the Tisa ranges from 80.0 m 3 /s to almost 4000.0m 3 /s.At the dam site, the difference between the upstream and downstream water levels ranges from 0.0 m to 5.0 m for 305 days in a year, while during the rest of the year the segment gates are open resulting in almost identical upstream and downstream water levels.The problem with the Tisa Dam is that high velocities occurring downstream of the dam, when water is discharging over the segment gate flaps, result in large water level fluctuations within the right bank area, which causes erosion of the right bank downstream of and in immediate vicinity of the emergency spillway.In addition to the large-scale bank erosion which is -according to the observations of the dam staff -the most intensive during low tailwater levels, the resulting flow pattern causes extensive river bed erosion at the location of interest, contributing to the instability of the riverbank, and thus affecting stability of the dam.Some of the eroded material is backfilled locally in form of an island in the middle of the affected area (referred to below as "erosion bay"), while the rest is moved downstream by the currents.All this together can lead to sediment displacement at this location, a gradual retraction of the right bank by succeeding bank slips which might reach and undermine the foundation, and finally destroy the geodetic bench installed nearby.The problem is generally caused by bed erosion in combination with bank slides within the right bank area resulting from high downstream flow velocities and water level oscillations.The paper describes in detail the application of a spatial (threedimensional) numerical flow model to highlight the causes of the problem and to find a suitable solution to the bank instability issue.Turbulent flow, complex flow patterns, dischargedependant physical processes, are characteristic to the flow downstream of the dam in combination with suspended and near-bed sediment transport associated with bed erosion and sediment deposition variable in space and time, thus requiring an application of a spatial (3D) mathematical model to simulate the ongoing processes as faithfully as possible in order to predict the outcome.Numerous methods are nowadays available for numerical modelling of flow, sediment transport, and morphological changes in river courses.However, many of the existing flow and sediment models exploit incomplete formulation of sediment processes.Instead of bedload and suspended-sediment transport modelling, some earlier models [1] operate with the total load simplification.Other models ignore bedload transport [2][3][4], changes in bed levels [3,4], or the transport of suspended-sediment [5][6][7].Some currently used models involve an incomplete description of sediment exchange processes [1,[4][5][6][7].It is well known that the method of accounting for the bed and near-bed processes has a significant impact on suspended-sediment transport, which is governed by sediment exchange mechanisms.Besides, the flow field is influenced by bed level changes.The Exner equation is utilized for computing bed-elevation changes in models described in [1,2,5,6,8].The bedloadlayer concept is exploited by Van Rijn [9] and Minh Duc et al. [10] to compute bedload transport, while bed level changes are determined by the total-load approach.Spasojević and Holly [11,12] employed the active-layer and active-stratum method to calculate bed-elevation changes and bedload movement.Olsen [13] used the difference in sediment continuity of bed cells to estimate bed-elevation changes.Finally, just a few of the present flow and sediment models operate with mixtures of sediment or exploit some aspects of the performance of Mathematical modelling of flow and sediment transport -the Tisa dam case study sediment-mixtures.A model for bed and near-bed processes relying on the Exner equation with bedload only, utilizing the bed armouring procedures of the ALLUVIAL model [14], is presented by Kassem and Chaudhry [7].Olsen [13] produced a model for solving the three-dimensional (3D) suspendedsediment transport equation.He applied the near-bed sediment concentration as a boundary condition.A discrepancy in sediment continuity of bed cells was utilized to compute the changes in bed elevation.The change in bed grain size distribution was determined by the budget method.Hung et al. [15] produced a model which follows closely the active-layer and activestratum methodology.It exploits the mathematical formulation and numerical solution found in the MOBED2 model to handle sediment mixtures, as introduced and applied by Spasojević and Holly [11,12].Practical application of the mathematical model for flow, sediment transport, and bed-elevation changes, can be found in [17][18][19][20][21].The current work exploits the Delft3D [16] software for flow pattern solution.Delft3D-Flow solves the Navier Stokes equations for an incompressible fluid, whereby shallow flow is assumed, and the Boussinesq assumptions are applied.The sediment transport and morphology modules support the calculation of both bedload and suspended load transport in the case of non-cohesive sediments based on the sediment mixture, active-layer, and active-stratum approach.

Review and analysis of available data.
scenarios selected for modelling

Analysis of available data for numerical model
The bathymetric data used for setting up the spatial numerical model of the right bank area at the location of interest are provided in the form of absolute coordinates X, Y, Z of points surveyed by echo sounder in 2010 and 2016 (Figure 2). Figure 2 reveals that a detailed hydrographic survey was carried out by echo sounder just for the bank section at risk (for the so-called "erosion bay") to meet requirements of the Delft3D model, whereas the downstream watercourse was surveyed at three cross-sections and two longitudinal sections.A detailed analysis of the results proves that the survey of the right bank at risk meets the model requirements in terms of both resolution and quality.Furthermore, the three additionally surveyed crosssections revealed a relatively uniform bathymetric chart of the downstream watercourse, where the most downstream cross-section is 665.10 m from the dam.An additional hydrographic survey was carried out in the zone of the stilling basin including the scour hole 50.0 m downstream of the dam; therefore, the measured data were used in the model instead of the blueprint status.The design specifies a stilling basin elevation of 65.00 m a.s.l.(Figure 3), while the measurements indicate the presence of mud cover over the concrete layer having an average thickness of 0.50 m.Also, a hump along the longitudinal midline stretches all along the stilling basin and even farther.This hump might be the remains of the cofferdam built during the dam construction, not removed in an appropriate way.Following the acquisition of data on the bathymetry and structure elevations, the final digital elevation model of the river bed in the area under study was obtained by interpolation, as shown in Figure 4.   Initial and boundary conditions in the form of relevant hydrographs and stage hydrographs are provided for the flow module.Since the upstream boundary conditions are set in this case by the dam itself (more precisely by the regulated discharge over the segment gate flaps and below the segment gates), discharge rating curves are defined as a function of discharge, headwater and tailwater levels, calibrated discharge coefficients, flap positions, and gate openings (Figure 5).The discharge rating curves are identical for all spillways and meet the requirements of the model in terms of upstream boundary conditions.Psamological maps obtained by a survey of the water intake area of the Novi Bečej -Banatska Palanka canal (located immediately upstream to the dam) were used in the module for the simulation of suspended sediment transport and morphological changes of the river bed.The input data include the following significant relations: a) a graph showing discharge vs concentration of suspended solids in the Tisa River, b) particle size distribution of the suspended load, c) mean particle size of suspended solids vs discharge of the Tisa River, d) hydrometric measurements corresponding to the cross-section just upstream of the dam performed in 1988, and e) envelopes of grain size curves.The data obtained on sediments provide information on the total suspended load in the form of mass concentrations, mean hydraulic particle size of suspended solids, and mean particle size (d 50 ) of bedload.Since the Delft3D sediment transport module is based on sediment mixtures, it requires the input of characteristic fractions which represent the grain size distribution of the bedload and suspended load.For that, the envelope grain size curve is used referring to the location of the dam itself, obtained by samples taken before construction of the dam.Due to the time and spatial discrepancy of the available sediment data, an additional analysis was carried out to get a better insight into the sediment grain size distribution on the site.By comparison of the envelope grain size curves with relevant psammological measurements carried out just upstream of the dam in 1988, the suspended load and bedload grain size curves were reproduced and adopted for the calculations based on the measured mean hydraulic grain size of the suspended load and the mean bedload particle size.

Scenario selected for modelling
The scenario adopted for investigation based on preliminary analysis of available data required: Creation of a three-dimensional flow model to aid the identification of the critical scenario in terms of hydraulic parameters i.e. to determine which of the dam operation regimes provides the most extreme flow pattern.First, the model is calibrated via a sensitivity analysis using the available data.Next, different regimes of dam operation are analysed based Mathematical modelling of flow and sediment transport -the Tisa dam case study on the number of active spillways, flow over the segment gate flaps or below the segment gates, characteristic discharges, and tailwater levels.Finally, the most demanding scenario for the right bank at risk in terms of flow conditions is adopted.Following the implementation of the flow model, the sediment transport model, having the capacity for dealing with the transport of suspended solids, bedload transport and deformation of the river bed, has been established.It is calibrated based on the available data and tested for the predefined set of selected critical flow scenarios.Finally, the regime of dam operation putting the right bank at the highest risk is identified based on preliminary testing of critical flow and sediment transport scenarios.Thus defined relevant conditions implemented into the model are expanded with the conditions introduced by river training works supposed to provide the solution.The final solution is the optimum one dealing with the problem in its full complexity (technical feasibility, cost-benefit ratio, etc.).

Methodology
Mathematical modelling of flow and sediment transport is used in the framework of this paper as a means for finding the reasons for the left bank erosion, which has taken place downstream to the dam.It has also helped in preparing the proposed solutions to overcome the problem.The Delft3D mathematical model, capable of modelling flow and transport of sediment mixture and morphological changes of the river bed in natural watercourses, has been used.The Delft3D model was applied to the case of the Tisa dam to check the basic hypothesis according to which the erosion bay at the right bank came about through permanent riverbed erosion due to the flow pattern created by flow from the nearby spillway.For that, the investigation targeted two issues: a) discovering the main reasons causing erosion of the right bank depending on the regime of dam operation, b) providing a solution for the described anomaly.At first, the mathematical model of flow, transport of sediment mixture, and river bed deformation, was set up based on available data.Next, the model was calibrated using the sensitivity analysis.Following the calibration, the model was applied to various flow scenarios to see how and to what extent different regimes of dam operation influence bank erosion.Next, according to the results of the previous analysis, the effectiveness of river training works in preventing further enlargement of the erosion bay was checked using the established model.The outcome of the analysis is the solution in form of the river training works that are optimum from both technical and economic points of view.The main feature of the presented analysis is the use of mathematical modelling to reproduce the pre-existing conditions (sediment transport, morphological changes of the river bed) corresponding to the regime of dam operation, and also to model the influences and efficiency of river training works.Having in mind the complexity of the analysed problem in terms of flow, sediment transport and morphological changes of the river bed at the locality of interest due to the presence of the dam itself, the reproduction of physical processes involved is possible only by using a model, either physical or mathematical.Mathematical modelling has been adopted for the analysis of the current case because of its robustness, efficiency, and accuracy.Used as aid, it has proven to be very efficient in reproducing various scenarios of dam operation with and without river training works installed, thus helping to prevent further development of the erosion bay.

Generation of 3D model
The spatial model is generated based on available bathymetry data.Since the dam -having a strictly controlled operation in form of regulated discharges along each spillway (discharge over the segment gate flaps or below the segment gates) -is located at the upstream section of the analysed watercourse, the dam itself sets the upstream boundary conditions.Consequently, a section just upstream of the dam at a distance of about 25.0 m, having a bed elevation of 66.0 m a.s.l., is adopted as the upstream boundary with water level readings corresponding to the registered discharge-head values of segment gates, providing appropriate spillway discharges.The most downstream cross-section surveyed, located 665.0 m downstream of the dam is adopted as the downstream boundary.It has been chosen for two reasons: on the one hand, it is far enough from the right bank at risk, which is why the downstream boundary condition is not expected to influence the conditions at the location of interest.On the other hand, the river bed morphology of the chosen boundary section is already available.
Once the boundary conditions are set, the implementation of the computational grid follows.This primarily means covering the entire physical domain with an appropriate curvilinear grid, whereas computation is carried out in transformed orthogonal coordinates.Due to the complex nature of the problem under study, the so-called ''domain decomposition'' method is used, where various zones of the modelled area are covered by grids of different resolutions.This is appropriate, partly because it enables subsequent local increase of grid resolution resulting in higher accuracy and stability of the model, partly because because the calculation efficiency is increased, i.e. coarser grid is applied to less demanding zones thus reducing computational time.Three different grid resolutions have been used in this model that is applied to three different segments of the modelled domain.The first segment includes the dam itself, the concrete paved plain area downstream to the dam (i.e. the stilling basin), and a 100.0 m long river course to the downstream of the stilling basin.The first segment is covered by a computational grid having 110 x 37 computational points with an average grid cell size of 8.0 m x 2.0 m.The second segment consists of the area downstream of the emergency spillway, to the right of the stilling basin.It is covered by a mesh having 27 x 63 computational points, with grid Matija Stipić, Ljubomir Budinski, Julius Fabian cells of average size 5.0 x 2.0 m.The third segment includes the watercourse section down to the downstream boundary of the model, covered by a computational grid of 30 x 27 computational points resulting in an average cell size of 20. 0 x 10.0 m. in this way, the most downstream area of the model, which has the lowest influence on the flow pattern, is covered by a grid of the lowest resolution, while the first and the second segments are covered by grids of much higher resolution to reproduce the involved flow and sediment transport processes with high quality, as accurately as possible.Moreover, since the model is 3D and requires vertical layering, a total of five vertical computational layers are defined over the entire model domain.Next, the initial test of the flow model, model of sediment transport, and the model of river bed deformation, is performed for relevant scenarios described earlier.

Flow model
The implementation of the existing conditions involves initial test of the flow model, sediment transport model, and the river bed deformation model, in order to establish appropriate values of parameters, particularly the value of the coefficient of turbulent diffusion and dispersion in all three directions, and the value of friction coefficient.Also, due to the changing regime of dam operation, the number of active spillways and the corresponding discharges related to the critical flow scenarios need to be established for the current conditions.This means that the number of spillways in operation, the way of water release (over the flap or under the gate), and the corresponding discharge in the spillway need to be varied to determine the most extreme flow conditions relevant for the computation of sediment transport and morphological changes of the river bed.
Based on available hydrological-hydraulic data and onsite observations, the ten-year minimum water level of 70.56 m a.s.l., taken from the annual records of the Republic Hydrometeorological Service of Serbia, has been adopted for the characteristic tailwater level.This water level was chosen for the most intensive bank slides that were observed at very low tailwater levels, corresponding to the biggest headwatertailwater difference.Preliminary analyses have shown that this interpretation is reasonable especially as the low tailwater could not resist the turbulence stimulating jets coming from active spillways.In addition to the tailwater level, two more parameters are required for definition of the flow model, namely the number of spillways in operation, and the corresponding discharges.The staff members running the operation of the Tisa dam observed that large scale bank slides occur when the seventh spillway (the one closest to the right bank) is in operation and water discharges over the flap of the segment gate, which has been taken as an initial condition to start with.Discharges have been varied in the range corresponding to low tailwater levels, i.e. up to 700.00 m 3 /s.Each particular discharge value was set by adjusting the headwater level, the position of the segment gate, and the corresponding value of the discharge coefficient.The number of spillways in operation has also been varied in parallel with the variation of the discharges in order to pinpoint relevant conditions generating the flow with the highest influence on the river bed and riverbank stability.Accordingly, the first case modelled is the seventh spillway in operation only with a minimum discharge of 26.82 m 3 /s over the segment gate flap at the tailwater level of 70.56 m a.s.l.The sensitivity analysis used for initial model tests resulted in the Manning's coefficient of n = 0.018 s/m 1/3 , whereas v t = 0.65 m 2 /s was adopted for horizontal turbulent viscosity.These results indicate an intensive flow over the segment gate flap and a rotational flow (with an island in its centre) in the erosion bay.Similar results were obtained by physical model tests carried out by the Department for Hydraulics of the Jaroslav Černi Institute for the Development of Water Resources.Only the visual presentation of the measured flow patterns was available to the authors.For the implemented discharge of 26.82 m 3 /s it is typical for the studied area that the mean velocity is pretty low (around 0.001 m/s) with a considerable increase in velocity (up to 0.080 m/s) at several locations: between the main flow and rotational flow, as well as at two locations between the island and the right bank.However, velocity is peaking at a position where the rotational flow releases water into the main stream, resulting in the flow rate of 0.28 m/s.In the first segment, at the point where the jet is leaving the spillway, the velocity is as high as 0.3 m/s.To the downstream along the stilling basin, it decreases to about 0.15 m/s, while in the remaining part of the stilling basin the velocity is about 0.0001 m/s, meaning that the water is almost still.Velocities are quite uniform in the third segment, amounting to approximately 0.03 m/s.Since the resulting flow pattern is entirely in line with the previous observations and the available data, the specified values of parameters -Manning's coefficient and the coefficients of horizontal turbulent viscosity -are taken as characteristic for all cases.
In the presented case, only the seventh spillway was in operation.Further cases were investigated additionally to see if the simultaneous operation of more spillways influences degradation processes within the erosion bay at the right bank.At low tailwater level -which was identified as a relevant condition for bank erosion -in ordinary circumstances, two or four spillways are likely to be simultaneously active in symmetrical arrangement as much as possible to avoid concentration of high flow velocities.Therefore, to get an insight into the impact of multi spillway operation on the scouring process, two additional cases were considered: parallel operation of two spillways (the sixth and the seventh) was considered, as well as the case of four active spillways (the second, the fourth, the sixth and the seventh).It was observed from simulated velocity patterns that the flow pattern within the erosion bay did not change in terms of the spatial arrangement of velocity vectors, i.e., the earlier determined rotational flow around the island, with a secondary eddy between the island and the jet coming from the seventh spillway, is present in these two cases as well.On the other hand, Mathematical modelling of flow and sediment transport -the Tisa dam case study there is an apparent increase in velocity magnitudes by about 5 ÷ 10 % over the entire area compared to the previous condition.
In general, the activation of additional spillways does not have a significant impact on the flow pattern within the erosion bay in the sense of spatial distribution and intensity.Accordingly, the case of four active spillways (the second, the fourth, the sixth, and the seventh) was identified as critical and representative, and was therefore selected for further investigation.
The impact of discharge on the flow pattern within the erosion bay was considered after analysis of the number of active spillways influencing erosion.For this, a set of three flow rates 60.0 m 3 /s, 90.0 m 3 /s, and 112.0 m 3 /s per spillway was considered with flow patterns produced for all discharges.The adopted flow rates correspond to low tailwater levels since low tailwater is considered a critical and relevant condition for the erosion of the right bank.Since the flow over the segment gate flaps was considered, the discharge in the Delft3D model was set by providing the value of the flow head and the discharge coefficient.As in the previous cases, the model reached the steady-state condition in about six hours, which was verified by checking the change of water levels and velocity readings over time in the entire calculation domain.A comparison of flow patterns and intensities of velocities is provided in Figure 6.It should be noted for the sake of clarity that all velocity vectors will be further shown herein only as depth-averaged (although the model generates three-dimensional velocity vectors).
The findings indicate that the spatial arrangement of velocity vectors remains unchanged over the entire domain, i.e. that discharge variations do not affect streamlines in terms of spatial orientation of velocity vectors in the stilling basin and the erosion bay at the right bank.The rotational flow around the island and the secondary eddy between the jet and the island are still the same as earlier in terms of vector directions and locations.However, the intensities of velocity vectors increased with an increase in discharge, especially in the earlier identified locations.This means that maximum velocities in the erosion bay are around 0.32 m/s at the discharge of 60.0 m 3 /s, 0.37 m/s at the discharge of 90.0 m 3 /s, and 0.42 m/s at the discharge of 112.0 m 3 /s.Having in mind the bathymetry of the erosion bay, the flow pattern is quite reasonable.Most of the flow is concentrated along the deepest middle section (having the lowest resistance), whereas shallow sections remain passive regardless of changes in discharge.The secondary eddy, the streamlines of which are shown in Figure 6a, significantly contributes to this trend.In general, the presented discharge variation analysis proves that an increase in discharge ranging between 60.0 ÷ 112.0 m 3 /s has no impact on spatial distribution of velocity vectors, that is, higher discharges do not affect the basic flow pattern within the erosion bay at the right bank.However, higher discharges increase velocities along the primary flow path by about 10% per each increase of 30.0 m 3 /s per spillway.Releasing water below the segment gates was the last numerical experiment related to the current conditions (without river training works applied).Since the trend of magnitudes and spatial distribution of velocity vectors concerning the number of spillways in operation and discharge variations has already been determined by the previous numerical experiments, the release of water below the segment gates with a discharge of 100.00 m 3 /s was addressed in the following experiment.All other conditions -such as the number and order of spillways in operation, the initial and boundary conditions, the tailwater level at 70.56 m a.s.l.-were adopted from the previous analyses.The resulting flow pattern for the case of flow below the segment gates is shown in Figure 7.Comparison of this scenario with the case of flow over the segment gate flaps indicates that in both scenarios the flow pattern retains its characteristic rotational flow around the island located in the middle of the erosion bay, and the secondary eddy between the jet coming from the seventh spillway is still present.The computed magnitudes of velocity vectors indicate that the maximum velocity along the primary flow path within the erosion bay is now 0.32 m/s.This is also reasonable having in mind the method of water release from the spillway and the physical processes involved.
In general, the analysed cases have shown that the spatial distribution of velocity vectors in the erosion bay at the right bank (see the flow patterns) is rather insensitive to any variation in discharge, way of water release through the spillways, or to the activation of further spillways in addition to the seventh spillway.The basic flow pattern stoutly remains the same, characterized by rotational flow around the island and a secondary eddy between the island and the jet coming from the seventh spillway.On the other hand, velocity magnitudes are directly affected by discharge variations, the way of water release, and the number of spillways in operation.

Modelling sediment transport and river bed deformations
The second part of the analysis related to the existing conditions is the modelling of sediment transport (suspended solids and bedload) and river bed deformations, depending on the already defined flow patterns.Although sediment transport computations were accomplished for all the previously analysed flow scenarios, for the sake of clarity, the following discussion will be limited to two critical scenarios only: flow over the segment gate flaps of four active spillways with a discharge of 112.00 m 3 /s, and flow below the segment gates of four active spillways with a discharge of 100.00 m 3 /s.As mentioned earlier, the Delft3D program package is equipped with a computational module for the calculation of suspended load transport, bedload transport, and river bed deformation resulting from sediment transport, based on sediment mixture.This means that sediment is not accounted for by its median grain size d 50 -which is a characteristic of less sophisticated numerical models -but grain size distribution is exploited, which is much more revealing for sediment transport modelling.This approach is more detailed, provides a more realistic description of sediment transport and morphological changes of the river bed.Given this, sediments are analysed and modelled by sediment fractions in the adopted modelling method, i.e. the numerical model requires grain size distribution for both the suspended load and bedload.In the Delft3D model, the grain size distribution of sediments is entered in form of an arbitrary number of fraction ranges, where a median grain size is computed for each defined fraction range as the geometric mean of the upper and lower boundary of the fraction range.Three fraction ranges were adopted following the analysis of grain size distribution in the suspension and on the river bed: fraction I1 ranges between 0.0003÷0.017mm, fraction I2 ranges between 0.017÷0.20 mm and fraction I3 ranges between 0.2÷1.00mm.It is important to notice that fraction I1 is present in suspension only, fraction I3 is present in bed only, whereas fraction I2 is found in both domains.As with the flow pattern model, initial and boundary conditions need to be set following the definition of fraction ranges for the module of sediment transport and river bed deformations.The average value of total suspended load concentration for flows of 1000.0 m 3 /s or less is estimated at 0.025 kg/m 3 based on the graph of suspended load concentration vs flow that is  Matija Stipić, Ljubomir Budinski, Julius Fabian adopted for sediment modelling.Since this value indicates the total concentration while the model uses sediment mixture ranges, the total suspended load concentration needs to be distributed between fractions I1, I2, and I3.Finally, the following suspended load concentrations defined by fraction ranges were obtained: I1 = 0.010 kg/m 3 , I2 = 0.015 kg/m 3 , and I3 = 0.0.It was assumed that these concentrations defined by fraction ranges were already present in the suspension, thus these values were both used as initial and boundary conditions for the upstream boundary.
On the other hand, the bedload mixture was defined in the form of proportions of particular fraction ranges as I1 = 0.0 %, I2 = 30 %, and I3 = 70 %.Following the boundary conditions that were set and the already chosen critical flow scenarios, the computation of suspended load and bedload transport was initiated.At last, the river bed deformation in the form of cumulative erosion/deposition was modelled for a single day period.The resulting distribution of a particular fraction range in the third layer at the bay, and the cumulative erosion/deposition for the case of water release over the segment gate flaps and below the segment gates (fraction I3 is completely absent from the suspension, hence the figure showing I3 distribution is omitted) are shown in Figure 9 and Figure 10, respectively.The results in both scenarios indicate that morphological processes have reached a steady-state, i.e., that equilibrium has been achieved in sediment transport processes within the specified period of modelling.An intensive flow in the stilling basin resulting from the release of water below the segment gates causes large scale erosion in this area promoting sediment intake from the bed of the stilling basin into the suspension.Graphs with concentrations and cumulative erosion/deposition, presented in Figure 10, clearly show that high load concentrations in the stilling basin are caused solely by bedload intake due to intensive flow from the spillways, while the so taken sediment is transported further downstream into the third segment.Also, the third fraction I3 is partially activated, which points to a very intensive flow in this area.Concerning the erosion bay at the right bank, due to bedload intake, increased concentrations of up to 0.014 kg/m 3 for fraction I1 and 0.005 kg/m 3 for fraction I2 in a stretch are indicated along the segment, while 0.0001 kg/m 3 is characteristic to the remaining part of the segment.This confirms the fact that, even if a higher sediment content from the main flow enters the erosion bay, it is transported by the flow like a high concentration stretch to the exit from this area.The third fraction I3 is not activated in the erosion bay.Concerning erosion/deposition processes within the erosion bay, locations identified as critical in the previous scenarios with the flow over segment gate flaps, are also effective in the case of flow under the segment gates, but of slightly lower intensity.
In general, the results from the sediment transport module indicate that the model of sediment transport processes and morphological changes is accomplished in the Delft3D environment for the studied location.Therefore, the entire sediment module truly reproduces spatial and temporal conditions of sediment transport and morphological changes of the river bed.Accordingly, it can be stated that both modules are properly calibrated at this stage and, as such, that they are capable of simulating potential impacts of any river training structure on hydraulic, morphological, and sediment transport conditions within the domain of interest.

Solution with river training works
The main criterion for assessing variou s potential solutions using river training works is their efficiency to reach the goal, i.e. to eliminate unwanted impacts of the flow downstream to the seventh spillway putting at risk the erosion bay and the right bank.Also, the proposed solution needs to reduce or eliminate permanent erosion of the river bed on the site of interest and prevent the transport of sediments from the erosion bay to the main flow.Solution options include analysis of the effect of (1) a longitudinal dike (parallel to the flow) on flow and morphological conditions at the right bank, proposed by the Public Water Management Company Vode Vojvodine, and (2) a groyne on the right bank in several variants.The longitudinal dike was investigated for the scenario with flow below the segment gates with the discharge of 100.0 m 3 /s at tailwater level 70.56 m a.s.l.The groyne was investigated also for flow below the segment gates for several discharges and tailwater levels: 100.0 m 3 /s at tailwater level of 70.56 m a.s.l., 150.0 m 3 /s at tailwater level of 72.0 m a.s.l., and 200.0 m 3 /s at tailwater level of 74.0 m a.s.l.

Option I -longitudinal dike
The first option considered as a possible on-site intervention for the elimination of adverse impacts of flow and sediment transport which put the stability of the right bank at risk is construction of a longitudinal dike (parallel to the main stream) between the island and the right bank.The proposed length of the structure is 140.0 m, whereas the crest level is set to 73.0 m a.s.l.Numerical simulation results are given in the form of depth-averaged velocities, the spatial distribution of sediment load concentrations, and cumulative erosion/deposition charts, as shown in Figure 11.
As mentioned earlier, the effect of longitudinal dike is considered for the scenario with flow below the segment gates at the discharge of 100.0 m 3 /s and for the tailwater level of 70.56 m.The first thing to note regarding the spatial distribution of flow is an almost complete elimination of flow on both sides of the island.Since the earlier present rotational flow around the island within the erosion bay is blocked, the flow is now passing by the island attacking more heavily the left side of the island oriented towards the main flow, and the downstream section of the erosion bay.Concerning sediment transport and morphological changes of the river bed, the model has shown that the sediment initially found in suspension in the erosion bay has settled down right there.This suggests that any mass Mathematical modelling of flow and sediment transport -the Tisa dam case study of soil, which eventually got into the erosion bay by potential bank slides triggered by pore pressure variations due to water level oscillations, would remain in the erosion bay area and would facilitate silting up of the scour hole.Concerning the deformation of the river bed, slightly higher erosion/deposition can be seen by the island on the side oriented towards the main flow and in the downstream section of the erosion bay at the point of local circulation.This supports the above mentioned observation that since the flow cannot uphold the current condition, it is now redirected to these locations with increased velocities, resulting in higher erosion and deposition rates ranging from -0.02 m to 0.02 m.However, a complete lack of deformation is observed in the remaining portion of the erosion bay, which was the primary purpose of the proposed structure.Moreover, the model demonstrates that the implementation of a longitudinal dike does not cause any safety risk to the downstream part of the river bank and to the watercourse.The second option considered is the installation of a groyne of a significantly smaller size compared to that of of the longitudinal dike.The groyne 36.0 m in length, with the crest at 73.0 m a.s.l., is considered as the first sub-variant presenting hydraulic conditions identical to those for the longitudinal dike.The crest elevation was selected following a preliminary hydrological analysis of tailwater levels.The selected crest level is equal to the average annual tailwater level.Numerical simulation results are given in Figure 12 as depth-averaged velocities, spatial distribution of sediment concentrations, and cumulative erosion/deposition.It can be seen in Figure 12 that the groyne is completely connected with the island, eliminating the flow As with the previous option, it can be seen from the resulting flow pattern that there is a complete lack of flow in the area adjacent to the bank, that is, along the sides of the groyne, thus indicating the same effect as the option with the longitudinal dike.The flow pattern in the remaining part of the erosion bay is similar to that in the previous option, i.e., higher velocities of the order of magnitude of approximately 0.20 m/s occur just along the side of the island oriented towards the main stream and in the downstream portion of the erosion bay.The remaining part of the modelled domain is unchanged, which is true for both the sediment transport and the river bed deformations.The first fraction I1 gradually settles around the groyne indicating ''still'' water in this zone, whereas the second I2 and the third I3 fractions show an identical pattern to that in the option with the longitudinal dike.The same applies to the cumulative erosion/deposition, i.e., these processes are localized at identical locations as in the previous option, which is the direct consequence of the modified flow pattern at this location.A groyne with the crest level kept at 73.0 m a.s.l.but with the tailwater level increased to 72.0 m a.s.l. and the flow rate increased to 150.0 m 3 /s per spillway is considered as the second sub-variant.This results in an enhanced flow pattern, where the groyne crest is above the water table while the island with a maximum level of 71.50 m a.s.l. is completely submerged.The results of the numerical simulations are given in Figure 13 in form of depth-averaged velocities, spatial distribution of Mathematical modelling of flow and sediment transport -the Tisa dam case study sediment concentrations, and cumulative erosion/deposition.It can be seen from the resulting flow pattern that the water is overtopping the island.The flow is exaggerated over the submerged island reaching velocities of up to 0.30 m/s.In general, the island overtopping caused the acceleration and enlargement of the secondary eddy initially formed between the island and the main flow.These conditions, together with the intensive main flow from the seventh spillway, might put this location at risk.Sediment transport and changes in river bed morphology follow the new flow pattern resulting in large scale erosion at the island.Due to the new conditions which disrupt the former balance between the flow and the morphological processes, an increased erosion of fractions I1 and I2 takes place here followed by the intake of eroded sediments into the main stream.This practically means that these hydraulic conditions lead to erosion of the island, which will last until a new balance is reached in hydraulic/morphological sense in this location.In the remaining part of the modelled area, conditions are further intensified by an increase in flow, but without significant changes in spatial distribution.
By the third sub-variant, the case of complete submersion of both the island and the groyne is investigated.The crest level of the groyne is reduced to 71.0 m a.s.l.To meet the described conditions, the tailwater level was set to 74.0 m a.s.l., and the flow rate was increased to 200.0 m 3 /s per spillway.The increase of discharges in the sub-variants is meant to increase flow velocities as well as sediment transport processes, which, in turn, should indicate the probable critical locations within the erosion bay.In terms of discharges and water levels,  Mathematical modelling of flow and sediment transport -the Tisa dam case study combinations of low discharges with low tailwater levels and high discharges with high tailwater levels are investigated.The regime of dam operation is followed.Theoretically, these are the most extreme scenarios.Combinations of low tailwater levels and high discharges were not analysed since their occurrence is not probable in practice.In full compliance with the project assignment, cases with a high total discharge of 800.0 m 3 /s were analysed in this set of numerical experiments.The results are given in Figure 14 in form of depth-averaged velocities, spatial distribution of sediment concentrations, and cumulative erosion/deposition.The resulting flow pattern, with both the groyne and the island submerged, point to a flow pattern similar to the one with low water levels without a groyne, but with the centre of the eddy shifted more towards the main flow, i.e., towards the jet coming from the seventh spillway.The flow enters the erosion bay from the downstream end, sets up a large eddy, and leaves the erosion bay at the upstream end.Although the flow pattern is similar to the one without the groyne as to spatial orientation of velocity vectors, a higher discharge through the spillways caused a significant increase in velocities.Maximum velocities of around 0.80 m/s are observed by the river bank at the downstream end of the erosion bay.A bit lower but still high velocity value of around 0.70 m/s is recorded in the location where the flow crosses the groyne as well as ten meters to the downstream.In the remaining part of the erosion bay, velocities are around 0.40 m/s.Based on the provided flow pattern it can be concluded that for tailwater levels at which the groyne and the island are completely submerged, the flow gains intensity within the erosion bay since nothing blocks the flow anymore above the crest of the groyne.The local flow velocity over the groyne is high, indicating that the height of the groyne plays a significant role in affecting the flow pattern within the erosion bay.
As to sediment loads, a higher concentration of the first two fractions can be observed in the erosion bay, while the third fraction is still absent from the suspension.A higher concentration of the first fraction I1 can be observed along the bank pointing to an intake of sediment from the main stream, and to an intake of sediment from the bottom due to intensified erosion.The average concentration in this area is around 0.030 kg/m 3 , while it has stabilized at 0.022 kg/m 3 in the remaining area.Large scale erosion of the bank is indicated by the spatial distribution of the second fraction I2, which also produced a stretch of increased concentration of around 0.025 kg/m 3 along the bank, whereas the concentration is considerably lower, around 0.0075 kg/m 3 , in the remaining area of the erosion bay.The stretch of increased concentration corresponds to the graph of cumulative erosion/deposition, indicating an increased river bed deformations adjacent to the bank, averaging at around 0.30 m.This inevitably indicates a risk to the bank stability not only at low tailwater levels and low discharges but also at high tailwater levels and relatively high discharges as well.

Discussion
The analysis of the current flow patterns at the location of the dam shows that the spatial orientation of velocity vectors within the erosion bay is neither sensitive to the discharge nor to the way in which water is released through the spillways (over the segment gate flaps or below the segment gates).However, flow intensity is sensitive to discharge and to the way water is released, while activation of spillways in addition to the seventh spillway does not change the intensity of the resulting velocity vectors significantly.Concerning the computation of river bed deformation, several critical points have been identified as to erosion/deposition along the bank and along the main stream, which point to the instability of both the river bed and the riverbank at this location.The implementation of a longitudinal dike demonstrates that the proposed longitudinal dike can eliminate adverse impacts of the flow pattern on the right bank, thus contributing to the stability of the river bed and the riverbank within the erosion bay.However, given the huge volume of the longitudinal dike having a total length of 140.0 m, it is obvious that alternative solutions need to be analysed to find a less costly solution.The second structure proposed is a groyne perpendicular to the right bank connected both to the bank and to the island.In comparison with the longitudinal dike, the primary advantage of the groyne having a total length of 36.0 m is its smaller size and its location.In general, this structure offers the same benefits concerning the flow pattern, sediment transport, and river bed deformation as the longitudinal dike does.Since the groyne is connected both to the bank and to the island, the resulting structure eliminates flow and morphological processes within the erosion bay.Besides, the layout of the structure promotes the natural silt up of the area around the groyne, which, in turn, reduces flow and increases stability not only of the bank but also of the river bed at the location of interest.The analyses lead to the conclusion that this type of structure is more effective from the technical and economic standpoints, compared to the longitudinal dike.

Conclusions
Using a 3D numerical flow and sediment transport model proved to be a powerful tool for identification of the reasons for riverbank erosion resulting in a slowly but permanently increasing erosion bay downstream of the Tisa Dam.Furthermore, it has aided the selection of the most suitable solution by predicting hydraulic effects of a longitudinal dike and a groyne on the already existing erosion bay.Finally, the groyne has been selected as a solution to the problem due to its favourable hydraulic effects and low cost.

Figure 1 .
Figure 1.DTD system map with location of Tisa Dam at Novi Bečej

Figure 2 .
Figure 2. Bathymetric chart of Tisa Dam location

Figure 3 .
Figure 3. a) Baseplot b) longitudinal section of the dam

Figure 4 .
Figure 4. Digital elevation model of the river bed

Figure 5 .
Figure 5. a) Discharge (Q) -water level (Zs) relationship for flow over flaps; b) discharge coefficient (µ) depending on position of the segment gate (a-gate opening, H-upstream depth)

Figure 7 .
Figure 7. Flow pattern in case of water release below segment gates with a discharge of 100.0 m 3 /s at tailwater level 70.56 m a.s.l.

Figure 9 .Figure 10 .
Figure 9. Distribution of suspended load concentration for fractions a) I1, b) I2, c) I3 and d) distribution of cumulative erosion/deposition for the scenario of discharge over the segment gate flaps

Figure 11 .
Figure 11.a) The resulting flow pattern, b) detailed flow pattern of the second segment, c) suspended load concentration at fraction I1, d) cumulative erosion/deposition for the longitudinal dike in place level at 73.0 m a.s.l. and tailwater level at 72.0 m a.s.l.

Figure 12 .
Figure 12. a) Resulting flow pattern, b) detailed flow pattern of the second segment, c) distribution of suspended load concentration at fraction I1, d) cumulative erosion/deposition for groyne in place with crest elevation at 73.0 m a.s.l and tailwater level at 70.56 m a.s.l.

Figure 13 .
Figure 13.a) Resulting flow pattern, b) detailed flow pattern of the second segment, c) suspended load concentration at fraction I1, d) cumulative erosion/deposition for the groyne in place with crest elevation at 73.0 m a.s.l. and tailwater level at 72.0 m a.s.l.

Figure 14 .
Figure 14.a) Resulting flow pattern, b) detailed flow pattern of the second segment, c) suspended load concentration at fraction I1, d) suspended load concentration at fraction I2, e) cumulative erosion/deposition for the groyne in place with crest at 71.0 m a.s.l. and tailwater level at 74.0 m a.s.l.