Estimation of Changes in Sediment Transport along the Free-Flowing Middle Danube River Reach

Abstract

The subject of this study is an approximately 300 km long Middle Danube River reach that spans from river kilometer (rkm) 1581 in Hungary to (rkm) 1255 in Serbia. The observed drying of floodplains in Hungary some thirty years ago initiated the hydrological studies. However, problems related to the navigation route maintenance of the Danube River and those in the water supply of irrigation and drainage canal networks are now present in the whole free-flowing middle Danube region. The study aims at investigating the correlation between the observed water level decrease and recorded incision of the river bed at gauging stations and the indirect estimation of the long-term sediment transport along the sand-bed reach based on the surveillance cross-sections’ data collected during regular monitoring surveys on the navigable Danube. It starts with hydrological analyses of the 70-year-long time series of water level and discharge yearly data and continues with morphological and correlation studies. It ends with the estimation of sediment transport. The decreasing trend in water levels and the increasing trend in cross-sectional areas are persistent. There is a linear correlation between the two. Longitudinal changes in sediment transport indicate the existence of both degrading and aggrading riverbed reaches.

1. Introduction

Large alluvial rivers carry large discharges. As such, they are valuable water resources that are used for multiple purposes (e.g., drinking and industrial water supplies, irrigation, electricity production in hydropower plants, and navigation). At the same time, they can carry significant amounts of sediments that play a role in shaping the riverbed and its planform. Since water is the carrying fluid, sediment dynamics depend on hydrodynamics, which are now modified by different human interventions and uses [1]. In the long run, morphological changes can significantly affect the hydrological cycle of surface and groundwaters through permanent decreases or increases in water levels in the channel and the groundwater table, which, in turn, can produce various imbalances in the environment [2]. This is why regularly monitoring water levels (and/or discharges), on the one hand, and river sediment, on the other, are essentially important. However, some large rivers, like the Danube River, run through several countries, and along some stretches, they are boundary rivers. While the fact that they are transboundary rivers has advantages in the flow of goods and services, the fact that they are also boundary rivers may cause problems in crisis times, especially regarding sediment monitoring. This is the case with the downstream part of the free-flowing Middle Danube River reach (from Bezdan to Novi Sad in Serbia, a part of which is a boundary with Croatia). There has been no sediment monitoring at all along this stretch of the river, not only during the crisis in the 1990ies [3]. The situation is similar to bed load monitoring in the Hungarian reach, where the bed load is sampled only in the gravel bed stretches at the most upstream gauging stations (GSs) of the country, while data from suspended sediment sampling, which is more frequent than bed load sampling, are of questionable quality because (1) the sample is usually not related to flood events when the suspended sediment concentrations are the highest, (2) the cross-section is not usually surveyed during the measurement campaign and (3) the discharge is not usually measured simultaneously with suspended sediment sampling. In contrast, water level monitoring was not and is not susceptible to the crisis to the same extent as the sediment monitoring, and there are only a few gaps in the stage time series.

Owing to the relationship between changes in the water level (stages) and morphological changes in the riverbed (either the incision or aggradation of the channel), trends in the time series of the stage data might be an indicator of trends in riverbed changes. Riverbed deformation results from an imbalance in the sediment transport capacity of the flow and sediment supply. Since water stage monitoring is simpler and cheaper than bathymetric surveys and sediment monitoring, analyzing its trends is a convenient initial step in the assessment of the average multiannual change in sediment transport.

The first sign which initiated analyses of long-term trends in the water stage and discharge data in the free-flowing Middle Danube River reach was an alarming drying of floodplain forests in the Gemenc area in Hungary in the late 1980ies and early 1990ies. The first analyses covered the time series from 1901 to 1992 at 17 GS along the Hungarian reach and 8 GS downstream of Budapest, including Dunaújváros at river kilometer (rkm) 1581, Dunaföldvár (rkm 1560), Paks (rkm 1531), Baja (rkm 1479) and Mohács (rkm 1447) [4]. The results have shown a decreasing trend of water levels at all these gauging stations. However, there was no attempt to correlate it with the suspected incision of the riverbed, most probably initiated by river training activities, mostly carried out in the second half of the 19th century and the first half of the 20th. Recently, Tamás et al. [5]. extended hydrological analyses to the entire free-flowing Middle Danube River reach (from Dunaújváros—rkm 1581 in Hungary to Novi Sad—rkm 1255 in Serbia) following the reports on water intake problems at different locations, out of which the Bezdan pumping station which supplies the Danube-Tisza-Danube irrigation and drainage canal network in Serbia [6]; and decreasing inundation frequencies in the nature reserve Kopački rit on the right bank floodplain in Croatia [7] were outstanding. The analysis covered a 70-year-long period from 1951 to 2019. During the first 30 years of this period, many river training works were completed along this reach to facilitate navigation, thus contributing to setting up one of the ten Transeuropean transportation corridors—the Danube-Rhine corridor. Despite the obligation of all Danubian countries to regularly monitor the river channel for operational purposes, there were only a few attempts to analyse the sediment budget along the river e.g., [8] until the Interreg project “Danube sediment” [9]. The project focussed on the adjustment of the sediment monitoring methodologies and the data in the Danube River basin.

