Remote Sensing of Sediment Discharge in Rivers Using Sentinel-2 Images and Machine-Learning Algorithms

The spatio-temporal dynamism of sediment discharge (Qs) in rivers is influenced by various natural and anthropogenic factors. Unfortunately, most rivers are only monitored at a limited number of stations or not gauged at all. Therefore, this study aims to provide a remote-sensing-based alternative for Qs monitoring. The at-a-station hydraulic geometry (AHG) power–law method was compared to the at-many-stations hydraulic geometry (AMHG) method; in addition, a novel AHG machinelearning (ML) method was introduced to estimate water discharge at three gauging stations in the Tisza (Szeged and Algyő) and Maros (Makó) Rivers in Hungary. The surface reflectance of Sentinel-2 images was correlated to in situ suspended sediment concentration (SSC) by support vector machine (SVM), random forest (RF), artificial neural network (ANN), and combined algorithms. The best performing water discharge and SSC models were employed to estimate the Qs. Our novel AHG ML method gave the best estimations of water discharge (Szeged: R2 = 0.87; Algyő: R2 = 0.75; Makó: R2 = 0.61). Furthermore, the RF (R2 = 0.9) and combined models (R2 = 0.82) showed the best SSC estimations for the Maros and Tisza Rivers. The highest Qs were detected during floods; however, there is usually a clockwise hysteresis between the SSC and water discharge, especially in the Tisza River.


