Spaceborne River Discharge From a Nonparametric Stochastic Quantile Mapping Function

The number of active gauges with open‐data policy for discharge monitoring along rivers has decreased over the last decades. Therefore, spaceborne measurements are investigated as alternatives. Among different techniques for estimating river discharge from space, developing a rating curve between the ground‐based discharge and spaceborne river water level or width is the most straightforward one. However, this does not always lead to successful results, since the river section morphology often cannot simply be modeled by a limited number of parameters. Moreover, such methods do not deliver a proper estimation of the discharge's uncertainty as a result of the mismodeling and also the coarse assumptions made for the uncertainty of inputs. Here, we propose a nonparametric model for estimating river discharge and its uncertainty from spaceborne river width measurements. The model employs a stochastic quantile mapping scheme by, iteratively: (a) generating realizations of river discharge and width time series using Monte Carlo simulation, (b) obtaining a collection of quantile mapping functions by matching all possible permutations of simulated river discharge and width quantile functions, and (c) adjusting the measurement uncertainties according to the point cloud scatter. We validate our method over 14 different river reaches along the Niger, Congo, Po Rivers, and several river reaches in the Mississippi river basin. Our results show that the proposed algorithm can mitigate the effect of measurement noise and also possible mismodeling. Moreover, the proposed algorithm delivers a meaningful uncertainty for the estimated discharge and allows us to calibrate the error bars of in situ discharge measurements.