This paper aims to present a comprehensive study of the sand bed free-flowing Middle Danube River reach, with hydrological and morphological analyses and the indirect estimation of the average multiannual sediment transport along the study reach. Given the fact that the suspended sediment monitoring on the Hungarian sub-reach is relatively scarce in space and time [3] and that there was no sediment monitoring at all on the Serbian side [3], the estimation of the sediment transport is, for now, based on analyses of the bathymetric data of surveillance cross-sections from the Directorate for Inland Waterways in Serbia and the Lower Danube-Valley Water Authority in Hungary. This procedure is standard in estimating reservoir sedimentation processes, e.g., the reservoir sedimentation at Shumburit Dam in Ethiopia [10]. A similar analysis was performed along the gravel bed reach of the Danube River between rkm 1800 and rkm 1826 for the period 1992–2001 by Rákóczi and Sass [11]. Martin and Ham [12] demonstrated the advantages of a morphological method based on available field surveys in evaluating the performance of different bed load transport formulae over larger spatial and temporal scales (decadal-scale) on a 65 km long gravel-bed reach of the Fraser River. Tang et al. [13] used bathymetric data collected in two years bounding 30 years (1977 and 2006) to analyze riverbed deformation, i.e., cumulative volume change in the channel of the largest Mississippi Delta distributary—the Atchafalaya River. Analysis of the bathymetric data revealed, among other facts, a discrepancy between the transport capacity of the flow and the inflow of sediments [13].

Sediment transport can also be successfully estimated from dune-tracking methods based on bathymetric surveys of channel longitudinal profiles recorded with a certain time shift [14,15]. Recently, Gaeuman and Jacobson [16] suggested a modified method based on routing local erosion and deposition volumes. It is applicable to channels with a complex bedform morphology when it is difficult to delineate individual bedforms. However, the conventional dune-tracking method and its modifications require on-purpose bathymetric surveys. This is why it is not possible to use these methods for the estimation of sediment transport in less monitored rivers like the free-flowing Danube River reach, as on-purpose surveys have not been performed there.

The research questions to which the answers are sought in this paper are as follows:

(1)

Are the previously observed decreasing trends in water levels stagnating? Did they revise, or are they continuing after 30 years from the initial studies?

(2)

Is there a correlation between the observed minimum water level and the recorded incision of the riverbed at gauging stations?

(3)

What are the long-term sediment transport characteristics along the investigated reach?

The paper is organized as follows. After the description of the study reach, the methods of hydrological and morphological analyses are presented, and the procedure for the estimation of the sediment transport is described. These are followed by a presentation of the results of homogeneity tests and trend analyses of water level time series and discharges, correlation analyses between the cross-sectional area change at GSs and the corresponding lowering of the water level, and the estimation of the average multiannual sediment transport within the sand bed reach with no tributaries. The relationships between all three phenomena are discussed in the fourth section. The key conclusions of the study are provided in the last section. They include an answer showing that the time series of water levels are inhomogeneous along the entire reach while the time series of discharges are homogeneous. It is found that the trend of lowering water levels and the deepening of the riverbed continues at the same rate and that there is a linear correlation between minimum water levels at GSs and percentage increases in the corresponding cross-sectional areas. At the upstream end of the sand bed reach, the 25-year sediment transport rate is positive, indicating the aggradation of the riverbed. In the remaining 86% of the investigated reaches, the 25-year sediment transport rate is negative, indicating the degradation of the riverbed.

2. Description of the Study Reach

The subject of this study is an approximately 300 km long reach spanning from rkm 1581 in Dunaújváros, Hungary, to rkm 1255 in Novi Sad, Serbia (Figure 1), and it belongs to the Middle Danube. The reach is situated between the Gabcikovo (rkm ~1800) and Iron Gate I (~rkm 1000) dams and hydropower plants. The Danube River is alluvial downstream of ~rkm 1600, although the transition from gravel bed to sand bed is around Paks in Hungary (rkm 1531). Since the Danube is a free-flowing river along this reach, in the course of creating the Danube waterway, now part of the European TEN-T network, there was extensive river-training work along the reach in the past. This included cutting off mostly overdeveloped meander bends and the installation of different stone structures after the 1870s (e.g., groynes and longitudinal stone toes). These river-training measures were known for increasing flow energy and sediment transport capacities, and they resulted in deepening erosion in the riverbed. Additionally, there were intensive dredging activities in the upstream part of the study reach between Dunaújváros (rkm 1581) and Paks (rkm 1531) in the second half of the 20th century until the beginning of the 2000s, which provoked further morphological changes in the river channel.