Introduction
Sediment is a substance that has been eroded from rocks by physical and chemical processes and is moved by gravity, water, wind, or by their combination. The term "fluvial sediment" refers to the sediment transported by rivers [1]. Fluvial sediment has many benefits, as it controls the global carbon cycle, morphological evolution of river channels, and delta development. Furthermore, it improves the agricultural productivity of the floodplain soils due to the deposition of nutrient-rich material [2]. In natural systems, climate and tectonics are the main driving factors of continental erosion [3] and thus of fluvial sediment availability [4]. However, the various anthropogenic activities (e.g., land use/land cover changes, construction of dams, mining activities, runoff, and erosion control) and climate change alter the quantity and timing of the fluvial sediment formation and transport. For instance, Martinez et al. [5] studied the sediment discharge of the Amazon River between 1995 and 2007, and the findings revealed that the suspended sediment concentrations rose by 20% due to changes in land use/land cover and global rainfall patterns. On the other hand, others [6,7] reported a considerable decrease in the sediment discharge of the Yellow and Mississippi Rivers due to the construction of sediment trapping systems and soil conservation activities.
Quantifying the amount of transported sediment in rivers plays an important role in various engineering works. For example, reservoirs should be designed to be large enough to accommodate the expected deposited sediment throughout their lifetimes, or the drop of and discharge, and the available historical data of discharge), with RMSE ranges between 4.1% and 13.8%. The river width is considered the most obvious hydraulic parameter that can be observed from space. Gleason and Smith [31] estimated river discharge at mass conserved reaches along many rivers by width measurements obtained from satellite images and at-a-station hydraulic geometry (AHG; Leopold and Maddock [12]). However, they discovered a log-linear relationship between the AHG coefficients and exponents of the cross-sections along mass conserved reaches of rivers, which is called at-many-stations hydraulic geometry (AMHG). In addition, if these coefficients and exponents are optimized to minimize the discharge difference between every pair of cross-sections, better discharge estimations will be obtained. Therefore, they estimated the reach discharge based on optimized AHG coefficients and exponents of several cross-sections. Durga Rao et al. [35] estimated the discharge of four Indian rivers at different discharge conditions using a time series of width data obtained by Landsat and IRS satellite images and the AMHG approach. The AMHG coefficients and exponents were optimized by a genetic algorithm, and the results showed good accuracies, especially during low and medium stages. Comparably, Mengen et al. [32] followed a similar approach to estimate the discharge of the Mekong River, utilizing Sentinel-1 data (active sensor) to obtain river width instead of passive sensor data.
In order to calculate or measure the suspended sediment discharge of a river, the sediment characteristics and its transport rate should also be known. The suspended sediment usually has a size range between 0.45 µm and 2000 µm [40]. In general, the bedload forms 5-20% of the total sediment transport of rivers [40], whereas the rest forms the suspended load and wash load. However, this percentage range is quite rough, as, in lowland rivers with abundant suspended sediment supply, the proportion of bedload could be less than 1% [41]. The suspended sediment concentration in rivers can be estimated by collecting water samples manually or automatically (e.g., passive samplers, automated pumping, or time-integrating samplers). The manual method can be implemented under the depth integrated approach, by collecting samples from the whole depth [42,43], or under the point integrated approach, when a sample is collected from a single point at the water surface [44]. Despite the great accuracy of the in situ methods for measuring suspended sediment concentration, they are very time consuming and labor intensive; moreover, their results are just limited to the given sampling station. As a result, many studies investigated the suspended sediment concentrations in rivers using remote-sensing images, exploiting the ability of suspended sediment to increase reflectance values at visible and near-infrared spectra [45][46][47][48][49][50]. The correlation between suspended sediment concentration and reflectance of satellite images was performed by semi-empirical approach [46,51,52], bio-optical modeling approach [53,54], or machine-learning algorithms [55][56][57].
Few studies attempted to retrieve sediment discharge using the remote-sensing method [2,3,5,[58][59][60]. Some of them estimated sediment discharge by combining in situ measurement of discharge and remotely sensed suspended sediment concentration [5,[58][59][60]; others used satellite gravimetry data to estimate sediment discharge of large rivers [3]; and fewer obtained both water discharge and suspended sediment concentration from satellite images to compute sediment discharge [2].
In Hungary, the medium-sized Tisza and Maros Rivers transport a large amount of suspended sediment (Tisza: 18.7 million t/y; Maros: 8.3 million t/y). The dominance of suspended sediment in their sediment budget is well reflected by the low proportion of their bedload transport (Tisza: 0.04%; Maros: 0.3%). Unfortunately, only 8 sediment gauging stations operate on the 962 km long Tisza and only one on the 750 km long Maros (Figure 1). In addition, the measurements are performed just once a month, regardless of the water stage. However, more frequent spatial and temporal data are needed to support environmental studies of the river system. Therefore, this study aims to calibrate and validate suspended sediment and discharge models for the Tisza and Maros Rivers by in situ and Sentinel-2 data. The derived models could be employed to estimate sediment discharge at ungauged periods. The detailed goals of this study are to: (1) derive the river discharge models based on AHG and AMHG methods, in addition to assessing the applicability and stability of discharge models Hydrology 2022, 9, 88 4 of 30 derived by combing AHG with machine-learning algorithms, which was not assessed in the literature; (2) calibrate the suspended sediment models based on machine-learning algorithms, such as artificial neural network (ANN), support vector machine (SVM), and random forest (RF), comparing their results and performance; (3) apply the best derived suspended sediment and discharge models to historical images of Sentinel-2 (ungauged periods) to produce a longer time series of sediment discharge of both rivers within the study area.
Hydrology 2022, 9, x FOR PEER REVIEW 4 of 32 Hydrology 2022, 9, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/hydrology the water stage. However, more frequent spatial and temporal data are needed to support environmental studies of the river system. Therefore, this study aims to calibrate and validate suspended sediment and discharge models for the Tisza and Maros Rivers by in situ and Sentinel-2 data. The derived models could be employed to estimate sediment discharge at ungauged periods. The detailed goals of this study are to: (1) derive the river discharge models based on AHG and AMHG methods, in addition to assessing the applicability and stability of discharge models derived by combing AHG with machinelearning algorithms, which was not assessed in the literature; (2) calibrate the suspended sediment models based on machine-learning algorithms, such as artificial neural network (ANN), support vector machine (SVM), and random forest (RF), comparing their results and performance; (3) apply the best derived suspended sediment and discharge models to historical images of Sentinel-2 (ungauged periods) to produce a longer time series of sediment discharge of both rivers within the study area.  The study focuses on the lower reaches of both rivers at their confluence area: the reach of the Tisza between Algyő and Szeged (201.5-172.8 f km) and the lower section of the Maros between Makó and the confluence (24.5-0 f km) (b). At the fluviometers (red circle), water stage and discharge (Q) data are collected daily, while suspended sediment concentrations (SSC) are measured monthly (yellow triangles). The water discharge and SSC measurement points are at the same place at the Makó fluviometer, while at Szeged, they are spaced ca. 1 km from each other.

Study Area
The presented research studied the suspended sediment transport of the Tisza and Maros Rivers in the eastern part of the Carpathian Basin, central Europe. The Tisza River is the longest tributary (962 km) of the Danube [16], and its catchment area (157,200 km 2 ) covers 19.5% of the Danube's catchment area [61]. The Maros/Mures River is the greatest tributary of the Tisza (length: 750 km; catchment area: 30,000 km 2 ) and one of the most dynamic rivers of the Carpathian Basin [62]. The study was performed at their confluence area ( Figure 1). The Maros joins the Tisza at 177 fluvial km (f km) of the Tisza, with a confluence angle of 23 • [47]. The Lower Tisza is wider (164 m on average) and considerably deeper (14.1 m) than the Lower Maros (width: 126 m; depth: 4.26 m; [18,63]).
The flood of the Maros usually precedes the flood of the Tisza by a few days; however, they could also occur simultaneously during long-lasting floods [64]. As a result of their low slope conditions, they could be impounded by the Danube or each other [47,62]. Usually both rivers have a flood period in early spring (March-April) caused by snowmelt and in early summer (June-July) due to summer rainfall. Usually, the early spring flood conveys more discharge (15% of the annual transported water [65] and sediment) than the early summer flood [16,62]. Sometimes the summer flood superimposes the spring flood (e.g., in 2006) because of the low slope (2.9 cm/km; [16,66]) on the long lowland section of the Tisza River, which results in high flood waves [62]. The floods of the Tisza usually cover the floodplain for 2-3 months; however, they extend for only 1-2 weeks in the Maros River due to its greater slope (Lower Maros: 28 cm/km; [16,66] [16]). There is a huge difference between the lowest and highest stages, which is 13.55 m on the Tisza and 7.28 m for the Maros [16].
Both rivers transport a substantial amount of suspended sediment. Based on in situ measurements , the Lower Tisza's mean annual suspended sediment load is 18.7 million t/y, and its mean concentration is 640 g/m 3 . However, in the confluence area with the Maros, the sediment conditions change considerably [67], as the Maros transports 8.3 million t/y suspended sediment (mean concentration: 1600 g/m 3 ) into the Tisza [16,47,68].
There are considerable spatio-temporal variations in the amount of transported suspended sediment along the Tisza River ( Figure 2). The main reasons for the spatial variations between two subsequent gauging stations are the aggradation of suspended sediment in the channel during low stages, which can be mobilized by small flood waves, and sedimentation on the floodplain [69]. The suspended sediment budget is probably also influenced by deposition in reservoirs and bank erosion, though they were not evaluated. According to the measurements on the Middle Tisza, at Szolnok and Kisköre [69,70], the suspended sediment concentration rapidly grows during the rising limb of the flood waves, and it reaches its maximum before the peak of the flood ( Figure 2). However, the suspended sediment transport conditions of even very similar floods are different [70], which makes the predictions difficult. June 1999 and 2001) (Data source: [69,70]).

In Situ Hydrological Data
The hydrological data were collected at three gauging stations. Daily water stage and discharge data were available from Szeged (173.6 f km) on the Tisza River and from Makó (24.5 f km) on the Maros River for the period between January 2015 and May 2021. Upstream of the confluence of the two rivers at Algyő (192 f km), only daily water stage data

In Situ Hydrological Data
The hydrological data were collected at three gauging stations. Daily water stage and discharge data were available from Szeged (173.6 f km) on the Tisza River and from Makó (24.5 f km) on the Maros River for the period between January 2015 and May 2021. Upstream of the confluence of the two rivers at Algyő (192 f km), only daily water stage data (January 2015-January 2020) were available. Therefore, here, the discharge was estimated based on a rating curve. On the Lower Tisza (ca. 1 km downstream of the fluviometer at Szeged (Figure 1; Q and SSC locations)) and the Maros (at Makó), water discharge and suspended sediment concentration (SSC) data were also collected monthly between February 2015 and November 2020. The water discharge was measured by acoustic doppler current profiler (ADCP) and usually used to calibrate the rating curves at the fluviometers. Concurrently, the SSC data were obtained by depth integrated suspended sediment samples collected using a pump. The locations of the SSC sampling points were slightly different from one campaign to another as a result of varying stages. Therefore, the campaigns were grouped into high and low stages, and the average distances between the sampling points (or a sampling point and the riverbank) for every stage condition were considered ( Figure 3). All data were obtained from the Lower Tisza Hydrological Directorate (ATIVIZIG).
Hydrology 2022, 9, x FOR PEER REVIEW 7 of 32 Figure 3. Locations of the 100 transects, which were used to develop the AMHG models at the vicinity of Algyő (a), Szeged (b), and Makó stations (c). The suspended sediment concentration was measured at Szeged (d) and Makó (e). The five yellow transects at three stations were averaged to develop the AHG power-law and AHG ML models. The red and green horizontal lines on the crosssectional profiles indicate the high and low stages; the red and green vertical lines (V1-V5) refer to the suspended sediment sampling locations during high and low stages.

Remote-Sensing Data
The Sentinel-2A and 2B satellite images (provided by the Copernicus open-access hub; [71]) acquired between February 2015 and November 2020 covering the study area were applied in the study. Cloudy images were excluded; therefore, altogether, 122 cloudfree images were utilized. As the in situ campaigns were performed monthly, without considering the satellite overpass dates, a time frame of ± three days between the in situ campaigns and the Sentinel-2 acquisition dates was applied, which resulted in 29 images Locations of the 100 transects, which were used to develop the AMHG models at the vicinity of Algyő (a), Szeged (b), and Makó stations (c). The suspended sediment concentration was measured at Szeged (d) and Makó (e). The five yellow transects at three stations were averaged to develop the AHG power-law and AHG ML models. The red and green horizontal lines on the cross-sectional profiles indicate the high and low stages; the red and green vertical lines (V1-V5) refer to the suspended sediment sampling locations during high and low stages.

Remote-Sensing Data
The Sentinel-2A and 2B satellite images (provided by the Copernicus open-access hub; [71]) acquired between February 2015 and November 2020 covering the study area were applied in the study. Cloudy images were excluded; therefore, altogether, 122 cloudfree images were utilized. As the in situ campaigns were performed monthly, without considering the satellite overpass dates, a time frame of ±three days between the in situ campaigns and the Sentinel-2 acquisition dates was applied, which resulted in 29 images for the Tisza and 27 images for the Maros.

Pre-Processing of the Satellite Images
All obtained satellite images were available in level 2A (i.e., atmospherically corrected), except those acquired in 2015 and 2016. Thus, the Sen2Core 2.6 plugin, integrated into SNAP 8.0 software [74], converted these images from level 1C to level 2A.
Afterward, all images were resampled into 10 m spatial resolution and clipped within the study area ( Figure 1) to reduce the processing time. The water pixels of the rivers were extracted from the images by applying the Normalized Difference Water Index (NDWI; [75]); furthermore, the OTSU thresholding method introduced by Nobuyuki Otsu [76] was applied to obtain the threshold of the NDWI automatically based on the histogram. The scikit-image library in Python was used to apply the OTSU method.

Water Discharge Estimation
As we aim to estimate sediment discharge from satellite images, various models were built to estimate the water discharge and the suspended sediment concentration based on in situ measurements and Sentinel-2 images. The entire methodology applied in the study to estimate sediment discharge was summarized ( Figure 4).

At-a-Station Hydraulic Geometry (AHG), Power-Law, and Machine-Learning Models
At-a-station hydraulic geometry (AHG) is applied at any given cross-section of a natural stream, and it refers to the power-law relationship (Equations (1)-(3)) between water discharge and some hydraulic parameters, such as water depth, width, and velocity (Leopold and Maddock, 1953).
where w is water surface width, d is mean water depth, v is mean water velocity, Q is water discharge, a, c, and k are numerical coefficients, and b, f, and m are exponents of the equations. The main differences between the derived equations for any cross-section are the values of the coefficients and exponents. It can be noticed that the multiplication of these coefficients (a.c.k) and the sum of the exponents (b + f + m) are equal to 1.
This study focuses on the relationship between water width and discharge, as the satellites, with their continuous observation of the ground, can be employed to produce a time series of water width for any river. The selected satellite images for the study area (i.e., Tisza River: 29 images; Maros River: 27 images) were utilized to produce a time series of the river width at the three gauging stations (Figure 1). Five widths spaced at ±10-20 m upstream and downstream of the original cross-section were measured and averaged to reduce the error, which could arise from a single width measurement (Figure 3a-c). The RivWidth V04 software [Pavelsky and Smith [77] (http://uncglobalhydrology.org/ rivwidth/; accessed on 2 January 2022)] was used to measure the river width of the Tisza and Maros Rivers (at every 10 m). Afterward, the calculated width and the corresponding discharge at the three stations were fitted by a power-law function in Python using the curve fit function in the SciPy library.

Pre-Processing of the Satellite Images
All obtained satellite images were available in level 2A (i.e., atmospherically corrected), except those acquired in 2015 and 2016. Thus, the Sen2Core 2.6 plugin, integrated into SNAP 8.0 software [74], converted these images from level 1C to level 2A.
Afterward, all images were resampled into 10 m spatial resolution and clipped within the study area ( Figure 1) to reduce the processing time. The water pixels of the rivers were extracted from the images by applying the Normalized Difference Water Index (NDWI; [75]); furthermore, the OTSU thresholding method introduced by Nobuyuki Otsu [76] was applied to obtain the threshold of the NDWI automatically based on the histogram. The scikit-image library in Python was used to apply the OTSU method.

Water Discharge Estimation
As we aim to estimate sediment discharge from satellite images, various models were built to estimate the water discharge and the suspended sediment concentration based on in situ measurements and Sentinel-2 images. The entire methodology applied in the study to estimate sediment discharge was summarized ( Figure 4). . Flowchart summarizing the entire methodology followed in the study to produce water discharge, suspended sediment concentration, and sediment discharge data. Figure 4. Flowchart summarizing the entire methodology followed in the study to produce water discharge, suspended sediment concentration, and sediment discharge data.
The AHG theory was combined with the machine-learning algorithms; thus, the artificial neural network (ANN), support vector machine (SVM), random forest (RF) algorithms were employed to produce predictive models of water discharge based on the time series of width. The train-test-split function in the Scikit-learn library was used to divide the data into 80% for calibration and 20% for validation.
The artificial neural network (ANN) consists of three layers: input layer, hidden layer, and output layer [78]. Each layer contains several cells (neurons), which receive a signal (impulse) from the preceding neuron to deliver their result after performing some transformations into the subsequent neuron. These neurons are connected by weights, which are randomly selected at the beginning of the training process, and then it updates at every iteration until reaching the optimum weights [79]. The neural network in the scikit-learn library in Python was used to build the discharge model. One input layer containing the width values, two hidden layers with 20 neurons each, and one output layer producing the estimated discharge were applied. The rectified linear unit (relu) activation function, which controls the weighted sum of the inputs to feed it to the subsequent neuron, was applied for the hidden layers. The "relu" function returns the same input value in all cases, except for the values ≤0, as it returns zero [f(x) = max (0, z)], thus it avoids the saturation problem in other activation functions (e.g., sigmoid and tanh). The limitedmemory Broyden-Fletcher-Goldfarb-Shanno (lbfgs) solver, which optimizes the weights of every neuron to minimize the cost function, was also applied to the hidden layers. The "lbfgs" updates weights by calculating the inverse of the Hessian matrix. The maximum number of iterations was set to 200.
The support vector machine (SVM) finds the best fit line (hyperplane), which passes through the maximum number of the fitted points [80]. However, the main difference between the SVR and the other regression models is that these models try to find the best fit line, which minimizes the square error between the fitted line and points. Meanwhile, the SVR gives some tolerance to the solution, as it tries to find the best fit line within a threshold value (ε), which increases the chance of obtaining better solutions [81]. The kernel is responsible for finding the hyperplane, especially in problems with higher dimensional space. This study utilized the SVR in the scikit-learn library in Python to produce the discharge model. The polynomial kernel (second degree) was applied; the ε was set to 0.01 and the regularization parameter (c) to 2.
The random forest (RF) was introduced by Breiman [82]. The core of the RF is the decision tree (DT) algorithm, as the RF consists of many DTs working in parallel. The mechanism of the RF is that the input features of data are divided into several DTs. Then, every tree gives its prediction based on the included features, and the final prediction of the RF is the average of all DTs results. The RF avoids overfitting, which could result in a single DT; moreover, it can enhance the accuracy of the model. Many hyperparameters can be tuned to achieve the best performance of the RF model, such as the number of estimators, the maximum feature, minimum sample split, and minimum sample leaf. Our study used the random forest regressor module in the scikit-learn library in Python to build the discharge model. The number of estimators was set to 300, the maximum feature to auto mode, the minimum sample split to 2, and the minimum sample leaf to 4. The criterion used to control the quality of the split was the mean square error (mse).
The machine-learning (ML) models were integrated by averaging their results to produce a combined model. The combination process was implemented in Python. All the settings and parameters implemented in these ML algorithms were applied after fine tuning the hyperparameters until reaching the best available results.

At-Many-Stations Hydraulic Geometry (AMHG) Models
Based on Leopold and Maddock [12], the water width at any river cross-section is correlated with the discharge by a power-law relationship with a specific coefficient a and exponent b (Equation (1)). However, Gleason and Smith [31] discovered that if more than one cross-section is considered, their coefficients a and exponents b are not randomly related to each other; instead, a log-linear relationship could represent their relationship. It was coined as "at-many-stations hydraulic geometry" (AMHG). Suppose coefficient a and exponent b are optimized to minimize the discharge difference between the crosssections. The reach discharge can be obtained with better accuracy than applying the classical hydraulic geometry AHG, as the estimation is based on several cross-sections [2].
In this study, 1 km long reaches at the three gauging stations were considered, where the stations are located at the middle of the reach ( Figure 3). Then, a time series of water width for 100 transects at each gauging station was produced by the RivWidth software (Figure 3a-c). For every transect, a power-law relationship was fitted between the water width and its corresponding discharge at a particular date (which is the same discharge for all cross-sections). The fitting process was performed in Python using the curve fit function in the SciPy library, which produced 100 a and b values for every gauging station. A log-linear relationship was established between the log (a) and b; then, their limits were identified.
The optimization process of a and b values was performed by the genetic algorithm (GA) [83]. The GA depends on biologically inspired operators, including selection, crossover, and mutation, to develop high-quality solutions to optimization problems. The GA has an objective function, which maximizes or minimizes its value by selecting individuals (solutions) from a population. Every solution is called a chromosome, which consists of a number of genes [84]. Initially, random chromosomes are drawn from the population; then, they are passed through selection, cross-over, and mutation for a user-defined number of generations to finally produce the best solution to the problem.
For every station, the GA was applied between every pair of cross-sections (every cross-section is paired with the 100 cross-sections), with an objective function of minimizing the discharge difference between every pair. The population bounds of a and b values were defined based on the log-linear relationship produced previously, and the number of chromosomes was set to 10. Every chromosome contains four genes, which are the a and b values for every cross-section. The GA starts the optimization process by selecting ten chromosomes randomly from the population pool. Then, these chromosomes are combined with the water width to check the validity of the selected chromosome, as it should produce a reasonable discharge value between the minimum and maximum discharge of the river. The chromosome that produces a discharge outside the discharge limit is replaced until reaching ten valid chromosomes. The river discharge bounds can be identified based on the historical discharge values for the river or estimated empirically by Equations (4) and (5) [31]. The minimum and maximum discharge for the Tisza River were set to 50 m 3 /s and 4000 m 3 /s, while those for the Maros were 10 m 3 /s and 850 m 3 /s, respectively.
where w min and w max are the minimum and maximum observed width of the river during the studied period. The ten chromosomes are then ranked according to their fitness (discharge conservation) in a process called selection. Afterward, each chromosome has a chance to crossover with another chromosome or to mutate. The crossover means that certain genes in the chromosome can exchange their values with genes in another chromosome, while the mutation means that a few genes can be randomly altered to make sure that additional solution is added, which can be helpful in case of bad solutions in the existing chromosomes. The crossover rate is set to 0.8 and the mutation rate to 0.1 [31]. The altered chromosomes are reranked again based on their new fitness, and the least fit five chromosomes are excluded and replaced with another five from the population space. The first generation is passed when the chromosomes pass these previous steps (selection, crossover, and mutation). In our study, we applied 20 generations for every cross-section pair.
The GA was applied to all possible pairs of cross-sections within the 1 km reach, which resulted in 990,000 (Equation (6)) estimates of a and b values for every cross-section.
(n 2 − n) × 100 (6) where n is the number of cross-sections (100 in our study). Then, the median a and b values for every cross-section were calculated and applied as the final optimized constants. The reach discharge was calculated by estimating the discharge for every transect (considering the optimized a and b and the water width); then, the median of the estimated discharges for every transect (100 discharge value) was considered as the reach discharge. The developed GA code in our study in Python is available at this respiratory: AMohsenMetwaly/Genetic-Algorithm (github.com; accessed on 2 January 2022).

Suspended Sediment Estimation
The ANN, SVM, and RF algorithms and their combination were applied to produce suspended sediment concentration (SSC) models for the Tisza and Maros Rivers. The models for the Tisza were built based on the in situ measurements of suspended sediment at Szeged and the corresponding 29 Sentinel-2 satellite images. Similarly, the Maros models were built based on the suspended sediment measurements at Makó and the corresponding 27 satellite images. The suspended sediment concentrations were available at five verticals (red or green lines; Figure 3). The whole cross-section was divided into five segments, in such a way that the location of the suspended sediment vertical (the midpoint between the high (red vertical) and low (green vertical) stages) lies approximately at the middle of every segment ( Figure 3). Then, the average of the reflectance values of the images at these five segments was extracted and correlated with its corresponding suspended sediment concentration using ANN, SVM, and RF algorithms. Bands 2, 3, 4, 5, 6, 7, 8, and 8a of the Sentinel-2 images were considered as independent variables. The reflectance and SSC data were divided; 80% of the data were selected for calibration and 20% for validation by the train-test-split function in the scikit-learn library.
Concerning the selected hyperparameters for every model, in the artificial neural network (ANN) model, one input layer containing the eight selected bands (8 neurons), one hidden layer containing 50 neurons, and one output layer were considered. The "relu" activation function and the "lbfgs" solver were applied for the hidden layer, and the maximum iterations were set to 20. For the support vector machine (SVM) model, the polynomial (fourth degree) kernel was applied; moreover, the ε was set to 0.1 and the regularization parameter (c) to 0.7. In the random forest (RF) model, the number of estimators was set to 600, the maximum feature to auto mode, the minimum sample split to 2, and the minimum sample leaf to 4. The criterion used to control the quality of the split was the mean square error (mse). The average of these ANN, SVM, and RF models was also calculated and considered a combined model.

Temporal Change of Suspended Sediment Discharge and Performance Evaluation
The suspended sediment discharge can be obtained by multiplying water discharge and SSC [85] (Equation (7)).
where Q s is the suspended sediment discharge (ton/day), k is a constant referring to units' conversion (0.0864 SI units), Q is water discharge (m 3 /s), and SSC is the suspended sediment concentration (mg/L). The best derived water discharge and SSC models were employed to produce a relatively long time series of the suspended sediment discharge at the three stations (Szeged, Algyő, and Makó). This was performed by applying these models for ungauged dates (122 images) covering the period between July 2015 and May 2021. The daily water discharge data (estimated by rating curves produced by the Lower Tisza Hydrological Directorate (ATIVIZIG) based on daily water stage measurements) at the three stations were used to evaluate the estimated discharge produced by the best discharge model.
The evaluation process between the predicted and measured datasets was performed by three statistical indicators: the coefficient of determination (R 2 ), root mean square error (RMSE), and average distance (AD) from the regression line (Equations (8)-(10)).
where y i is the observed value,ŷ i is the predicted value by a given model, y is the mean of the observed values, and N is the total number of observations.

Temporal Change of Water Discharge and Suspended Sediment Concentration
The rating curve for the Algyő fluviometer (Tisza) was produced based on the daily water stage and discharge data between January 2020 and May 2021, with a coefficient of determination (R 2 ) of 0.98 and RMSE of 62 m 3 /s (Figure 5a). The derived power-law equation was employed to estimate the discharge at Algyő for the ungauged periods (February 2015-December 2019) based on the measured daily water stages. The temporal change of the water discharge and the SSC (2015-2020) based on the monthly in situ measurements is depicted for the Tisza (Figure 5b) and Maros Rivers (Figure 6a). The average discharge of the Tisza at Algyő (575.2 m 3 /s) is always lower than at Szeged (702.6 m 3 /s), as the Szeged gauging station is located downstream of the confluence of the Tisza and the Maros. These discharge values are approximately four times greater than the mean discharge of the Maros River (Makó:157.3 m 3 /s). However, the average SSC of the Maros River (129.14 mg/L) was approximately two-fold the that of the Tisza River (58.5 mg/L).
The largest flood during the studied period occurred in 2016. The Tisza had its peak discharge in March (Algyő: 1640 m 3 /s; Szeged:1760 m 3 /s), but the peak on the Maros (465 m 3 /s) was recorded in June. Usually, the early spring and summer floods are accompanied by high suspended sediment concentrations; however, their amount can be very different. In addition, the SSC peak in the rivers may lag before or after the discharge peak. On the Tisza, for example, during the early spring flood in 2016, the SSC peak was 248 mg/L (Szeged), but during the early summer flood in 2019, it was only 138.4 mg/L. Another typical phenomenon is that the highest SSC is not necessarily concurrent with the discharge peak. For instance, during the early spring flood of 2016, the SSC reached its peak in February, though the discharge peak was measured in March (Figure 5c,d). A similar temporal pattern (hysteresis) was observed in the early summer flood of 2019. On the other hand, on the Maros River, the highest SSC and the discharge peak are likely to occur simultaneously (Figure 6a), as it occurred in the June 2016 flood. This was also confirmed by the Pearson correlation test, as the correlation coefficient (r) between the average SSC and the discharge was 0.88 for the Maros River, and it was only 0.72 for the Tisza River.
Usually, both rivers have low stages and discharge between August and December. In the Tisza River, the lowest water discharge of the studied period was 141 m 3 /s (Szeged, August 2015), and it was only 50.6 m 3 /s in the Maros (Makó, October 2019). In low-stage periods, the SSC usually drops significantly, reaching its lowest values in both rivers, as in October 2019, the Tisza carried only 8.28 mg/L suspended sediment, which was very similar to the data of the Maros (8.19 mg/L).
The in situ SCC measurements were carried out along cross-sections ( Figure 3); thus, they provide data on the spatio-temporal changes of the SSC distribution within the channel (Figures 5b and 6a). The data reflect that in the cross-section of the Maros at Makó, the SSC had a more homogeneous distribution than the Tisza River at Szeged. In the case of the Maros, the suspended sediment concentrations at the five verticals were usually more similar than in the Tisza River. However, in the Maros, slight differences could be observed during floods (e.g., June 2016, March 2017, and February 2019) but without any dominant pattern. For instance, sometimes, along the right bank, there were higher SSC than along the left bank (e.g., February 2019) or vice versa (e.g., April 2018). However, the highest SSC concentrations usually occur in the middle, which is related to the highest velocity and stream power. On the other hand, along the cross-section of the Tisza River, a moderate to significant contrast was found in the SSC, especially during floods (e.g., March 2016, March 2017, and April 2018). Unlike the Maros River, the Tisza River showed a dominant SSC pattern: along the left bank, usually a higher SSC was measured than along the right bank. This is related to the confluence location upstream of the Szeged gauging station.

Estimating Water Discharge
To predict water discharge based on the width measured by Sentinel-2 images, the AHG power-law, AHG machine-learning (ML), and AMHG water discharge models were developed; additionally, their results were compared to define the most suitable method for estimating sediment discharge.

Estimating Water Discharge
To predict water discharge based on the width measured by Sentinel-2 images, the AHG power-law, AHG machine-learning (ML), and AMHG water discharge models were developed; additionally, their results were compared to define the most suitable method for estimating sediment discharge.

Estimating Water Discharge
To predict water discharge based on the width measured by Sentinel-2 images, the AHG power-law, AHG machine-learning (ML), and AMHG water discharge models were developed; additionally, their results were compared to define the most suitable method for estimating sediment discharge.

AHG Power-Law Models
The average width of the five transects ( Figure 3) at a particular date and the corresponding measured water discharge were fitted by power-law equations for Szeged, Algyő, and Makó gauging stations (see Appendix A, Figure A1). The estimated discharges by the AHG power-law models were compared to the measured discharges, and the R 2 , RMSE, and AD values were calculated (Figure 7). The performance of the AHG power-law model at the three stations was similar; however, it had a slightly better performance at the Szeged station (R 2 = 0.6) than at Algyő (R 2 = 0.57) or Makó (R 2 = 0.54) stations. Usually, the models gave good predictions for the discharges below the bankfull (Szeged: 1628 m 3 /s, Algyő: 1600 m 3 /s, and Makó: 412 m 3 /s). However, worse predictions were obtained for discharges close or exceeding the bankfull discharge (e.g., floods in March 2017 at Szeged and Algyő, and in June 2019 at Makó), where the models significantly overestimated the discharges. Although the AHG power-law models of Szeged and Algyő stations had the highest R 2 , they produced high RMSE (Szeged: 754 m 3 /s; RMSE/mean observed: 1.07; Algyő: 1072 m 3 /s; RMSE/mean observed: 1.8) due to their significant overestimation of high discharges (e.g., March 2017 and July 2020). The average width of the five transects ( Figure 3) at a particular date and the corresponding measured water discharge were fitted by power-law equations for Szeged, Algyő, and Makó gauging stations (see Appendix A, Figure A1). The estimated discharges by the AHG power-law models were compared to the measured discharges, and the R 2 , RMSE, and AD values were calculated (Figure 7). The performance of the AHG powerlaw model at the three stations was similar; however, it had a slightly better performance at the Szeged station (R 2 = 0.6) than at Algyő (R 2 = 0.57) or Makó (R 2 = 0.54) stations. Usually, the models gave good predictions for the discharges below the bankfull (Szeged: 1628 m 3 /s, Algyő: 1600 m 3 /s, and Makó: 412 m 3 /s). However, worse predictions were obtained for discharges close or exceeding the bankfull discharge (e.g., floods in March 2017 at Szeged and Algyő, and in June 2019 at Makó), where the models significantly overestimated the discharges. Although the AHG power-law models of Szeged and Algyő stations had the highest R 2 , they produced high RMSE (Szeged: 754 m 3 /s; RMSE/mean observed: 1.07; Algyő: 1072 m 3 /s; RMSE/mean observed: 1.8) due to their significant overestimation of high discharges (e.g., March 2017 and July 2020).

AHG Machine-Learning Models
The SVR, RF, and ANN algorithms were calibrated by 80% of the data and validated by 20%. A combined model considering the average of their predictions was also derived. The R 2 and RMSE for every model at the three stations were produced based on the validation dataset ( Table 1). The performance of the machine-learning models (SVR, RF, and ANN) was promising, as the R 2 ranged between 0.65 (SVR model-Szeged station) and

AHG Machine-Learning Models
The SVR, RF, and ANN algorithms were calibrated by 80% of the data and validated by 20%. A combined model considering the average of their predictions was also derived. The R 2 and RMSE for every model at the three stations were produced based on the validation dataset ( Table 1). The performance of the machine-learning models (SVR, RF, and ANN) was promising, as the R 2 ranged between 0.65 (SVR model-Szeged station) and 0.88 (combined model-Makó station); moreover, the RMSE for Tisza models ranged between 154.12 m 3 /s (RF model-Szeged station) and 98.92 m 3 /s (SVR model-Algyő station), and for Maros models between 33.65 m 3 /s (SVR model-Makó station) and 25.37 m 3 /s (RF model-Makó station). The models with the highest R 2 were considered the best performing models for every gauging station; thus, the combined model was considered the best model for Szeged (R 2 = 0.83) and Makó (R 2 = 0.88) stations, and the ANN for Algyő (R 2 = 0.79) station. Table 1. The results of the evaluation metrics (R 2 and RMSE) for estimating water discharge using support vector regression (SVR), random forest (RF), artificial neural network (ANN) algorithms and their combined model based on the validation dataset (20%) for Szeged, Algyő, and Makó stations. The models in bold refer to the best AHG machine-learning model at every station based on the R 2 values. The R 2 and RMSE of the AHG power-law and AMHG methods were added for the purpose of comparison. The best performing models were applied to predict the discharge for the entire dataset (including the calibration and validation data) at every gauging station (Szeged and Algyő: 29 images; Makó: 27 images). This enabled the comparison of the various adopted methods (i.e., AHG power-law, AHG machine-learning (ML), and AMHG methods) to retrieve water discharge. The predicted and actually measured discharges, in addition to the R 2 , RMSE, and AD data, are depicted in Figure 8. The estimated discharges by the AHG ML models were more consistent with the measured discharges than the AHG power-law models. This was especially the case with high discharges, which were usually overestimated by the AHG power-law models; meanwhile, the AHG ML models gave closer estimations to the measured discharges. This was also reflected in the RMSE results, as in the case of Algyő station it dropped by 80% (from 1072 m 3 /s to 181 m 3 /s), at Szeged by 75% (from 754 m 3 /s to 192 m 3 /s), and by 58% (from 110 m 3 /s to 46 m 3 /s) in case of Makó station compared to the AHG power-law models. Based on the evaluation metrices (Figure 8d), it can be noted that the AHG ML models performed best at Szeged station, followed by Makó and Algyő stations.

AMHG Models
The "at-many-stations hydraulic geometry" (AMHG) method was also tested based on the 100 transects at every gauging station to estimate water discharge. The relationship between the logarithm of the coefficient log (a) and the exponent b based on the 100 transects at every gauging station, and their optimized values by the genetic algorithm (GA), are depicted ( Figure 9). All the gauging stations showed a log-linear relationship between log (a) and b. However, the best relationship was achieved at Szeged (R 2 = 0.88), followed by Algyő (R 2 = 0.86) stations, and the worst was for the Makó station (R 2 = 0.47). The optimized coefficients and exponents of the 100 transects at every station seem to be clustered around approximately the middle of the log-linear regression line.   Table 1 applied on the calibration and validation dataset. The R 2 , RMSE, AD values and the linear fitting line are also indicated for every station (d).

AMHG Models
The "at-many-stations hydraulic geometry" (AMHG) method was also tested based on the 100 transects at every gauging station to estimate water discharge. The relationship between the logarithm of the coefficient log (a) and the exponent b based on the 100 transects at every gauging station, and their optimized values by the genetic algorithm (GA), are depicted ( Figure 9). All the gauging stations showed a log-linear relationship between log (a) and b. However, the best relationship was achieved at Szeged (R 2 = 0.88), followed by Algyő (R 2 = 0.86) stations, and the worst was for the Makó station (R 2 = 0.47). The optimized coefficients and exponents of the 100 transects at every station seem to be clustered around approximately the middle of the log-linear regression line.   Table 1 applied on the calibration and validation dataset. The R 2 , RMSE, AD values and the linear fitting line are also indicated for every station (d). The estimated discharges by the AMHG method were compared to the measured discharges, and the R 2 and RMSE for every station were represented ( Figure 10). The performance of the AMHG models outweighed the AHG power-law models, as the R 2 was improved by ≈5%, on average, and the RMSE was declined by ≈65%, on average. Mean-  The estimated discharges by the AMHG method were compared to the measured discharges, and the R 2 and RMSE for every station were represented ( Figure 10). The performance of the AMHG models outweighed the AHG power-law models, as the R 2 was improved by ≈5%, on average, and the RMSE was declined by ≈65%, on average. Meanwhile, their performance was closer or slightly lower than the AHG ML models, as the R 2 was 8% (on average) lower than that in the AHG ML models (average R 2 = 0.7), and the RMSE was 19% higher (AHG ML: average RMSE: 139.67 m 3 /s). Although the AMHG models overcame the significant overestimation of the AHG power-law models for high discharges, they usually underestimated the real discharges (e.g., March 2017 and June 2019 at Szeged and Algyő stations; June 2016 and June 2019 at Makó station) but with lower rates. The AMHG gave good predictions for low and medium discharges. The AMHG model at the Szeged station was the best (R 2 = 0.67), followed by Algyő (R 2 = 0.67) and Makó (R 2 = 0.55) stations. The same order was observed in the case of the AHG power-law models.

Estimating Suspended Sediment Concentration
The spectral signatures of water in the Tisza (Szeged station) and Maros (Makó station) Rivers were analyzed to estimate their suspended sediment concentrations. The spectral signature based on the average of reflectance during low and high stages is depicted (Figure 11a). Furthermore, it was also represented for low and high stages separately for both rivers, as it is expected that they would produce different signatures due to the significant variations in SSC during these hydrological situations (Figure 11b). Gen- Based on the performance in terms of R 2 , RMSE, and AD of the AHG power-law, AHG ML, and AMHG models at the three gauging stations, it can be noticed that the AHG ML and AMHG were the most reliable methods to retrieve the discharge, while the AHG power-law was the least accurate method. Usually, the water discharge estimations at Szeged were better than at Algyő, and they had the lowest quality at Makó. This could be inferred from the gradient of R 2 and AD, as they showed the best values at Szeged and the worst at Makó (Figures 7-10). It reveals the significant role of the river channel size and cross-sectional shape in the applied water discharge estimation methods.

Estimating Suspended Sediment Concentration
The spectral signatures of water in the Tisza (Szeged station) and Maros (Makó station) Rivers were analyzed to estimate their suspended sediment concentrations. The spectral signature based on the average of reflectance during low and high stages is depicted (Figure 11a). Furthermore, it was also represented for low and high stages separately for both rivers, as it is expected that they would produce different signatures due to the significant variations in SSC during these hydrological situations (Figure 11b). Generally, the reflectance of the Maros River was higher than the Tisza River; furthermore, the reflectance at high stages in both rivers was higher than at low stages (Figure 11a,b). Both rivers showed low reflectance at the blue band (B2), which refers to the turbidity of water, as clear water has high reflectance at this band [2]. Furthermore, a decreasing trend in the reflectance can be noticed between the red edge 2 (B6) and SWIR 2 (B12) bands. The highest reflectance occurred between the green (B3) and red edge 1 (B5) bands. The Tisza River showed the highest reflectance at the green (B3) band, while that in the Maros River was at the red edge 1 (B5) band (Figure 11a). Similarly, the green (B3) band had had the highest reflectance at the low stages, and the red edge 1 (B5) at the high stages of both rivers, which indicates the usefulness of the green (B3) band for detecting low SSC and red edge 1 (B5) for high SSC. The R 2 , RMSE, and RMSE/mean observed of the SVR, RF, ANN, and combined SSC models for the Tisza and Maros Rivers were tabulated (Table 2). Generally, all the derived models showed good results, as the R 2 ranged between 0.76 and 0.9, and the RMSE was between 24.9 mg/L and 15.23 mg/L. However, the SSC models of the Tisza River showed slightly lower accuracies than the Maros models, as the average R 2 of the Tisza models (0.8) was lower than the Maros models (0.85) by 5% (based on the average of all derived models). However, as the Maros River had greater SSC, the average RMSE of the Maros (21.8 mg/L) was slightly higher than the Tisza (16.47 mg/L). The best derived SSC model for the Tisza River was the combined model, followed by the ANN, SVR, and RF models. Meanwhile, the RF was the best in the Maros River, followed by the combined model, ANN, and SVR models. Hence, the combined SSC model in the Tisza and the RF model in the Maros were used to predict the sediment discharge (Qs).
The measured and predicted SSC by the best models were compared for both rivers based on the validation dataset ( Figure 12). The predicted SSC was quite consistent with the measured one for both rivers, especially for the Maros River, as the R 2 in the Maros was higher than the Tisza by 8% (based on the best derived models for both rivers). Despite the accurate estimations of the derived models, both models tended to slightly un- The R 2 , RMSE, and RMSE/mean observed of the SVR, RF, ANN, and combined SSC models for the Tisza and Maros Rivers were tabulated (Table 2). Generally, all the derived models showed good results, as the R 2 ranged between 0.76 and 0.9, and the RMSE was between 24.9 mg/L and 15.23 mg/L. However, the SSC models of the Tisza River showed slightly lower accuracies than the Maros models, as the average R 2 of the Tisza models (0.8) was lower than the Maros models (0.85) by 5% (based on the average of all derived models). However, as the Maros River had greater SSC, the average RMSE of the Maros (21.8 mg/L) was slightly higher than the Tisza (16.47 mg/L). The best derived SSC model for the Tisza River was the combined model, followed by the ANN, SVR, and RF models. Meanwhile, the RF was the best in the Maros River, followed by the combined model, ANN, and SVR models. Hence, the combined SSC model in the Tisza and the RF model in the Maros were used to predict the sediment discharge (Q s ). Table 2. The R 2 , RMSE, and RMSE/mean observed for the applied machine-learning models at the Tisza and Maros Rivers for estimating suspended sediment concentration (SSC) by support vector regression (SVR), random forest (RF), artificial neural network (ANN) algorithms based on the validation dataset (20%). The bold fonts refer to the best performing model for a given river based on the R 2 value of the applied model. The measured and predicted SSC by the best models were compared for both rivers based on the validation dataset ( Figure 12). The predicted SSC was quite consistent with the measured one for both rivers, especially for the Maros River, as the R 2 in the Maros was higher than the Tisza by 8% (based on the best derived models for both rivers). Despite the accurate estimations of the derived models, both models tended to slightly underestimate the actual measured SSC, as the mean estimated SSC based on the validation dataset in the Tisza was 38.5 mg/L, contrary to 40.7 mg/L mean measured SSC (percent bias: −5.31%). Similarly, the mean estimated SSC in the Maros was 86.6 mg/L, contrary to the 89.5 mg/L mean measured SSC (percent bias: −3.18%). Based on the best performing models, the SSC for ca. 20 km long sections of the Tisza and Maros Rivers were calculated. As an example, the spatial distributions of the SSC during high (5 March 2020) and low stages (6 September 2020) were presented ( Figure  13a,b). During the high-stage situation (Figure 13a), the average SSC of the Tisza (106.22 mg/L) was 6.3 times greater than at the low stage (16.96 mg/L) (Figure 13b). Similarly, the average SSC of the Maros at the high stage (211 mg/L) was 5.2 times greater than at the low stage (40.8 mg/L). In both hydrological situations, the Maros River had higher SSC than the Tisza River; however, the concentration difference decreased during the high stages, as on the given days, the Maros (211 mg/L) had 2 times higher concentrations than the Tisza (106.22 mg/L), and their difference increased to 2.4 times at low stage (Maros: 40.8 mg/L; Tisza: 16.96 mg/L). The spatial distribution of the SSC was more homogeneous along the Tisza River's cross-sections during low stages than the high stages. For instance, during high stages, the highest SSC could be noticed at the middle of the Algyő gauging station ( Figure  13(aI)), and it decreased gradually toward the riverbanks; in the meantime, it had almost the same concentrations along the entire cross-section during low stages (Figure 13(bI)). Concerning the Maros River, the SSC was almost the same along the cross-sections during high or low stages (Figure 13(aIV,aV)), Figure 13(bIV,bV); however, the SSC of some pixels along the riverbanks were underestimated due to the tree shadow, which had a drastic impact on this river due to its limited width (average width 126 m). The SSC spatial distribution at the confluence area showed a more complicated spatial distribution than the rivers themselves, as, during high stages, the Maros is usually impounded by the Tisza due to higher momentum of the Tisza. Thus, it extended longitudinally with limited width (Figure 13(aII,aIII); meanwhile, during low stages, the Maros water extended laterally wider than at the high stages (Figure 13(bII,bIII)). However, it mixed faster, as it showed a very limited longitudinal extent compared to the high stages. Based on the best performing models, the SSC for ca. 20 km long sections of the Tisza and Maros Rivers were calculated. As an example, the spatial distributions of the SSC during high (5 March 2020) and low stages (6 September 2020) were presented (Figure 13a,b). During the high-stage situation (Figure 13a), the average SSC of the Tisza (106.22 mg/L) was 6.3 times greater than at the low stage (16.96 mg/L) (Figure 13b). Similarly, the average SSC of the Maros at the high stage (211 mg/L) was 5.2 times greater than at the low stage (40.8 mg/L). In both hydrological situations, the Maros River had higher SSC than the Tisza River; however, the concentration difference decreased during the high stages, as on the given days, the Maros (211 mg/L) had 2 times higher concentrations than the Tisza The spatial distribution of the SSC was more homogeneous along the Tisza River's cross-sections during low stages than the high stages. For instance, during high stages, the highest SSC could be noticed at the middle of the Algyő gauging station (Figure 13(aI)), and it decreased gradually toward the riverbanks; in the meantime, it had almost the same concentrations along the entire cross-section during low stages (Figure 13(bI)). Concerning the Maros River, the SSC was almost the same along the cross-sections during high or low stages (Figure 13(aIV,aV)), Figure 13(bIV,bV); however, the SSC of some pixels along the riverbanks were underestimated due to the tree shadow, which had a drastic impact on this river due to its limited width (average width 126 m). The SSC spatial distribution at the confluence area showed a more complicated spatial distribution than the rivers themselves, as, during high stages, the Maros is usually impounded by the Tisza due to higher momentum of the Tisza. Thus, it extended longitudinally with limited width (Figure 13(aII,aIII); meanwhile, during low stages, the Maros water extended laterally wider than at the high stages ( Figure 13(bII,bIII)). However, it mixed faster, as it showed a very limited longitudinal extent compared to the high stages.

Estimated Suspended Sediment Discharge of the Tisza and Maros Rivers Based on the Derived Models
After evaluating the performance of the three applied methods (i.e., AHG powerlaw, AHG ML, and AMHG) to estimate the water discharge from the Sentinel-2 satellite images, the best performing AHG ML models were applied to predict the water discharge at Szeged, Algyő, and Makó gauging stations between July 2015 and May 2021. Similarly, the combined SSC model was utilized to predict the SSC in the Tisza River (Algyő and Szeged stations), whereas the RF was applied to estimate the SSC of the Maros River (Makó station).
A time series of the suspended sediment discharge (Qs) at Algyő, Szeged, and Makó stations between July 2015 and May 2021 based on 122 Sentinel-2 images was created ( Figure 14). Although the SSC in the Maros River was usually 2-3 times greater than in the Tisza River, the suspended sediment discharge of the Tisza River was 2-3.5 times greater than that of the Maros River. This is because of the elevated discharge of the Tisza River (four-fold) compared to the Maros River. Usually, the suspended sediment discharge of the Tisza upstream of the confluence (at Algyő average Qs= 2584 t/d) was greater than that of the Maros (at Makó average Qs = 1264 t/d). However, downstream of their confluence, the suspended sediment discharge at Szeged increased considerably (Qs= 4385 t/d). On the other hand, in some hydrological situations, the Qs at Makó was greater than the Qs at Algyő (e.g., June and August 2016, July 2017, July-August 2018, June-July 2019). This usually occurs during low stages, when the SSC in the Maros River increases due to small flood waves compared to the very low SSC in the Tisza River.
Usually, the temporal changes in suspended sediment discharge followed a similar pattern as water discharge; thus, the sediment discharge peaks were observed during

Estimated Suspended Sediment Discharge of the Tisza and Maros Rivers Based on the Derived Models
After evaluating the performance of the three applied methods (i.e., AHG powerlaw, AHG ML, and AMHG) to estimate the water discharge from the Sentinel-2 satellite images, the best performing AHG ML models were applied to predict the water discharge at Szeged, Algyő, and Makó gauging stations between July 2015 and May 2021. Similarly, the combined SSC model was utilized to predict the SSC in the Tisza River (Algyő and Szeged stations), whereas the RF was applied to estimate the SSC of the Maros River (Makó station).
A time series of the suspended sediment discharge (Q s ) at Algyő, Szeged, and Makó stations between July 2015 and May 2021 based on 122 Sentinel-2 images was created ( Figure 14). Although the SSC in the Maros River was usually 2-3 times greater than in the Tisza River, the suspended sediment discharge of the Tisza River was 2-3.5 times greater than that of the Maros River. This is because of the elevated discharge of the Tisza River (four-fold) compared to the Maros River. Usually, the suspended sediment discharge of the Tisza upstream of the confluence (at Algyő average Q s = 2584 t/d) was greater than that of the Maros (at Makó average Q s = 1264 t/d). However, downstream of their confluence, the suspended sediment discharge at Szeged increased considerably (Q s = 4385 t/d). On the other hand, in some hydrological situations, the Q s at Makó was greater than the Q s at Algyő (e.g., June and August 2016, July 2017, July-August 2018, June-July 2019). This usually occurs during low stages, when the SSC in the Maros River increases due to small flood waves compared to the very low SSC in the Tisza River.

Discussion
Although sediment discharge monitoring in rivers is an important issue from several aspects (hydrology, geomorphology, ecology, contamination control), many rivers are still ungauged, and the number of gauged rivers is being decreased because these measurements are labor intensive and expensive [19]. In the meantime, remote sensing with its large-scale coverage and frequent revisiting could be employed to monitor rivers effectively. Therefore, this study aimed to estimate suspended sediment discharge by Sentinel-2 images, considering the medium-sized Tisza River and its tributary Maros River in central Europe as a case study. Therefore, various methods were applied to estimate the water discharge and the suspended sediment concentration (SSC); thus, the suspended sediment discharge could be calculated.

Suspended Sediment and Water Discharge Conditions in the Rivers Based on In Situ Measurements
The measured SSC in the Maros was twice as much as that in the Tisza River upstream of their confluence because the catchment of the Maros consists of highly erodible rocks, and there are intensive mining activities in Transylvania [86]. The temporal change of the SSC was in good agreement with the water discharge (Tisza River: r = 0.88; Maros River: r = 0.72), as the highest SSC were detected during floods and the lowest during low stages (Figures 5b and 6a). However, a hysteretic behavior in the relationship between the SSC and Q was observed in both rivers, especially in the Tisza River; the SSC peak usually precedes the water discharge peak (clockwise hysteresis). This was also observed at the Middle Tisza (Szolnok and Kisköre) based on five years (1998)(1999)(2000)(2001)(2002) of measurements of the SSC and water discharge conducted by Csépes et al. [69]. A similar phenomenon was reported in many rivers, such as the Save River, France [87], Rybárik Basin, Slovakia [88], Chabagou Creek, China [89], and Holbeck Basin, England [90]. This clockwise (positive) hysteresis could develop due to the availability of sediments in the channel and in the neighboring areas at the beginning of the flood wave [87], which can be reflected in a higher SSC in the rising limb of the storm hydrograph than the falling limb.
The lateral distribution of the measured SSC at the cross-section of the rivers usually corresponds to the velocity distribution [91]. This was the case for the Makó and Algyő Usually, the temporal changes in suspended sediment discharge followed a similar pattern as water discharge; thus, the sediment discharge peaks were observed during floods and the lowest values during low stages. Usually, the Q s peak at the Maros River was concurrent with that in the Tisza River; however, in some cases, it preceded (e.g., February 2019 and February 2021) or followed (e.g., June 2016 and March 2017) the peak Q s of the Tisza. As the Szeged station is located on the Tisza River downstream of the confluence of the Tisza and the Maros, here, sediment discharge differed significantly compared to the upstream sections. As a result, sometimes, the Q s peak at Szeged station shifted forward (e.g., March 2017) or backward (e.g., February 2021), the Q s at Algyő following the Q s of the Maros River.

Discussion
Although sediment discharge monitoring in rivers is an important issue from several aspects (hydrology, geomorphology, ecology, contamination control), many rivers are still ungauged, and the number of gauged rivers is being decreased because these measurements are labor intensive and expensive [19]. In the meantime, remote sensing with its largescale coverage and frequent revisiting could be employed to monitor rivers effectively. Therefore, this study aimed to estimate suspended sediment discharge by Sentinel-2 images, considering the medium-sized Tisza River and its tributary Maros River in central Europe as a case study. Therefore, various methods were applied to estimate the water discharge and the suspended sediment concentration (SSC); thus, the suspended sediment discharge could be calculated.

Suspended Sediment and Water Discharge Conditions in the Rivers Based on In Situ Measurements
The measured SSC in the Maros was twice as much as that in the Tisza River upstream of their confluence because the catchment of the Maros consists of highly erodible rocks, and there are intensive mining activities in Transylvania [86]. The temporal change of the SSC was in good agreement with the water discharge (Tisza River: r = 0.88; Maros River: r = 0.72), as the highest SSC were detected during floods and the lowest during low stages (Figures 5b and 6a). However, a hysteretic behavior in the relationship between the SSC and Q was observed in both rivers, especially in the Tisza River; the SSC peak usually precedes the water discharge peak (clockwise hysteresis). This was also observed at the Middle Tisza (Szolnok and Kisköre) based on five years (1998)(1999)(2000)(2001)(2002) of measurements of the SSC and water discharge conducted by Csépes et al. [69]. A similar phenomenon was reported in many rivers, such as the Save River, France [87], Rybárik Basin, Slovakia [88], Chabagou Creek, China [89], and Holbeck Basin, England [90]. This clockwise (positive) hysteresis could develop due to the availability of sediments in the channel and in the neighboring areas at the beginning of the flood wave [87], which can be reflected in a higher SSC in the rising limb of the storm hydrograph than the falling limb.
The lateral distribution of the measured SSC at the cross-section of the rivers usually corresponds to the velocity distribution [91]. This was the case for the Makó and Algyő gauging stations, thus upstream of the confluence. However, downstream of the confluence of the Tisza and Maros Rivers (at Szeged), the SSC was the highest at the low velocity left bank, and it decreased toward the right (concave) bank. This unusual situation occurred due to the special spatial configuration of the gauging station just 4.1 km far from the confluence. As the Maros joins from the left side, usually with higher SSC, at the gauging station, the waters of both rivers are still mixing, as the complete mixing requires a much longer distance (e.g., 15-150 km, according to Jirka [92]). This result is consistent with Mohsen et al. [47], who investigated the spatio-temporal dynamism of the mixing process in the confluence of the Tisza and Maros Rivers, reporting that the water of the Maros could be extended to the Szeged gauging station in certain hydrological situations (i.e., Q Makó /Q Szeged = 0.28-0.51; and their slope difference is ∆S Tisza-Maros = −17.7-8.2).

Performance of the Water Discharge Models
According to our study, the combination of the hydraulic geometry theory and remote sensing was very useful, as following the best described method, the water discharge could be estimated with an accuracy of Tisza: 186.5 m 3 /s (RMSE/mean observed = 0.29) and Maros: 46 m 3 /s (RMSE/mean observed = 0.29). Thus, it could be a useful tool to estimate the water discharge along ungauged sections. However, selecting the appropriate reach is critical for obtaining representative discharges. This was also reported by Flores et al. [2], as they applied the Bayesian AMHG method to obtain water discharge at various reaches (i.e., braided, relatively straight, and low sinuosity); the accuracy varied, but the best accuracy was reached at relatively straight and low sinuosity reaches. Although the three tested reaches in our study were relatively straight (Figure 3 In addition, the Tisza has V-shape-like cross-sections, whereas the Maros has a trapezoidal or a near rectangular shape (Figure 3). From the point of view of the width-discharge relationship, the V shape channels provide better hydraulic geometry than rectangular channels [12]; therefore, the water discharge estimation could be more precise. On the other hand, the bankfull level is considered one of the main limits of the hydraulic geometry theory (discharge-width relationship), as the geometrical characteristics of the cross-section of the channel change above this level due to covering the floodplains. Therefore, the derived models usually underestimate (e.g., AHG ML and AMHG method) or overestimate (e.g., AHG power-law method) its actual value, regardless of the selected cross-section, due to losing the hydraulic geometry of the original cross-section.
The performance of the AMHG method outweighed the AHG power-law method, as the average R 2 for the AMHG models was 0.62 in contrast to 0.57 of the AHG power-law models (based on the three tested reaches). Similarly, the RMSE dropped from 172 m 3 /s in cased of AMHG models to 645 m 3 /s in case of AHG power-law models (73% difference). Similar results were reported by Durga Rao et al. [35], who applied the AHG powerlaw and AMHG methods to estimate water discharge for four main rivers in India using Landsat and Indian remote sensing satellites (IRS) images. According to their data, the Nash coefficients for the AHG power-law and AMHG models were 0.696 and 0.803, respectively. This difference could be explained by the usage of the hydraulic geometry of many crosssections along the reach in the case of the AMHG method, while only one cross-section was considered in the AHG power-law method. In addition, the AHG method is considered an inherent part of the AMHG method. We developed a novel method by combining the AHG method with machine-learning algorithms, and it gave promising results on water discharge estimation. The results of the AHG ML models were comparable or even better than the AMHG models, as the average R 2 for the three gauging stations based on the AHG ML method was 0.7 in comparison to 0.62 in the case of the AMHG method; moreover, the RMSE dropped from 172 m 3 /s (AMHG) to 140 m 3 /s (AHG ML). The AHG ML models were very stable, as they showed good estimations when they were applied to estimate water discharge for the time series of 122 images (2015-2021), which were totally new data or unseen during the training process (see Appendix A, Figure A2).
The superiority of the AHG ML method can be interpreted by the nature of machine learning itself, as it was mainly developed to give the best possible predictions by employing general-purpose learning algorithms to uncover patterns in data that are frequently large and unmanageable [93]. Unlike the conventional statical regression methods, such as the AHG power-law method, the ML algorithms usually give ambiguous structures, as they learn from the provided dataset [93]. Based on this finding, the AHG ML method can substitute the hydraulic geometry of the multiple cross-sections in the AMHG method; thus, it can avoid the complex optimization process of the coefficients and components in the AMHG, and consequently, larger cross-sections along the river can be monitored.
It is important to note that the accuracy of the derived water width from the satellite images affects the estimated water discharge by the hydraulic geometry method significantly. This mainly depends on the spatial resolution of the images and the extraction method of water pixels. The automatic thresholding of the NDWI using the OTSU method proved its ability as a good tool to extract the river water pixels effectively, as the RMSE for the estimated water width by the satellite images at Szeged gauging station was 8.17 m and 9.9 at Makó gauging station, which was very close to the actual pixel size of the Sentinel-2 images. Thus, applying this method to finer spatial resolution images would produce good estimates for the water width of a river.

Performance of the Suspended Sediment Models
Based on the spectral signatures of water of the studied rivers, the green (B3), red (B4), and red edge 1 (B5) bands seem to be the most appropriate bands for retrieving SSC from Sentinel-2 images, as they showed the greatest reflectance throughout the spectrum. Similar results were reported by many studies [46,50,[94][95][96]. Machine-learning algorithms produced good SSC models, as the combined model of the Tisza gave R 2 of 0.82, and the RF gave R 2 of 0.9 in the Maros River. The accuracy of the developed models in our study was comparable to other machine-learning-based SSC models developed in other studies. For instance, the extreme learning machine SSC model in the Missouri-Mississippi Rivers system showed an R 2 of 0.9 [56]. Umar et al. [97] developed an SSC model based on an RF algorithm for the confluence of the Missouri-Mississippi Rivers with an R 2 of 0.72. On the other hand, the performance of the machine-learning-based SSC models outweighed the statistical regression models, which reveals the usefulness of machine learning for dealing with complex regression models. For example, Larson et al. [98] produced an SSC linear regression model for the Maumee River with an R 2 of 0.56, and a similar SSC linear regression model was produced for the Mississippi with an R 2 of 0.62.
Although the three tested ML algorithms gave satisfactory results, which was also reported in previous studies [99][100][101], the RF algorithm produced the best SSC models, followed by the ANN and SVR. This finding was consistent with Larson et al. [95]. The superiority of RF and ANN algorithms over the SVR algorithm arises from the nature of these algorithms, as both the RF and ANN approaches compute several models and then average them to create a single model, but the SVR requires extensive parameter adjustment before producing a model [95]. Our study also indicated that the aggregation of more than one algorithm could enhance the estimation accuracy significantly.
Despite the significant temporal dynamism of the hydrological condition and the SSC in the rivers, especially during floods, the selected timeframe of three days between the in situ campaigns and the acquisition data of the images did not affect the accuracy of the derived models significantly. Additionally, averaging the surface reflectance of the pixels within each segment for the five segments at the cross-sections at Szeged and Makó gauging stations could also improve the accuracy of the models by mitigating the effects of sun glint, surface wave, and the slight shifts of the location of the verticals away from the actual sampling location [56].

Suspended Sediment Concentration and Discharge Conditions of the Studied Rivers Based on Sentinel-2 Images
The 3-5 days temporal resolution of the Sentinel-2 images allowed more frequent monitoring (122 images; 2015-2021) of the rivers than the monthly in situ monitoring (59 in situ campaigns); however, cloud cover is the main limitation of this space-based alternative. On the other hand, the availability of the large-scale coverage of the SSC provided by the space-based alternative, which could cover the entire river, is a great advantage, which is unobtainable through the in situ campaigns. The temporal change of the estimated SSC and water discharge was consistent with the measured one. During floods, in most cases (Tisza: 66%; Maros: 66%), the SSC peak is concurrent with the water discharge peak; however, clockwise (Tisza: 22%; Maros: 16%) and anticlockwise (Tisza: 11%; Maros: 16%) hysteresis was also observed. Usually, the SSC peak of both rivers occurs simultaneously. Although the spatial resolution of the Sentinel-2 images is 10 m, thus the Maros was covered by 12 pixels on average and the Tisza by 16 pixels, it was very useful to reveal the spatial dynamism of the SSC along the studied area. Based on these SSC estimations, we disclosed that both rivers have more SSC spatial variability along the rivers' cross-section during floods than at low stages due to the significant amount of the transported suspended load and the elevated stream power. Thus, the highest concentrations were usually found along with the fastest stream current and decreased gradually apart from this line. In the meantime, several mixing patterns were observed in the confluence area, which was studied in detail by Mohsen et al. [47].
Based on our findings, we offer a combination of remote-sensing and machine-learning algorithms to monitor the suspended sediment discharge in rivers. It is a useful tool to capture the spatio-temporal dynamism of the Q s . Although this study was applied on a medium-sized river (Tisza River) and its tributary (Maros River), the derived models produced good estimations of Q s , thus better estimations would be expected with largersized rivers.

Conclusions
Sediment discharge is affected by various parameters of the catchment. However, under the climate change and human impact, the changes in the spatial and temporal pattern of the sediment household could negatively affect the fluvial system and the environment. Sediment discharge is usually monitored by in situ measurements of water discharge and SSC; however, their spatio-temporal limits hinder efficient monitoring. Therefore, in this study, a remote-sensing alternative for sediment discharge monitoring was offered using Sentinel-2 images. The water discharge was estimated by AHG powerlaw, AHG machine-learning, and AMHG methods; additionally, the SSC was estimated by SVC, RF, ANN, and combined models, and their performance was compared. The best performing models were employed to estimate the suspended sediment discharge of two medium-sized rivers (2015-2021).
Based on the results of this study, we concluded that: • The hydraulic geometry of the river channel in terms of the discharge-width relationship is very useful for estimating water discharge from space; however, the morphology of the selected reach (i.e., the shape of the cross-sectional profile) significantly affects the accuracy of the estimation.

•
The automatic thresholding of the NDWI using the OTSU method was very useful, as it estimated the river width with an accuracy close to the pixel size. Thus, good channel width estimations could be obtained by applying this method on finer spatial resolution images, which would enhance the accuracy of water discharge estimations. • Applied in a traditional way, the performance of the AMHG method outweighs the AHG power-law method, as was also reported in the former literature; however, our novel AHG machine-learning method performed as well as the AMHG method due to the capability of machine-learning algorithms for deriving complex regression models. • Among the three tested algorithms for obtaining suspended sediment concentration based on Sentinel-2 images, the RF and ANN algorithms gave the best SSC estimations; thus, they are highly recommended for further studies. We also recommend combining algorithms to enhance the estimation accuracy significantly.

•
The accuracy of the estimated Q s by this space-based alternative is limited by: (a) river width and the spatial resolution of the utilized images; (b) selected cross-section shape, as some shapes have more suitable hydraulic geometry than others; (c) bankfull discharge, as the models usually under-or overestimate the real discharge above this level due to losing the hydraulic geometry of the original cross-section; (d) shadow, sun glint, and surface wave in water pixels affecting the reflectance significantly, and consequently, the predicted SSC.
The combination of remote-sensing images and machine-learning techniques could efficiently monitor the suspended sediment concentration and discharge of rivers, as it could produce robust regression models applied on large-scale and frequent images. Although Sentinel-2 images have a spatial resolution of 10 m, they could provide good estimations on suspended sediment discharge, though finer spatial resolution images could achieve better estimations. Applying this method to many rivers worldwide could improve our understanding not only of the spatio-temporal variations of suspended sediment discharge but also of the transport of suspended contaminants (e.g., microplastic and mining waste), as their transport might be related to the natural sediment load.