10.1029/2021WR030277 2 of 27 been developed for estimating river discharge based on Manning's equation for open-channel flow by employing spaceborne measurements of river water level and width, channel slope, and flow depth. For example, Durand et al. (2014) presented the Metropolis-Manning (MetroMan) algorithm to estimate river bathymetry, the roughness coefficient and discharge based on input measurements of river water level and slope using the Metropolis algorithm in a Bayesian Markov chain Monte Carlo scheme. Garambois and Monnier (2015) implemented Manning's equation to retrieve the river low flow bathymetry, roughness, and discharge. Durand et al. (2016) introduced the Mean-annual Flow and Geomorphology (MFG) algorithm, which uses the so-called wide-channel approximation leading to a form of Manning's equation that approximates river depth as the difference between water surface elevation and the cross-sectional average river bathymetry. Tourian et al. (2017) proposed an algorithm to estimate the average river depth from spaceborne river water level, width, and slope measurements using the two Manning type models for discharge estimation developed by Bjerklie et al. (2003) and Dingman and Sharma (1997).
The upcoming Surface Water and Ocean Topography (SWOT) satellite mission will simultaneously measure the river width, height, and slope for major rivers. The challenge of finding the best model for estimating discharge using these quantities has been addressed in several studies e.g., Durand et al. (2014), Garambois and Monnier (2015), Gleason and Smith (2014), Bjerklie et al. (2018), Tuozzolo, Langhorst, et al. (2019), Tuozzolo, Lind, et al. (2019), Altenau et al. (2017), Bonnema et al. (2016), and Sichangi et al. (2018). The proposed algorithms present different variants of the McFLI paradigm; however, all of them invoke at least one flow-law and mass conservation . All of the aforementioned techniques follow almost the same approach for estimating discharge. The river is divided into several reaches which represent the basic hydraulic unit of the river. Then the river discharge for a reach is represented by a given flow law (like Manning's equation). Therefore, the river discharge is expressed by a number of hydraulic variables such as the cross-section area of the channel, an empirical measure of flow resistance, the hydraulic radius, and the water surface slope of the river cross-section . After acquiring the measurements, the McFLI techniques search for the best fit for the cross-section area, flow resistance, and discharge through the time in a way that the mass and momentum are conserved. Apart from the river height, width, and slope, SWOT mission provides discharge by applying different McFLI (e.g., BAM and MetroMan;Hagemann et al., 2017) and data assimilation techniques (e.g., Hierarchical Variational Discharge Identification (HiVDI; Larnier et al., 2021) algorithm, and SIC4dVar (Oubanas et al., 2018).
Due to parameter equifinality in the flow-law inversion, the performance of SWOT discharge algorithms highly depend on the availability of prior values (Gleason & Durand, 2020;Hagemann et al., 2017). As Gleason and Durand (2020) state: the more we know about a river to start with, the better the final estimate of discharge will be. One viable choice for the prior information is either estimated or measured discharge data in the gauged and semigauged river basins. Several studies have presented approaches to develop empirical rating curves between simultaneous ground-based discharge observations and spaceborne measurements of a hydraulic quantity of the corresponding river reach like water level (e.g., Coe & Birkett, 2004;Getirana, 2012;Papa et al., 2010;Tarpanelli et al., 2019;Zakharova et al., 2006) or river width (e.g., Brakenridge et al., 2005;Elmi et al., 2015;Pavelsky, 2014;Smith et al., 1996;Tarpanelli et al., 2013). For developing river width-discharge model, the location of river centerline and cross-sectional river width can be automatically obtained from satellite images benefiting from the softwares like RivaMap (Isikdogan et al., 2017) and RivWidthCloud , which is the adapted version of RivWidth algorithm (Isikdogan et al., 2017;Yang et al., 2019) to the Google Earth Engine platform (Gorelick et al., 2017).
Since the number of active gauges has been decreasing in the recent years, obtaining prior information from gauges is no longer an option for most of the basins. However, the legacy discharge measurements of inactive gauge stations can still be used. Developing a reliable empirical relationship between the new remote sensing data (e.g., height, width) and legacy in situ discharge allows to extend the number of gauges with available prior discharge information (Tourian et al., 2013). Recently, Feng et al. (2019) compared the discharge estimates of the conventional rating curve and those from BAM technique, which performs using width observations solely. They concluded that BAM method can outperform the conventional rating curve as long as river discharge shows a stationary behavior. Gleason and Durand (2020), in their review paper, categorized the space-based rating curve technique as the calibration of local channel hydraulic relationships which is appropriate to apply in gauged, 10.1029/2021WR030277 3 of 27 semigauged and regionally ungauged basins. The primary goal of this group of algorithms is to provide discharge as-accurate-as-possible. Leopold and Maddock (1953) observed, for the first time, the strong dependency of cross-sectional river flow width, depth and velocity values, and the corresponding river discharge. Accordingly, they described three fundamental power-law functions in hydraulic geometry where Q is the river discharge, W is the river flow width, D is the river water mean depth, and V is the river water mean velocity. Since the discharge Q is the product of river width, depth, and velocity, a + c + k and b ⋅ f ⋅ m should equal to one. Generally, the corresponding coefficient and exponent in each power-law equation are empirically estimated from simultaneous measurements of discharge and river width, depth, or velocity. The measurements can be repeated at a single river cross-section, at-a-station hydraulic geometry (AhG), or for some fixed frequency of discharge between cross-sections either downstream or on other rivers, downstream hydraulic geometry (dhG) (Gleason, 2015). The AhG relations as functions of cross-section hydraulic (average velocity) and geometry (bankfull width and maximum depth) have been proven qualitatively (Ferguson, 1986) and quantitatively (Dingman, 2007).
Since the AhG exponents and coefficients describe the channel as a regular geometry (Figure 1a), modeling error may occur due to the lack of measurements in all percentiles. Since the number of measurements at moderate flow are usually dominant in the scatterplot (Figure 1a), AhG coefficients best represent the middle part of the point cloud. This may cause significant error in estimating discharge especially at high and low flow, at both ends of the scatterplot.
To compensate this error, Mersel et al. (2013) defined two linear trends in the width-height relationship, the socalled Slope-Break method which is further investigated by Schaperow et al. (2019). The same assumption can be made for defining the river width-discharge rating curve ( Figure 1b). Accordingly, the point cloud is divided into two (or more) parts and then the separate power-law models are defined for each part. By applying this technique, the mismodeling will be reduced. However, the performance is sensitive to the number and location of Slope-Break points. The river cross-section is modeled through (a) a power-law rating curve, (b) two rating curves using a Slope-Break point, (c) a nonparametric mapping function. Since our aim is to estimate river discharge using river width measurements, the corresponding relationship in Equation 1 will be redefined as Q = aW b .
In this study, our approach is to rely on a nonparametric mapping function. Generally, the nonparametric algorithm seeks the best fit for each point cloud after breaking up the scatterplot into several small parts to generate a mapping function. As a result, the obtained mapping function represents the river cross-section geometry more efficiently. We suggest to derive such a mapping function through matching quantiles of river width and discharge realizations, obtained through a Monte Carlo simulation. In fact, this can be seen as a multistage rating curve (Reitan & Petersen-Øverleir, 2009) in a way that the river width-discharge relationship can be varied in each percentile. We validate our method over 14 different river reaches along the Niger, Congo, Po, and Mississippi rivers.
To elucidate the method's performance, the results over four river reaches are provided in full detail. Our results show that the proposed algorithm can mitigate a possible mismodeling and compensate the effect of erroneous measurements. Moreover, the proposed algorithm delivers a meaningful uncertainty for the estimated discharge and allows us to calibrate in situ discharge uncertainty.

Data Sets and Case Studies
The proposed method is illustrated in full detail on four river reaches in Africa and Europe in Figure 2. For further assessment, the estimated discharge over 10 river reaches in Mississippi river basin ( Figure 5) will analyzed. Case Studies and Data for Implementation and Illustration The first river reach ( Figure 2a) is part of the Niger River at the confluence of Benue River near the city of Lokoja in the south-central of Nigeria. The second river reach ( Figure 2b) is also selected along the Niger River. The capital of Mali, Bamako, is located at the bank of this river reach upstream from Koulikoro. Part of the Congo River near the cities Kinshasa and Brazzaville is the third case study (Figure 2c) of this paper. The last case study is part of the Po River near Borgoforte in Italy (Figures 2d).
The main data sets used for developing river width-discharge models are ground-based river discharge measurements and average river reach widths, obtained from the MOD09Q1 product of modis imagery. In this section, only those parts of the time series that are used during the training and validation periods (specified in Table 1) are presented. River Discharge From Ground-Based Stations In situ river discharge measurements are available for each case study to develop the discharge estimation models during the training period and also for assessing the estimated river discharge during the validation period (Figure 3). River discharge observations are collected from two databases: Global Runoff Data Centre and Agenzia Interregionale Fiume Po, lit. "The Interregional Agency for the Po River", see Table 1.
The uncertainty of in situ river discharge measurements are usually not provided in the data sets. Since estimating the uncertainty of the gauge-based discharge measurements is not straightforward, we consider here a multiplicative error of 10% of the signal as the measurement uncertainty in the first iteration.

River Width From Satellite Imagery
In order to obtain time series of average river width variation for each river reach, we apply the region-based classification algorithm introduced by Elmi et al. (2016). In this method, a Markov Random Fields (MRF) model was developed to consider all types of available information in the image stack, including pixel intensity, and spatial and temporal interactions between pixels. Then the Maximum A Posteriori (MAP) solution for the Bayesian framework was found in order to determine the most probable water mask at each epoch. Since the high computational effort of finding a global solution was a serious concern, the problem was reshaped as an energy minimization, for which we developed an undirected two-terminal graph and defined the max-flow solution for the graph by augmenting all possible paths between the two terminals. The final residual graph offered the binary water mask representing the MAP solution. Moreover, the uncertainty of the river mask was defined by measuring the marginal probability of all nodes in the final residual graph. The assigned labels for pixels with marginal probabilities higher than 10% in both water and land masks were retained, which left a third region that contains pixels with marginal probabilities <10% in both masks. This is the uncertain region, for which a proper label cannot be defined based on the available information. Elmi (2019) considered the area of the third region as the uncertainty of the river mask.
The dynamic river masks can be converted to effective river width (W e ) as introduced by Smith et al. (1996), by dividing the river reach water area by the length of the reach. For simplicity, from now on, we use the term river width instead of effective river width. River width time series of all four case studies together with the estimated uncertainty are presented in Figure 4.

Case Studies and Data for Further Validation
We further validate the proposed method over 10 river reaches ( Figure 5) in the Mississippi river basin, for which the width data without uncertainty are generated through a rather more straightforward approach on Google Earth Engine. For selecting the case studies, we only consider river reaches with a natural behavior. As a results, the river reaches that their behavior is regulated by dams, channel locks, or dikes are excluded. We select 10 representative reaches with different types of river width-discharge relationship to ensure a thorough validation of our proposed methodology (Table 2). Such an assessment also assures the transferability to other locations for a wider community. For these reaches, river discharge is obtained from USGS gauge stations (Table 2).
Time series of river width for each river reach is obtained from the global surface water (GSW) data set. Benefiting from the entire archive of the Landsat 5, 7, and 8 imagery acquired, the European Commission's Joint Research Centre in the framework of the Copernicus Program developed a consistent monthly global surface water (GSW) data set over the years 1984-2020 (Pekel et al., 2016). Figure 6 presents the time series of simultaneous monthly river width and discharge for all selected river reaches in the Mississippi river basin. In each river reach, 70% of simultaneous measurements are selected for developing the river discharge estimation models and the 30% remaining are used for validating the estimated discharge. The training and validation periods are separated using a dashed line in each time series of Figure 6. The GSW data set does not provide any information about the uncertainty of its water masks; therefore, we assume a multiplicative width uncertainty equal to 10% of the signal for the width time series.

Nonparametric Regression for Estimating the River Discharge
The remote estimation of river discharge is typically done via developing an empirical relationship between coincident measurements of gauge-based discharge and space-based river water level or width during a training period. When the model is established, the discharge can be determined using just the space-based measurements and the model. Since this empirical model characterizes the morphology of the river cross-section, it should be a one-to-one (bijection) and increasing function. The measurements typically generate a wide point cloud in their scatterplot ( Figure 7b) because of the measurement uncertainty. As a result, they cannot be directly used as a mapping function. Fitting a parametric rating curve to the scatterplot is the most straightforward solution for this problem Nonparametric regression techniques become increasingly popular in modern statistics since, as Hazelton (2005) stated, they let the "data speak for themselves" in determining the form of the mapping function instead of forcing upon them strong assumptions about the number and relationship of the model parameters (Hazelton, 2005). However, computation complexity and a less understandable mapping model are the costs of achieving a more realistic regression function. The ideal river width-discharge mapping function has the following characteristics: 1. Since a bijective relationship between river discharge and width must be established, the derived mapping function should be monotonically increasing 2. Since the mapping function projects the physics of river cross-section morphology, the mapping function should be piecewise continuous 3. Since the measurements in both data sets are contaminated by uncertainty, the model should be robust against noisy measurements 4. Since the uncertainty of river discharge and width measurements is available, the methodology should provide information about the uncertainty of the discharge estimates In the literature, a number of algorithms like smoothing kernel, Gaussian process, local polynomial, spline smoothing, and regression trees are offered to develop the lookup table and mapping function (Härdle, 1990).  However, to satisfy the aforementioned conditions, several constraints and assumptions need to be imposed. Tourian et al. (2013) suggest an alternative solution for developing a mapping function for estimating the river discharge using satellite altimetry measurements.
Discharge gauges and spaceborne sensors measure different physical quantities of the water flow repeatedly which leads to different statistical distributions due to the different sampling rate of the sensors. However, if the river cross-section morphology is statistically well represented (using a scatterplot of simultaneous measurements), the measurements of one data set can be used to estimate the counterpart value and its distribution. In order to obviate the need for simultaneous measurements, Tourian et al. (2013) present a statistical approach based on the matching of quantile functions of in situ discharge and altimetry water level measurements. In other words, they characterize the river section morphology by mapping the quantile functions of the measurements rather than the value of coinciding measurements.
The quantile function, which is the inverse of the cumulative distribution function (CdF) of a time series, express the probability that a certain value is exceeded in the time series. Since the quantile functions of river discharge and width have a same x axis (cumulative probability), it is possible to connect their y axis directly. The obtained scatterplot can be used as a mapping function between river discharge and width as long as there is only one corresponding value available for each probability in the quantile functions. The resulting mapping function satisfies the first two aforementioned conditions for an ideal mapping function. Since the quantile function has a monotonically increasing behavior in its entire domain, it is a one-to-one function. Therefore, we can conclude that the derived quantile mapping function is also a monotonically increasing and one-to-one function. Generally, nonparametric methods are more vulnerable to overfitting. However, overfitting error is not a concern in this study since the nonparametric model is attained through the quantile function of the measurements. As the quantile function has a monotonically increasing behavior in its entire domain, the effect of noisy measurements is reduced significantly and the model would not bend toward any specific measurements. The performance of the mapping function will be affected by the noise in the measurements, since the uncertainty of measurements is not considered in determining the CDFs. Moreover, unlike parametric methods, describing the probability estimation of the unknown parameters in the nonparametric methods is not straightforward. To cope with this issue, the proposed algorithm in this study will approximate the uncertainty of the mapping function by using random sampling referred to as Monte Carlo methods. The uncertainty of the both measurement types (river discharge and width) in the training period is available, therefore random samples can be generated using the value and the standard deviation of any given measurement. The flowchart in Figure 8 describes the procedure of the proposed algorithm.
The algorithm performs the following steps for obtaining the stochastic quantile mapping function (Figure 8) 1. Generating a number of Monte Carlo samples for both time series according to their measurement uncertainties 2. Deriving the collection of quantile mapping functions by matching all possible permutations of river width and discharge quantile functions 3. Estimating the mean quantile mapping function together with the uncertainty for each percentile from the collection of the quantile mapping functions 4. Evaluating the performance of the derived model by comparing the estimated and measured discharge performing a 3σ test 5. Terminating the algorithm if the variation of RMSE from the previous step is less than the defined ϵ otherwise performs the next step 6. Updating the measurement uncertainties according to the result of the 3σ test and returning to the first step The algorithm cannot estimate discharge accurately at the early iterations because prior knowledge is unavailable. In the first iteration, a multiplicative uncertainty of 10% of the signal is considered as the discharge uncertainty. In each further iteration, the uncertainties of the input discharge are updated by considering the residuals between the measured and estimated discharge. The iteration for improving the input discharge uncertainty is terminated if the algorithm reaches a convergence in terms of variation of RMSE between the estimated and measured discharge. After termination, the river width-discharge mapping function that leads to the minimum possible discharge RMSE during the modeling period is obtained. Moreover, the standard deviation of the mapping function stack at each percentile is introduced as the uncertainty of the discharge estimate. The proposed algorithm is explained in detail in Appendix A.
In the remainder of the paper for simplicity rating curve refers to the parametric power-law relationship with one coefficient and exponent, deterministic quantile mapping function refers to the algorithm by Tourian et al. (2013) and the proposed algorithm in this study is called stochastic quantile mapping function.

Result and Validation Over the Niger, Congo, and Po River
The proposed methodology will be applied to the case studies introduced in Section 2 to evaluate the performance of the stochastic quantile mapping function for estimating the river discharge using the river width measurements. Apart from validating the estimated discharges, we will compare the performance of the proposed algorithm to the conventional rating curve technique, as well. For each river reach, the measurements during the training period (Table 1) are used for developing the models since the rating curve technique depends on the simultaneous observations. The parametric rating curve is modeled by a power-law model between river discharge and width measurements. In order to estimate the model parameters and the variance-covariance matrix we followed the Gauss-Helmert (Gh) formulation, described in Appendix B. The developed models for each case study are plotted over the scatterplot of the simultaneous measurements in Figure 9.
The performance of parametric models depends highly on the scatteredness of the point cloud. In the first and second case study (Figures 9a and 9b), the parametric rating curve (red line) and the nonparametric models (black and blue lines) similarly pass through the point cloud of measurements in the low discharge range. In the medium discharge range, parametric and nonparametric models deviate from each other because there are fewer data points to constrain them. Since the parametric model is fitted for the entire point cloud, the sparse measurements in the medium and high discharge values do not have a significant contribution in defining the model parameters.
On the other hands, the nonparametric model minimizes the measurement errors in each percentile, therefore the mapping function in each percentile only depends on the corresponding measurements.
Simultaneous measurements of the Congo River reach (Figure 9c) show a wide point scatter. This part of the river passes through urban areas and the river shoreline is restricted in Brazzaville city. Therefore, the river does not behave naturally and the correlation between the measurements decreased significantly. Since the surrounding area of the river reach is vast and flat, the river reach variations cannot reflect the river discharge variations which is measured at a gauge station located on the downstream of the river in Brazzaville city. These reasons cause a phase shift between the measurements which leads to the wide point cloud.
In the fourth river reach (Figure 9d), the width measurements does not represent the dynamics of the river properly due to a lack of measurements in the high percentiles and also a high noise level in the low percentiles. Since the parametric rating curve model passes through the point cloud in the medium discharge, it estimates discharge at the low and high end with a large error. In this case, the deterministic quantile mapping function approach can also not provide proper estimates because the width measurements are not well distributed through all the percentiles. As a result, the deterministic model (blue line) only interpolates the values for those percentiles without any width measurements. On the other hand, the stochastic quantile mapping function can overcome to this issue by matching the stack of various realization of the time series in each iteration.
A comparison between the deterministic (blue lines) and stochastic (gray line) quantile mapping functions in all case studies shows that the deterministic approach can only provide reliable width-discharge model if enough measurements in both data set, which are not corrupted by noise, are available in all percentiles (like Niger River reaches). In case of lacking a measurement in any percentile, the model estimates discharge by interpolating between the nearest available measurements (Figure 9d). Moreover, since measurements in the deterministic approach are considered error-free, the algorithm is very sensitive to the level of the noise in the measurements. As a result, the approach cannot alleviate the effect of noisy measurements in river width-discharge modeling. Among these models, only rating curve and stochastic quantile mapping function technique can provide the discharge uncertainty ( Figure 10).  The confidence envelopes for the parametric models are estimated by propagating the variance-covariance matrix of the model parameters (Equation B13). The confidence envelope describes how precisely the fitted line passed through the point cloud. Since a power law is fitted to the point cloud, the confidence envelope is typically small for the low values and get larger by increasing the values. As a result, the obtained model uncertainty cannot reflect the model performance for different parts of the distribution. On the other hand, the uncertainty for the stochastic quantile mapping function provides a more realistic evaluation of the precision of the estimated discharge, since each point in the mapping function has its own uncertainty describing the availability and accuracy of discharge and width measurements at each percentiles.
In Figure , the measured discharge is compared to the estimated discharge obtained by the rating curve and the stochastic mapping function techniques during the validation period. Since the correlation between river width and discharge measurements in the first two river reaches is high, the derived models estimate the river discharge accurately. Therefore, the residual values (Figures 11a and 11b) in most epochs are <0.05 km3/day. However, the models are not able to estimate discharge accurately in the wet season because of lack of width measurements. As expected, both parametric and nonparametric models present an acceptable performance in the first two case studies. Comparing the discharge uncertainties estimated by both models show that the parametric models only provide an artificial envelop around the estimated discharge. The stochastic information provided by the rating curve models is quite misleading, since the pink envelope is narrow at low discharge and gets wider with increased discharge. As a result, the discharge uncertainty in high discharge is extremely large despite the relatively small difference between estimated and measured discharge. On the other hand, the uncertainty estimates by the proposed algorithm offer a measure of the accuracy of the estimated discharge. Therefore, the uncertainty envelope only gets larger when the residual between estimated and measured discharge increases.
In the Congo River reach (Figure 11c) due to the phase shift in the training period, both models also estimate the discharge with a phase shift especially in the wet seasons which leads to a large residual. Nevertheless, in the dry Figure 11. Comparison between estimated discharge via the parametric rating curve models and nonparametric stochastic quantile mapping functions and measured river discharge during the validation period. season, the residuals are small, typically <0.5 km 3 /day. Because of the high level of the noise, the model uncertainty obtained by the parametric rating curve is too large in all percentiles.
The comparison between the time series for the Po River reach (Figure 11d) shows that the estimated discharge values cannot follow the measured discharge time series properly. The main reasons for this poor performance are the narrow river width with respect to the MOdIS pixel size and the lack of width measurements in the upper part of the quantile function. As a result, the derived model estimates the river discharge with a large residual in the wet season. Also, the large uncertainty of the width measurements causes a large uncertainty envelope for the estimated value. The size of the uncertainty envelope gets larger in the wet season since the model uncertainty is large on that part. The estimated discharge using the parametric model (red line) can only estimate discharge accurately at the low discharge since the model has been developed based on the simultaneous measurements which are concentrated in the low discharge.
To assess whether the estimated discharge values follow the same statistical distribution as the measured discharges, Q-Q plots of measured and estimated discharge are generated ( Figure 12). The estimated values generate exactly the same statistical distribution as measured values, if the points in the Q-Q plot align along the diagonal. The Q-Q plots show that both parametric rating curve and the stochastic quantile mapping function technique can reproduce the same distribution as the measured values in the early percentiles for the first three case studies. However, from the middle percentiles, the deviation from the identity line is noticeable especially in Congo and Po River reaches. Accordingly, in Congo River reach, both models tend to overestimate the peaks when the discharge is larger than 5 km3/day. The phase shift between the measurements prevents the model from matching the corresponding percentiles correctly. In the Po River reach, the estimated discharge via the rating curve model cannot reproduce the distribution of the discharge due to the large mismodeling. On the other hand, the discharge estimates obtained from the proposed algorithm follow the identity line in the low and medium discharge values. However, the discharge estimates deviate from the identity line in high percentile since the model underestimates the discharge in latest percentiles. The main reason for this large deviation is the absence of width measurements in the high peaks which leads to a mismatch between the statistical distribution of river width and discharge. As it is mentioned in Allen et al. (2020) and Nickles et al. (2019), the mismatch between the width and discharge distributions is the main reason for the large error at high discharge. In our examples, except the Po River reach, the statistical distributions of river width and discharge sparsely cover large values, so the mismatch between the quantile functions does not occur. We suggest analyzing the histogram of in situ discharge and width measurement distributions since the proposed algorithm is sensitive to the mismatches between the quantile distributions.
The metrics in Table 3 assess the performance of the rating curve and the stochastic quantile mapping functions quantitatively. The statistical indicators for the Niger River reaches show that both models can estimate discharge accurately. In both case studies, the normalized Root Mean Square Error   The performance of the methods significantly decreases over the Congo River reach. The phase shift between river discharge and width time series during the training period is the main reason for the relatively poor performance (normalized RMSE 20% ). Considering the very low NSE value (0.17) for the rating curve model, there is no significant difference if we replace the estimated discharge with the mean measured discharge of the training period. However, the nonparametric model performs better in this example, since the NSEs increased to 0.31. Although this value is not high for the NSE, but it is minimally acceptable with respect to the complex relationship between river width and discharge.
The statistics for Po River reach show the poor performance of the models due to the narrow width and unusual fluctuations. The normalized RMSE is too large (52%) for the rating curve model because of the huge mismodeling in the high peaks. Consequently, the NSE coefficient for the rating curve techniques is negative, since this parameter is very sensitive to the large residuals. The magnitude of the RMSE is decreased significantly by applying the proposed method, since the mismodeling is compensated in the estimated values. As a result, the NSE coefficients increased from −8.9 for the rating curve technique to 0.13 for the proposed method. Although the NSE value obtained by the proposed algorithm is not impressive but the improvement from a negative to positive value shows the contribution of the mismodeling. Finally, the probabilistic performance of the rating curve and stochastic quantile mapping function model is compared by the means of the region-wise scatterplot ( Figure 13) The probabilistic performance of estimated discharge via rating curve and stochastic quantile mapping function can be assessed through the region-wise scatterplot of the obtained residuals and uncertainties ( Figure 13). In such an assessment, the scatterplot of residual ( ) and estimated uncertainty (σ Q ) is separated into three different regions: (a) optimistic region (yellow), (b) realistic region (white), and (c) pessimistic region (red). Those estimates located in the optimistic region carry optimistic uncertainties leading to ∕ Q larger than 3. In the pessimistic region given an unnecessary large uncertainties, the aforementioned ratio will be smaller than 1/3. Within the white region, all estimated discharge passes the 3σ test successfully and also their uncertainty is not unnecessarily large. In an ideal situation, the majority of points should be located in the white region where the residual and uncertainty are both less than the standard deviation of signal that defines the gray boundary. In the gray region either the obtained residual or the estimated uncertainty is larger than the standard deviation of the signal. Accordingly, the dark red region is a region with pessimistic uncertainty larger than the standard deviation of signal.
Over Niger-Lokoja, the rating curve technique (Figure 13a) delivers a pessimistic uncertainty for 49% of measurements: 33% in the red region and 16% in the dark red region. The rating curve technique can only deliver realistic uncertainty for 50% of measurements, with 6% in the gray region. The distribution of the discharge estimates in Figure 13b shows that our proposed algorithm can estimate discharge with more realistic uncertainty as the number of measurements in the white region increases from 44% to 56%. Accordingly, the population of the measurements in the dark red and red is reduced to 7% and 30% from 16% and 33%, respectively. Over the Niger-Koulikoro, we also observe an increase in the number of measurements with realistic uncertainty. There, the number of estimated discharges in the white region is improved from 41% to 55%. These improvements along with very similar NSE and RMSE from both methods over the Niger-Lokoja and the Niger-Koulikoro (Table 3) highlight the fact that the advantage of our proposed method should not be sought in standard performance metrics.
In the Congo River reach, the rating curve technique fails to locate measurements inside the white region since the estimated uncertainty is extremely large. As a result, 68% of the measurements are in the dark red region with uncertainty estimates that are three times larger than the residual and also larger than the standard deviation of the signal. The quantile mapping function drastically improves the performance since the population of the measurements in the dark red region is reduced to 6%. Although still 13% of the measurements lie in the pessimistic region, the majority of the points are moved to the white (37%) and gray (42%) regions. In the Po River reach, the performance of the rating curve technique is not reliable since none of discharge estimates is located in the white region. Due to the obtained large uncertainty for the discharge estimates, most of the points locate in the upper part specifically in the dark red region of the scatterplot (Figure 13g). In this example, the stochastic quantile mapping function improves the uncertainty estimate by reducing the population of points in the dark region from 62% to 29% and having 18% of points in the white region.
Comparing the distribution of points in the scatterplots and the magnitude of the discharge values (represented by the color of the points) shows a direct relationship between the magnitude of the discharge and the uncertainty estimates via the rating curve model. In all four cases, the estimated uncertainty is relatively small in low discharge, and as the magnitude of the discharge increases, the discharge uncertainty determined by the model also increases. Since in parametric models the discharge uncertainty is a function of the model parameters and the input data, such dependence is unavoidable. However, this undesirable pattern is not found in the scatter plots of the proposed algorithm. Therefore, it can be concluded that the uncertainty of the estimated discharge by the stochastic quantile mapping function is independent of discharge magnitude.

Result and Validation Over the Mississippi River Basin
With the promising results in the previous sections, the question remains if the proposed method is applicable to other river reaches and how its performance depends on the availability of proper width uncertainty information.
To this end, we assess the performance of the method on 10 river reaches in the Mississippi river basin (MRB) ( Figure 5). The training and validation data set have been described in Section 2 (Table 2 and Figure 6).
The developed rating curve and stochastic quantile mapping models for each river reach are plotted over the scatterplot of the simultaneous measurements in Figure 14. It should be mentioned that since the GSW data set provides monthly water masks, we simulate monthly discharge measurements by accumulating daily discharge values. Comparing the rating curve models and stochastic quantile mapping functions shows that only in two river reaches (MRB-a, h), parametric and nonparametric models almost behave similarly and pass through the entire point cloud. The models for MRB-c reach present an identical behavior over the whole range of the scatterplot expect where the discharge is larger than 700 m 3 /s. In some river reaches (Figures 14d, 14e, 14i, and 14j), two models are almost identical in the low discharge but show a large deviation in high discharge values. The rating curve and stochastic quantile mapping technique present completely different river width-discharge models in Figures 14b, 14f, and 14g due to the presence of a wide point cloud in low discharge.
In Figure 15, the measured discharge is compared to the estimated discharge obtained by the rating curve and the stochastic quantile mapping technique during the validation period. The length of the validation period is different in each case study according to the availability of in situ discharge. In the majority of the river reaches,   To assess whether the discharge estimates via both techniques follow the same statistical distribution as the measured discharge, the Q-Q plots of measured and estimated discharge are generated ( Figure 16). Apart from MRB-c river reach (Figure 16c), the stochastic quantile mapping function can almost reproduce the same distribution as the measured discharge. However, the similarity decreases in the high discharge values in some examples (Figure 16d, 16f, and 16g) mainly due to the lack of the training data in high discharge percentiles. On the other hand, rating curve technique fails to generate the same distribution like the measured discharge. In the majority of examples, the rating curve models tend to underestimate discharge as in most percentiles the red points are located below the diagonal. The better results of stochastic quantile mapping functions are confirmed by NSE and RMSE metrics listed in Table 4.

Application to SWOT Discharge Estimation
SWOT discharge will be provided in two separate forms: gauge unconstrained and gauge-constrained, for which the in situ discharge data will be assimilated to SWOT discharge estimates (Stuurman & Pottier, 2020). The gauge-constrained discharge data will be the SWOT discharge product for a broader community and end users. However, one of the challenges in assimilating the gauge data with SWOT's estimated data is to have an appropriate uncertainty for the gauge data. The proposed method in this study provides realistic stochastic information for the in-situ data, as well as a realistic uncertainty for discharge estimates. Figure 17 shows the obtained realistic uncertainty from the proposed methodology for in situ discharge values (green line). The error envelop for each model can be considered as a confidence measure for the estimated discharge uncertainty function. The gray lines represent the commonly assumed multiplicative uncertainty (10% of discharge) within the community   for in-situ discharge, which is our starting uncertainty value within our iterative procedure. We note that with the intercept added to the green line, which can be interpreted as a constant noise level, the uncertainty of the in situ discharge is no longer a pure multiplicative uncertainty. In Niger-Lokoja, we suggest a lower uncertainty than 10% of discharge for high discharge values, while a high uncertainty of about 50% seems to be a realistic uncertainty for low discharge values. Similarly, for the Congo in situ discharge measurements, our method proposes a drastic reduction in uncertainty for high discharge values. While for the Niger-Koulikoro, the change is mainly an additional constant noise level, for the Po River case, no constant noise level is introduced and instead a lower ratio is proposed. The error envelops are relatively narrow in the Niger River reach examples since the simultaneous width and discharge measurements distributed in all percentiles and the modeling error did not occur. In the Congo example, the error envelop is extremely large mainly due to the phase shift between the width and discharge measurements in the training period and also large uncertainty of the width measurements. In the Po River reach, the wide error envelop indicates a low confidence to the estimated uncertainty function in the high discharge. Noisy width data and also lack of width measurements in the high discharge can be the reason for the large error envelop. In a few examples of MRB river reaches (Figures 17f, 17i, and 17j) the constant noise level is significantly large as the river width is highly uncertain when the river gets too narrow. In examples like Figures 17e, 17g, 17k, and 17n, the green and gray lines have almost the same slope but the green lines show a constant noise level.
In order to implement the proposed method, only time series of width measurements together with in situ discharge are required. To this end, within the first year of the SWOT mission, using the collected time-varying width, and surface water height supplemented by available data from satellite imagery, a suitable uncertainty level for the in situ river discharge can be obtained. A realistic in situ discharge uncertainty guarantees properly assimilating gauge data with SWOT discharge estimates to obtain a better gauge-constrained product.

Summary, Conclusions, and Outlook
Developing a single-stage rating curve between the ground-based discharge and spaceborne river water level or width is the most straightforward approach to obtain river discharge using space-based and ground-based measurements. However, this does not always lead to promising results, since the AhG power-law equations describe a river section with a regular geometry. This assumption may cause a large modeling error. Moreover, rating curves do not deliver a proper estimation of discharge uncertainty as a result of (a) the mismodeling and (b) the coarse assumptions made for the uncertainty of inputs.
In order to address the aforementioned limitations, we have introduced a nonparametric approach to estimate river discharge from spaceborne river width measurements. The algorithm proposes a stochastic quantile mapping function scheme by performing the following procedures, iteratively: (a) generating realizations of river discharge and width time series using the Monte Carlo simulation, (b) obtaining a collection of mapping function by matching all possible permutations of simulated river discharge and width quantile functions, (c) adjusting the measurement uncertainties according to the scatter of point cloud. The algorithm's estimates are improved in each iteration by updating the measurement uncertainties according to the difference between the measured and estimated values. To do so, first, we examine whether the estimated values are in the 3σ confidence interval of the measured values. The algorithm iterates by updating the measurement uncertainties until it converges. In fact, by performing the above procedure, instead of a single quantile mapping function, a stack of quantile mapping functions is generated by propagating the river discharge and width measurements according to their corresponding uncertainties. The magnitude of the mapping function collection explains how the model can estimate the river discharge at any given percentile.
We have employed the proposed method over 14 river reaches with different dynamic behavior along Niger, Congo, Po, and different rivers in the Mississippi river basin. Evaluating the estimates shows that the performance of the proposed method is superior to the single-stage rating curve technique especially in challenging cases. The rating curve technique can only present satisfactory results in river reaches (e.g., Niger and MRB-a river reaches) where the simultaneous measurements are distributed regularly over all percentiles. The performance of the rating curve technique decreases if the measurements are contaminated by noise; however, the stochastic quantile mapping function can mitigate the effect of noisy measurements. As a result, the NSE values in Congo and MRB-b river reaches improved from 0.17 and −0.67 to 0.31 and 0.45. Describing the channel as one with a regular geometry is a strong assumption that makes the rating curve technique vulnerable to variation in the river width-discharge relationship at different discharge percentiles. The unsatisfactory performance of the rating curve technique in cases like Po and MRB-e, f can be attributed to mismodeling due to irregular geometry. Dealing with the mismodeling, the stochastic quantile mapping function can perform significantly better than the rating curve for estimating discharge in challenging examples as the NSE value has improved from from −9 to 0.13 for Po River reach and from −1.19 and −0.5 to 0.63 and 0.37 for the MRB-e, f river reaches. The results over the Mississippi river basin with a multiplicative uncertainty model for the river width data show that the performance of the proposed algorithm is less influenced by the uncertainty of width data set and more with the data availability.
The stochastic quantile mapping function is a nonparametric data-driven method that does not need to characterize a river section geometry through a limited number of model parameters. This, in fact, is the main reason for the obtained better performance. However, the performance of the technique highly depends on the availability of enough measurements in both data sets. It means that if a certain percentile is not represented in the quantile functions of both variables, the model will rely on interpolation for estimating discharge on that percentile.
Moreover, the proposed algorithm allows us to obtain a realistic measure for discharge uncertainty. Estimating the uncertainty for the conventional rating curve model is straightforward since the measurement uncertainties are introduced to the adjustment model as elements of a weight function. However, the obtained covariance matrix for the model parameters is influenced by blunders and mismodeling. Moreover, evaluating the model uncertainty with a limited number of parameters leads to a coarse understanding about the accuracy of the model over its domain. This is highlighted in the region-wise scatterplots shown in Figure 13, where rating curve model fails to locate the measurements in the white region. On the other hand, within the proposed stochastic nonparametric method, at any given point of the mapping function, the uncertainty of the model depends on the distortion of the quantile function realizations in the corresponding percentile. Therefore, the model uncertainty in any percentile will explain the availability and accuracy of river discharge and width measurements in that specific percentile.
Estimating the uncertainty of gauge-based discharge measurements is not straightforward since a lot of parameters like change in channel dimensions, river bank erosion and sensor miscalibration have a contribution in the error budget. Therefore, according to Harmel et al. (2006), the uncertainty of the in situ discharge is conventionally ignored or replaced by an arbitrary margin of safety (e.g., a 10% multiplicative uncertainty). In this study, as a side product of our proposed algorithm, we showed that the optimal uncertainty for discharge does not match with the assumption of a multiplicative uncertainty. In fact, following our algorithm, the input uncertainty of discharge can be recalibrated. Our results show that the low discharge values would require a larger uncertainty than 10% for a better estimation of Monte Carlo realizations. Moreover, since SWOT discharge will be provided in two separate forms of gauge unconstrained and gauge-constrained, an improved in situ discharge uncertainty leads to an appropriate assimilation of gauge data with SWOT discharge to obtain a better gauge-constrained product.
In fact, both discharge and spaceborne river width are contaminated by certain noise spectrum that are time dependent. However, in this study, we have considered a simplified time invariant noise following a Gaussian distribution. Therefore, the results can be further improved by incorporating the time dependent uncertainties with non-Gaussian distribution within the simulation and using, which will be the focus of future work.

Appendix A: Implementation of the Stochastic Quantile Mapping Function Algorithm
The algorithm starts with generating 100 Monte Carlo samples for both time series, leading to 10,000 mapping functions, when matching all possible permutations ( Figure A1e).
The collection of quantile mapping functions in Figure A1e characterizes the uncertainty of measurements in both axes. Since we consider a multiplicative uncertainty (10% of the signal) for discharge measurements, the higher the discharge, the larger the uncertainty. On the other hand, the uncertainty of river width measurements decreases by increasing river width. As a result, the set of mapping functions is wide along the x axis at the left side of the plot because of the large uncertainty in the river width measurements. At the right side of the plot, the mapping function collection gets larger along the y axis with increasing discharge value. To evaluate the performance of the derived quantile mapping functions, the river discharge from width measurements during the training period has been compared to the measured discharge. Figure A2a presents the comparison between the estimated and measured discharge for the training period. For each river width measurement in the training period, river discharge is estimated several times using all quantile mapping function realizations ( Figure A1e). Each time series of estimated discharge is presented in gray and the estimated discharge via the mean mapping function is presented in black. Considering multiplicative uncertainty for the discharge measurements may lead to a huge difference in the estimated discharge values in the high peaks. Therefore, the difference between estimated and measured discharge gets larger with increasing discharge in Figure A2b. The statistical distribution of the residuals ( Figure A2c) almost follows a normal distribution with zero mean. The quantile-quantile plot (Q-Q plot) of the measured and estimated discharge is generated by plotting the quantiles of measured discharge against the quantile of estimated discharge ( Figure A2d). Since most of the points are on the identity line (blue line), both estimated and measured discharge represent a same statistical distribution. As a result, we can consider the mean quantile mapping function as an unbiased estimator for predicting discharge using the river width measurements. The influence of considering multiplicative uncertainty for the discharge measurements is also obvious in the Q-Q plot, as the vertical gray bars get larger with increasing discharge. To reduce this effect, we will modify the discharge uncertainties according to the obtained residuals. In order to assess the performance of the developed model, the estimated discharge will be compared to the measured discharge. To this end, we apply a 3σ rejection criterion to the residuals. This test reveals whether or not the mean of the estimated discharges is inside the confidence interval (±3σ) of the measured observations.
The result of the 3σ test for all the estimated discharges during the training period is plotted in Figure A3a. In this scatterplot, the accepted points are mostly located around the identity line, since these points have a small residual. Most of the rejected points are gathered in the bottom left side of the scatterplot. Although the magnitude of residuals is small in this part of the scatterplot, most of them were rejected due to the insignificant uncertainties. This pattern is more obvious in Figure A3b in which the residual between measured and estimated discharges are plotted against the defined uncertainty for the measured discharges. Points with a residual smaller than 0.1 km 3 / day are separated into two parts. At the origin of the scatterplot, the estimations are accepted since the residuals are very tiny or zero. By moving along the y axis, the population of gray points have increased because the confidence interval gets larger. On the other hand, the majority of the estimated discharges are located outside the accepted intervals by moving along the x axis. This analysis shows that considering a multiplicative uncertainty for the discharge measurements is not an appropriate choice. To overcome this issue, we develop a new function for the discharge uncertainty by assuming that the obtained residuals are the critical value for accepting in the 3σ test (the process, adjustment of the measurement uncertainties, in Figure 8). The same analysis has been done for the estimated and measured river widths.
The scatterplot of the measured and estimated river width values ( Figure A4a) shows a good agreement. In the bottom left side of the plot, the point cloud is wide around the identity line, because the width measurements are less accurate when the river is narrower. The point cloud in the scatterplot of residual between estimated and measured river width and the uncertainty of the width measurements ( Figure A4b) does not present any correlation between the residuals and measurement uncertainties. For most of the points, the residual is <100 m since the measurement uncertainty is between 30 and 180 m. Therefore, all the estimated river width values are accepted in 3σ test. Therefore, updating the river width uncertainty is not required in Figure A2. (a) is the time series of estimated and measured river discharge. The difference between measured and estimated discharge is plotted in (b). (c) is the histogram of residuals. The Q-Q plot of the measured and estimated discharge is plotted in (d). The unit of all values is km 3 /day. the next iteration. For further evaluation, we repeat the test by reducing the significance level to 2σ and 1σ. When the confidence interval is reduced to 2σ, just five points are rejected by the test. The number of rejected points increases to 26 when the confidence interval is reduced to 1σ. This analysis shows that the defined uncertainty for the river width measurements is more realistic than those assigned for the river discharge measurements.
After updating the discharge uncertainties, the algorithm returns to the first step and repeats the procedure. The iteration will be terminated, if the RMSE of the residuals changes by less than ϵ, which means no further improvement in discharge and width estimation.
The procedure of deriving quantile mapping function is repeated after modifying the uncertainty of measurements. Figure A5 shows that after the first iteration, the magnitude of ΔRMSE is at the level of 10 −5 . After two more iterations, the magnitude of ΔRMSE reaches a relatively constant level around the value 10 −7 , meaning that no significant improvement can be obtained after the iteration 4 ( Figure A5). Therefore, it can be concluded that the algorithm reaches a robustness level after adjusting the uncertainty of the measurements.
The comparison between two sets of mapping functions ( Figure A6) shows a significant difference. The collection of quantile mapping functions in Figure A6b characterizes different patterns in both ends of the plot because of different stochastics for discharge measurements. In the last iteration, more possibility is given to the river discharge realizations to vary especially in low discharge; therefore, the set of quantile mapping functions grows faster when the discharge is lower than 0.1 km 3 /day. However, the collection eventually grows slower than the first iteration at medium and high discharge. The main difference between the two plots is the reduction of the quantile mapping function set in the right tail of the curves. For analyzing the performance of the new mapping function, adjusted observations will be compared to the measured discharge again.
The modified quantile mapping functions can predict the discharge more accurately, as the magnitude of the residuals are reduced significantly (Figures A7a and A7b). Using the mean quantile mapping function, the maximum residual has been reduced from 0.4 to <0.1 km 3 /day and the RMSE of the residuals reduced from 0.06 to 0.05 km 3 /day. The histogram of the residuals ( Figure A7c) presents almost a zero mean normal Figure A3. The scatterplot of estimated and measured discharge is in (a). In (b), the defined measurement uncertainties are plotted versus the residuals. In both plots, the residual of the red points exceeds three times its defined standard deviation. The unit of all values is km 3 /day. Figure A4. The scatterplot of estimated and measured river width in (a). The measurement uncertainties are plotted versus the residuals in (b). The unit of all values is km 3 /day. distribution (the mean of the residual is −0.006 km 3 /day). Since the magnitude of the residuals has been reduced in the final iteration, the histogram presents a normal distribution with a smaller standard deviation (−0.059 km 3 /day) from the first iteration (−0.051 km 3 /day). Finally to evaluate whether the mean of estimated discharges is inside the confidence interval of the measured discharge, the 3σ hypothesis test will be performed again.
The result of the 3σ test presents a significant improvement in estimating the discharge after the iterations as the number of rejected points (red points) is reduced from 111 to 45 points. Most of the rejected points are located in the left side of the point cloud, where the discharge value is small. Although the residuals of the rejected points are small, they are not located inside the confidence interval of the measured discharge since the assigned uncertainties are small, too. The scatterplot in Figure A8b explains the reason why the number of accepted points is increased in the final iteration. Since a constant value is used for the measurement uncertainties, the point cloud does not accumulate in the origin of the plot. Therefore, most of the points with a small residual are accepted in the 3σ test. The comparison between the discharge uncertainty in the first and last iteration ( Figure A8c) shows that considering a 10% multiplicative uncertainty is not representative of the real uncertainty. In the proposed algorithm, the function of discharge uncertainty is modified at each iteration according to the difference between estimated and measured discharge. As a result, the obtained function in the last iteration represents a realistic estimate of the measured discharge uncertainty. Therefore, we can claim that our method recalibrates the discharge uncertainty based on the obtained residual. Figure A5. The evolution of RMSE through iteration process.

Appendix B: Gauss-Helmert Adjustment Model
It is assumed that a power-law relationship is established between river discharge and width in a natural river section where Q is in situ discharge measurement, W is the average width of the river reach, and a and b are the unknown model parameters. Once the model parameters are defined, we are able to predict river discharge by using the river width measurements and estimated model parameters.
The observations (river discharge and width) in both sides of Equation B1 are corrupted by noise. Therefore, to estimate the model parameters (a and b), the squared sum of inconsistencies must be minimized. One approach is to reformulate the above equation as a nonlinear conditional equation with unknowns. This model is called the general model of adjustment or Gauss-Helmert model are equal to zero. Therefore, the initial values for the Tylor series (a 0 , b 0 ) can be obtained by solving the equation using the ordinary least squares. Now, we write the Taylor series expansion Figure A8. The scatterplot of estimated and measured discharge is in (a). In (b), the defined measurement uncertainties are plotted versus the residuals. In both plots, the residual of the red points exceeds three times its defined standard deviation. The comparison between initial and modified discharge uncertainty function is presented in (c). The unit of all values is km 3 /day. ( , , , ) = ( 0, 0, In the end we can write the equation like where ω is the vector of misclosures. e is the vector of observation inconsistencies, A is the design matrix, and δξ is the vector of unknowns. Our aim is to determine the optimal values for the vectors e and δξ by minimizing the target function e T Pe, where P is the weight matrix which includes the uncertainty of measurements By replacing the equivalent of the vector ̂ from the Equation B7 into the Equation B9, we write the solution of this linear system By estimating ̂ and ̂ , and replacing them into Equation B7, the vector of observation inconsistencies (̂ ) will be updated too. After each iteration, the initial values of the unknowns and the observations will be updated. The iterations will be continued until the convergence is achieved for the norm of ̂ and ̂ . Apart from the model parameters, we can estimate the variance-covariance matrix of the unknown parameterŝ and by estimating the posterior variance of unit weight ( ) where n is the number of observations, P is the weight matrix, and ̂ contains the final estimated inconsistencies of the observations. Therefore, the variance-covariance matrix of adjusted unknown parameters will be defined aŝ̂=̂̂.
Calculating the variance-covariance matrix of the adjusted unknown parameters helps us to provide an uncertainty measure for the estimated river discharge through the error propagation of Equation B1 2 = 2 2 + ln( ) 2 + ( ln( )) 2 2 2 + ( −1 ) 2 2 , in this equation, a and b are the model parameters. σ a , σ b , and σ ab are the elements of matrix ̂̂ and σ W is the uncertainty of the measured river width (W). The goodness of the line regression is usually measured by assessing the discrepancy between the observed and estimated values by the model (the vector ̂ ).

Data Availability Statement
The data set containing the river width and estimated discharge time series of the 14 case studies are be publicly available online at hydrosat.gis.uni-stuttgart.de.