To the best of the authors’ knowledge, no comprehensive analysis of the changes in the water regime or, at the same time, the morphological parameters of the entire free-flowing alluvial reach of the middle Danube River has been performed to date. To fill this gap, available time series of the water level and discharge data for the entire free-flowing alluvial reach of the Danube River and the morphological evolution of cross-sections at corresponding gauging stations (GSs) are analyzed in the present study. Moreover, the digital database of surveillance cross-sections from the state institutions authorized for hydrographic surveys along the international waterway is used to assess variations in the long-term sediment transport along the reach.

3. Materials and Methods

3.1. Methods of Statistical Analysis of Hydrological Data

In the first phase of the hydrological study, analyses were performed using time series of water levels and discharges for the 70-year-long period between 1950 and 2019 for the entire free-flowing alluvial reach of the River Danube downstream of Budapest, Hungary, and upstream of Slankamen, Serbia (~300 km). This was considered necessary, as the older studies published on this topic used different periods, older datasets, and different/shorter river reaches, e.g., in [17]. The analysis included the official water level and discharge data of the Hungarian Hydrological Forecasting Service and the Hydrometeorological Service of the Republic of Serbia for a total of ten gauging stations (Figure 1b): Dunaújváros (rkm 1581), Dunaföldvár (rkm 1560), Paks (rkm 1531), Baja (rkm 1479) and Mohács (rkm 1447) in Hungary, and Bezdan (rkm 1425), Apatin (rkm 1402), Bogojevo (rkm 1367), Bačka Palanka (rkm 1299), and Novi Sad (rkm 1255) in Serbia, similar to Tamás et al. [5].

The first step in this analysis was a test of the homogeneity of the time series. Computations were made in XLStat (Addinsoft, 2023) [18]. A two-tailed, two-sample Kolmogorov–Smirnov test at a confidence level of α = 0.05 was used. The homogeneity test was first performed for values of the yearly characteristics of water levels and discharges for all stations. Consequently, the trend analysis was performed for the annual characteristic values to find trend parameters. Given the fact that there are no large tributaries between the upstream end of the reach and Apatin, it is presumed that discharges upstream of the Drava River confluence, which is downstream of Apatin, do not change significantly. That is why annual discharge data were analyzed only at two gauging stations, Baja and Bezdan, which are located upstream of the confluence of the Drava River. The discharges for the two stations were also tested for normality using XLStat [19].

For further analyses, as the intention was to connect water level changes with morphology, the two stations at the upstream end of the reach (Dunaújváros and Dunaföldvár) were excluded, as the bed material of the river there contains a relatively large portion of gravel, and according to previous studies [20] dredging played a role in the lowering of the water levels on this part of the reach. The two stations at the downstream end of the reach (Bačka Palanka and Novi Sad) were also excluded as the backwater effect and sediment deposition upstream of the Iron Gate I Dam (~rkm 1000) can influence water levels there [14]. The trend of the yearly characteristic values of water level data at the six remaining stations (Paks, Baja, Mohács, Bezdan, Apatin, and Bogojevo) was examined. All trend analyses were performed using linear regression with the least squares (LSQ) method in MS Excel.

3.2. Bathymetric Data Collection

Bathymetric data that were used as a basis for the estimation of the sediment transport were obtained from the databases of the Directorate for Inland Waterways—Plovput in Serbia, the former Environmental and Water Management Research Institute (VITUKI), and the Lower-Danube-Valley Water Directorate (LDVWD) in Hungary. These included 23 surveillance cross-sections in Hungary at an average distance of approximately 2000 m and 46 surveillance cross-sections in Serbia at an average distance of approximately 1000 m. All cross-sections were distributed along the sand bed reach, spanning from Baja (rkm 1479 + 000) to the proximity of the confluence of the Drava and Danube Rivers (rkm 1387 + 200) (Figure 2).

Plovput is the state institution responsible for hydrographic surveys on international inland navigable waterways in Serbia. Plovput’s typical cross-section single-beam surveying equipment includes the following:

• Survey Echo Sounder with a depth measurement precision of 1 cm +/− 0.1% depth;
• DGPTS-RTK positioning with a measurement precision of +/−20 cm;
• Data-processing and post-processing applications.

VITUKI was a coordinator for hydrographic surveys on the Danube River until 2012, and the LDVWD is in charge of river management and navigation routes on the Danube reach downstream of rkm 1560 + 000. The hydrographic survey equipment of VITUKI included a single-beam ultrasonic echosounder with an accuracy of 2 cm and a total station on the shore. The accuracy of the total station was less than 2 mm at 200 m in the robotic survey mode. The equipment of the LDVWD consisted of a single-beam ultrasonic echosounder and an RTK GNSS receiver. The accuracy of both instruments was 2 cm.

The echosounder used to survey cross-sections in Serbia was mounted on a 5 m long boat, specifically equipped for this activity. All collected data (which contained the x, y, and z coordinates of the measured points) are in the State Geographical Coordinate System. Before the commencement of a survey, the coordinates of boundary points were entered into the software. Survey lines ran along predefined cross-sections. Data on location and depth (x, y, z) were transferred to the software for hydrographic surveys. The software synchronized the data constantly. Thus, the boat position was supplied in real time. The speed of sound in water was entered according to the information provided by an SVP (sound velocity profiler) device. A differential GPS station was mounted at the reference point with known geographic coordinates. A station was connected to the boat with the radio signal, sending information on differential correction and providing the required accuracy for the performed survey. Depth information was obtained by using the time necessary for ultrasound waves to travel from the echosounder to the river bottom and back. Data on depth and location were synchronized in real time, and information was stored in ASCII format in the form of x, y (position), and z (depth) coordinates. After completing the survey, quality control was performed, spikes were removed, and data were stored in the database with cleaned x, y, and z coordinates for each cross-section.

VITUKI surveyed cross-sections with a special motor boat. The total station on the occupied stance on the shore followed the active signal, coming from the echosounder beacon-receiver placed on the same rod with fixed height and without horizontal offset. The processing result was a 3D point (x, y, and z coordinate) list.

The LDVWD surveyed the bed with satellite positioning (RTK GNSS) and single-beam ultrasonic depth measurement equipment based on a specific cross-section plan prepared in advance. To convert the satellite positioning coordinates into a Cartesian spatial coordinate system, original positioning in the global system and the officially accepted Field Transformation Procedure were used.

3.3. Analysis of Morphological Changes

Since Hungary and Serbia use different surveying systems, it was necessary to make all data comparable. After the conformance of the data in the two countries, the cross-sections were visualized and processed using the AutoCAD software (Map 3D 2011) and Plovput’s in-house AutoCAD routines.

To establish a standard for comparisons along the stretch being studied, the next step involved defining a reference. To achieve this, Étiage Navigable (EN) was selected, which is not a reference level but a reference water surface elevation profile along the reach, belonging to a relatively low discharge and hydraulically steady state. EN was originally determined by the Danube Commission in 2012 for each gauging station. This requires that each cross-section has its own reference level. However, as EN is significantly below the bankfull stage, it was decided to raise it by 2 m. Once this was accomplished, the area of each cross-section below the newly established EN + 2 m was calculated.

The correlation between the decrease in water levels and the incision of the riverbed was analyzed using cross-sectional areas of gauging stations along the sand bed reach upstream of the Drava River confluence. The analysis was conducted for the selected years during the 60-year-long period (years that had hydrographic surveys performed in both countries). In total, eight different years were found, starting from 1975, and the first available hydrographic survey data in Serbia from 1964 were also included in the analysis of the GSs on the Serbian side. The rate of cross-sectional area change was analyzed first, followed by a correlation between the percentages of the cross-sectional increase and minimum water levels. The percentage of the cross-sectional increase was expressed with reference to the first year with available bathymetric data. All trend and correlation analyses were performed using linear regression with the LSQ method in MS Excel.

3.4. Indirect Method for Estimation of the Sediment Transport

Sediment transport was inferred indirectly from the 69 surveillance cross-sections from Baja to the confluence of the Drava River. The study is limited to this shortened sand bed reach, with no large tributaries with massive sediment inputs. Surveillance cross-sections are from the databases of the three previously mentioned institutions. The procedure for defining the cross-sectional area is described in the previous subsection. The estimation of the sediment transport should be based on data from the two years when hydrographic surveys were made in both countries. Since there was no overlap in the last 25 years, we decided to use data from surveys in two consecutive years in the two countries. These were datasets collected in 1996 and 2021 by Plovput and datasets collected in 1997 and 2022 by VITUKI and LDVWD, respectively. The datasets were used for the estimation of 25-year sediment transport changes along Serbian and Hungarian reaches. The areas of 69 cross-sections were used to calculate channel volumes between each pair of cross-sections along the study reach (sr):

$${\forall }_{sr,i}\left(t\right)=0.5\left(A_i\left(t\right)+A_{i+1}\left(t\right)\right)∆x,i=1,68$$

(1)

The long-term sediment transport between two cross-sections is thus

$$Q_{s,i}=\frac{{\forall }_{sr,i}\left(t_1\right)-{\forall }_{sr,i}\left(t_2\right)}{∆t}{\rho }_s\left(1-\lambda \right)=\frac{∆{\forall }_{sr,i}}{∆t}{\rho }_s\left(1-\lambda \right),i=1,68$$

(2)

where ${\rho }_s$ is sediment particle density and $\lambda$ is a porosity of the bed material.

4. Results

4.1. Results of the Homogeneity Tests for the Hydrological Data

Homogeneity tests for water level yearly characteristics were performed by Tamás et al. in 2021 [5], and the results are presented in Figure 3a–c only for the Baja GS, as all gauging stations along the investigated reach showed very similar behavior. Both minimum and mean water stages and levels were highly inhomogeneous, with large D-statistic values (>0.3) and very low probabilities (p) (less than 0.01; at some GSs, these were not even traceable; see

Table 1. Results of the homogeneity test of water level yearly characteristics [5].

Gauging Station	Minima			Means			Maxima
	D	p	Homog.	D	p	Homog.	D	p	Homog.
Dunaújváros	1.000	<0.0001	NO	0.943	<0.0001	NO	0.343	0.033	NO
Dunaföldvár	0.971	<0.0001	NO	0.857	<0.0001	NO	0.257	0.197	YES
Paks	0.714	<0.0001	NO	0.543	<0.0001	NO	0.159	0.837	YES
Baja	0.457	0.001	NO	0.400	0.007	NO	0.114	0.976	YES
Mohács	0.457	0.001	NO	0.429	0.003	NO	0.086	1.000	YES
Bezdan	0.371	0.016	NO	0.371	0.016	NO	0.114	0.976	YES
Apatin	0.429	0.003	NO	0.371	0.016	NO	0.114	0.976	YES
Bogojevo	0.429	0.003	NO	0.429	0.003	NO	0.171	0.683	YES
Bačka Palanka	0.257	0.255	YES	0.357	0.038	NO	0.121	0.976	YES
Novi Sad	0.229	0.320	YES	0.171	0.683	YES	0.114	0.976	YES

). For the sake of the completeness of the results and the explanation for why the evaluation of the sediment transport was limited to one part of the investigated reach, the D-statistics and p-values are taken from [5] and presented again in Table 1.

The homogeneity test was also performed for discharges, and the results are presented in

Table 2. Results of the homogeneity test of yearly discharge characteristics (taken from [5]).

Gauging Station	Minima			Means			Maxima
	D	p	Homogen	D	p	Homogen	D	p	Homogen
Baja	0.139	0.896	YES	0.143	0.879	YES	0.190	0.569	YES
Bezdan	0.143	0.867	YES	0.171	0.683	YES	0.114	0.976	YES

, as in [5], but only for the two gauging stations, as was explained in Section 3.1. They are graphically presented only for the Baja GS (Figure 3d–f), as they are very similar for all other GSs. Unlike the water level data, the yearly characteristic discharges at the investigated gauging stations are homogeneous for all three characteristics—minima, means, and maxima (all D-statistic values are less than 0.2, and values of probability (p) are much larger than 0.70). In the same study ([5]), it was shown that the annual discharge minima, means, and maxima showed no trend at the downstream end of the study reach over the last 70 years.

The results of the discharge normality test are visualized in Figure 4. Although the frequency of the measured data (green columns) shows slight asymmetry, the test undoubtedly indicates that all three yearly discharge characteristics follow normal distributions at all four gauging stations.

4.2. Results of the Hydrological Data Trend Analysis

The results of the trend analysis of the annual characteristic values of water levels are presented in Figure 5, based on data from [5]. The data show a linear decreasing trend at all gauging stations; the trendline slope is negative along the reach except in Novi Sad at the downstream end, where the free-flowing Danube enters the Iron Gate I reservoir. It is readily apparent that the rate of decrease in minimum water levels decreases in the downstream direction. The decreasing rate is more than three times greater on the gravel bed (approximately 3.9 cm/year) than on the sand bed sub-reach, where the average decreasing rate upstream of the confluence of the Drava River is 1.1 cm/year. Along the sand bed reach between Baja and Bogojevo (just upstream of the Drava River confluence), it almost linearly decreases from 1.3 cm/year to 0.9 cm/year. Based on these data, the total change in water level minima (ΔZ) for the analyzed 70-year period ranges from 1.45 m at the beginning of the sand bed reach (Paks) to approximately 0.60 m at Bogojevo (Figure 6). Between Mohács and Apatin, the total drop in water levels changes almost linearly from 0.76 m to 0.70 m. The corresponding average yearly water level decrease (δZ) is around 2 cm/year at Paks, 1.3 cm/year at Baja, and approximately 1 cm/year between Mohács and Apatin. Downstream of the Drava River confluence, it is less than 0.9 cm/year.

Variation of trends for the mean values on the gravel-bed sub-reach is almost the same as for minima, while that on the sand-bed reach is similar, except the rates are higher. Since means include the influence of water level maxima, any anomaly in the water level maxima reflects on them. It is visible at Mohács (Figure 5), where the rate of decrease is approximately 25% greater than expected (0.0146) had not there been this anomaly in the water level maxima.

The annual characteristic values of discharges on the reach between Baja and Apatin did not change significantly during the investigated period, as annual maxima show a slightly increasing trend, while mean values and minima fluctuate around almost constant values (Q_min around 1100 m³/s and Q_mean around 2300 m³/s) [5].

4.3. Correlation between the Minimum Water Level and Cross-Sectional Area Changes at GSs

The evolution of cross-sections at four GSs on the sand bed reach from Baja to Apatin is presented in Figure 7. Visual inspection of the cross-sections shows that they are more or less stable starting in 1985 (dark blue line). The most stable cross-section is in Baja (Figure 7a), where the amplitude of changes is the smallest. The cross-section in Apatin is the most unstable one, having the largest amplitudes in bed level changes (Figure 7d). The analysis of the cross-sectional area change compared with the first survey supports these visual observations (

Table 3. Change trends of the cross-sectional area with reference to the first survey, A_Y/A_ref − 100 = aYear − b; the reference years (ref) are 1975 for GSs in Hungary and 1964 for GSs in Serbia, except for Apatin, where 2004 is taken as the reference.

Gauging Station	River km	a	b
Baja	1478.70	0.08	152.36
Mohács	1446.90	0.24	489.29
Bezdan	1425.59	0.59	1172.1
Apatin	1401.90	1.92	−3890.6

). This revealed that the Apatin cross-section became stable after 2004 upon additional river training work. Thus, 2004 was used as the reference year for this GS. It is readily apparent that the rate of change (value of the parameter, a) increases in the downstream direction.

The results of the analysis of rates of change in the cross-sectional area between the two successive surveys (ΔA/Δt) at the four considered GSs are presented in Figure 8. They confirm the visual observations. The lowest rates are at Baja GS (Figure 8a). They are less than 20 m²/year or 0.75% of the cross-sectional area below EN + 2 m, except between 2004 and 2007 (23.5 m²/year) when there was a large flood wave in the spring of 2006 with very long overbank flow durations along the entire Danube River course. The most prominent morphological changes at the Mohács GS happened between the first three surveys—an aggradation of approximately 50 m²/year (1975–1985) and a degradation of around 65 m²/year (1985–1988) (Figure 8b). These changes correspond to 2% and 2.6% of the cross-sectional area each. After that, ΔA/Δt did not exceed 30 m²/year, except, again, between 2004 and 2007, when it was around 35 m²/year. Changes at Bezdan GS (Figure 8c) were less than 10 m²/year before 1997, i.e., less than 0.5% of the cross-sectional area. Between the following two surveys, the cross-section aggraded and degraded for equal amounts of approximately 45 m²/year or 2% of A. After the large flood wave of 2006, the changes in the cross-sectional area increased, and between the last two surveys, they exceeded 20 m²/year, or 0.8% of A. Finally, Figure 8d quantitatively confirms that the GS in Apatin is the most unstable in the reach. This section is one of 17 critical sectors for inland navigation in Serbia. The cross-sectional area aggraded at a rate of 137 m²/year between 1985 and 1988, which is 6.9% of A. The aggradation/degradation between other surveys was always greater than 50 m²/year, or more than 3% of A.

Finally, the correlation between the percentages of the observed cross-sectional area increase (morphological variable) and the minimum water level (hydrological variable) is presented in Figure 9. The cross-sectional area increase is calculated compared with the reference year, i.e., the first year with available bathymetric data, as indicated in the caption of Table 3. There is a linear correlation between the two variables, and the decrease in the low water level follows the incision of the riverbed. Values of regression coefficients are provided in

Table 4. Correlation coefficients in the linear regression equation, Z_min = a (A_Y/A_ref − 100) + b, where the reference years (ref) are 1975 for GSs in Hungary and 1964 for GSs in Serbia.

Gauging Station	River km	a	B	R2
Baja	1478.70	−0.124	8.21	0.64
Mohács	1446.90	−0.031	8.02	0.09
Bezdan	1425.59	−0.032	8.10	0.62
Apatin	1401.90	0.004	7.92	0.02

. Figure 9 confirms what was already presented in Figure 5 and Figure 6, i.e., that the highest rate of water level decrease is in Baja (rkm 1478.70) and that Z_min decreases at a rate of 0.12 m² with a percentage increase in the cross-sectional area, with a coefficient of determination of R² > 0.6 (Table 4). The rate of Z_min decrease is almost the same as in Mohács (rkm 1446.8) and Bezdan (rkm 1425.59). It is approximately 0.03. However, the coefficient of determination in Bezdan (R² = 0.62) is significantly greater than that in Mohács (R² ≈ 0.09). The trend in Apatin (rkm 1401.90) decreases for the entire period (1964–2016), i.e., the slope of the trendline is positive but very small (a = 0.004). Moreover, the coefficient of determination is very low (R² ≈ 0.02). These values altogether indicate that there is no correlation between the morphological and hydrological variables at this particular GS.

4.4. Estimation of the Sediment Transport

Variations in the cross-sectional area (A) along the study reach in the two consecutive years 1996/1997 and 2021/2022 (as explained in Section 3.4) are presented in Figure 10. Generally, the cross-sectional area (A) increased in 52 cross-sections (75%), while it decreased in the remaining 17 cross-sections (25%). The longest stretch with the observed decrease in A is between Baja and Mohács (11 cross-sections). The percentage decrease reduces in the flow direction from 21% just downstream of Baja to 1% just upstream of Mohács. The remaining six cross-sections are on the Serbian reach. The percentage reduction is generally less than 3%. However, there are two locations with 12 and 17.5% reductions (just upstream of Bezdan and some 20 km upstream of the Drava River confluence, respectively). The stretch with the most pronounced increase in A is between Bezdan and Apatin, with a maximum enlargement of A of approximately 30%.

These longitudinal profiles of A were used to estimate the 25-year long-term sediment transport, presented in Figure 11. The largest multiannual sediment transport rate (Q_s) (around 65∙10³ t/year) is on the Hungarian side between Baja and Mohács, where the aggradation of the riverbed is evident (Figure 10). The most intensive degradation of the riverbed is downstream of Bezdan and resulted from a multiannual sediment transport rate of approximately 45∙10³ t/year. Generally, downstream of Mohács, the sediment transport rate is negative, meaning there is a continuous uptake of sediments from the riverbed, i.e., the incision of the riverbed. The uptake rates are mostly greater than 10∙10³ t/year, and between some surveillance cross-sections, they reach approximately 25∙10³ t/year.

5. Discussion

The minima and means of water levels are inhomogeneous along the entire reach except at Bačka Palanka and Novi Sad, where the Danube approaches and enters the Iron Gate Reservoir and where the effects of the backwaters are or could be felt (Table 1). The non-homogeneity of water levels results from their strong relationship with riverbed changes, which are the consequence of intensive river training work along this stretch of the river. On the other hand, the yearly maxima are homogeneous, as the shape of the main-channel riverbed at high flows merely influences the water level, and it is supposed that the riverbed incision mainly affects the riverbed of the main channel and not the floodplains. Unlike water level data, the yearly characteristic discharges at the investigated gauging stations are homogeneous for all three characteristics—minima, means, and maxima. The homogeneity of the discharge data (Figure 3d–f) can be explained by the fact that the shape of the riverbed does not affect water quantity, which arrives from upstream, and the fact that there are no large tributaries downstream of Szob (rkm 1700).

The trend analysis of a 70-year long time series for the annual characteristic water level values of the entire free-flowing Middle Danube reach [5] confirmed the decreasing trends in all three characteristics (minima, means, and maxima) found in [4,17] along the Hungarian stretch until 1992 and 2005, respectively. Inspecting the variation in the rates of change of Z_min and Z_av along the studied reach reveals that rates along the gravel bed reach are up to three times larger than those along the sand bed reach. Between Mohács and Bogojevo, the ΔZ/Δt is almost the same for Z_min, with an average value of 1.1 cm/year. Downstream of the Drava River confluence (in Bačka Palanka) and at the entrance to the Iron Gate I reservoir (in Novi Sad), the decreasing rate rapidly reduces because of the deposition of sediments brought by the Drava River and the backwater effects of the Iron Gate I Dam. Moreover, under the influence of Iron Gate I backwaters, the incision terminates and changes into aggradation, an usual physical process where a river enters a reservoir, lake, or sea. Thus, this is a positive answer to the first research question.

The second research question about the correlation between the observed minimum water levels and the recorded incision of the riverbed resulted from the fact that the water level decrease is usually a consequence of a channel incision. The riverbed incision is a continuous process in rivers with numerous cut-offs. Rectification of the river course by cutting-off meanders and disconnection of narrower side channels along anabranching river reaches was a river engineering technique applied in the 19th century aiming at facilitating navigation, flood control and agricultural production along European rivers that played a role in demographic, economic, cultural and urban development, such as the Rhine and Danube Rivers. Dredging and groyne construction might amplify the incision process, both of which happened in the second half of the 20th century along the Middle Danube. The comparison between the longitudinal profiles of the sand bed free-flowing Middle Danube River reach in 1996/1997 and 2021/2022 (Figure 12) shows that the riverbed incision continued after 2013 when a report on long-term morphological development was published [21] (this was approximately 110 years after the start of cutting-off work on the Middle Danube). Long-term analysis of longitudinal riverbed profiles of the Upper Rhine River has shown, at one location, a drop in thalweg as large as 7 m 120 years after Johann Tulla’s work [22]. Thus, a recorded maximum thalweg drop of approximately 6 m in a bend upstream of the Drava River confluence (Figure 12) agrees with observations in other rectified rivers.

In addition to actual longitudinal riverbed profiles, large-scale bed slopes at the beginning and end of the considered time interval are drawn in Figure 12. These readily demonstrate that the large-scale bed slope of the study reach is decreasing, i.e., that the river channel is still responding to changes from 150 and 60 years ago. The average bed slope reduced from 0.094‰ to 0.083‰.

To the best of the authors’ knowledge, there is no information on the correlation between the minimum water level and the incision of the channel at GSs on the Danube River. In this study, a linear relationship between the Z_min and the corresponding percentage of the cross-sectional area increase due to a channel incision after the reference year of observation was found at three out of four considered gauging stations. The exception is the Apatin GS, where the channel dynamically changes with the growth and washing of a large sand bar. This location is one of the critical sectors for navigation along the Serbian/Croatian border reach. The rates of Z_min lowering decrease in the downstream direction, which is consistent with the results of the trend analysis in Figure 5 and Figure 6 and in [6].

Since the Danube River channel on the investigated reach was predominantly of the meandering type in the 19th century, it was rectified with numerous cut-offs (Figure 4.1.12 in [21] (p. 42)). Only on the Hungarian side were there more than 20 cut-offs between Dunaföldvár and the Hungarian–Serbian border [23]. Thus, an alternating degradation and aggradation pattern could be expected. This pattern is further amplified by the construction of river-training structures starting in 1876 and dredging in the second half of the 20th century. The answer to the third research question, i.e., what are the long-term sediment transport characteristics along the investigated reach, is sought through the indirect estimation of the sediment transport rate between Baja and the Drava River confluence, as there are no large tributaries along the reach that could affect riverbed morphology. Both figures—Figure 10 (with the raw data) and Figure 11 (with the results of the sediment transport estimation)—show an expected pattern. The stretch between Baja and a location some distance upstream of Mohács aggrades (the main-channel cross-sectional area reduces—Figure 10 and positive values of Q_s in Figure 11, which are more than four times greater than the negative ones along the reach exposed to degradation). It is expected since it is located downstream of a stretch with numerous cut-offs and intensive dredging in the second half of the 20th century (the area close to Dunaföldvár). The short river reach around Baja can be considered stable (Figure 7a and Figure 8a), possibly owing to the very wide floodplain (~6 km), where large discharges with large velocities (high energy level) can spread. This area is known as Gemenc from ~1503 to ~1469 rkm. Downstream of Gemenc, the river channel is narrow until it reaches Mohács. Close to Mohács, a reach with numerous cut-offs and river-training structures built to aid navigation, it starts again (see the map on the right in Figure 4.1.12 in [21] (p. 42)). As such, it continuously degrades with a mean 25-year long-term sediment transport of approximately 10 t/year. The most intensive washing of the riverbed material during the analyzed period (around 50 t/year) is throughout the length of 5 km downstream of Bezdan (even the translation of the name of this small town—bottomless—suggests that there is a permanent deepening of the riverbed at this location). This is logical, as Bezdan is at the upstream end of the 25 km rectified river channel. The sediment transport rate ceases toward the Drava River confluence. Another quantitative drop in the washing rate is close to Apatin, where the sediment transport drops from 13.5 t/year to 3 t/year on average.

During the investigated 250-year-long time interval, several large floods with peaks exceeding 7000 m³/s happened (in the summer of 2002, the spring of 2006, the winter of 2010, and the biggest one in the summer of 2013). Among them, the one with undoubtedly the most prominent effect on the morphological changes was that of the spring of 2006, when there were 48 days with discharge larger than the bankfull discharge of 4500 m³/s and 5 days with a discharge larger than 7500 m³/s (see Figure 8, block VI, for the time interval between 2004 and 2007, when there was a significant deepening of the cross-sections along the entire reach). Since large amounts of sediment are carried during extreme flood events, it would be interesting to see whether and how it affected the sediment transport rate along the study reach and, consequently, sedimentation and erosion patterns. As was already mentioned, there were only two almost-coinciding years with bathymetric surveys in both countries at the disposal for the analysis—1996/1997 and 2021/2022. Thus, it would even be useful to analyze asynchronous data to infer to what extent floods alleviate or deteriorate the incision of the riverbed in the main channel and how much time is needed after the flood event to reach dynamic equilibrium. From the standpoint of this analysis, it would be advisable if bathymetric surveys were performed along the entire reach after each large flood wave.

6. Conclusions

A comprehensive hydrological–morphological study of the 300 km long free-flowing Middle Danube River reach between rkm 1581 (Dunaújváros), Hungary, and rkm 1255 (Novi Sad), Serbia, led to the following conclusions:

• The data series of minimum and average water levels are inhomogeneous and show a continuous decreasing trend. They are both the consequence of the permanent riverbed changes initiated by the rectification of the meandering river channel in the 19th century, and they amplified with the consequent construction of groynes and longitudinal stone toes since 1876 and intensive dredging after the 1950s.
• The rate and type of cross-sectional area change (an increase or a decrease) depend on the location of the cross-section along the river course. The slowest decrease is downstream of reaches with cut-offs and intensive dredging. More intensive changes in the cross-sectional area are at the upstream ends of rectified main-channel reaches.
• The correlation between low water levels and the percentage increase in the main channel cross-sectional area at official gauging stations in the last 50 years is linear. Thus, the hypothesis about the continuous lowering of water levels caused by the incision of the main channel is correct. The large-scale longitudinal riverbed slope decreased in the last 25 years by 0.01‰. A maximum water level decrease of approximately 1 m is on the upstream end of the study reach, and it reduces when approaching the Drava River confluence.
• In the absence of systematic sediment transport monitoring, regularly surveyed cross-sections of the navigable river can be used to estimate the long-term averaged sediment transport, i.e., to infer the average rates of aggradation and degradation of the channel along the reach.
• Along the sand bed reach of the free-flowing Middle Danube with no large tributaries, the maximum value of the positive sediment transport in the last 25 years is approximately 30% greater than the maximum absolute value of the negative transport.
• The maximum absolute value of the negative transport is near Bezdan, at the upstream end of another rectified river sub-reach. The sediment transport ceases toward the confluence of the Drava River.
• In 86% of the study reach, the long-term sediment transport for the studied period is negative. Thus, sediment balance changes also prove that decreasing water levels are a good indicator of an ongoing severe riverbed erosion process.