# Combining Landsat 5 TM and UAV images to estimate river discharge with minimal groundbased flow velocity and water level observations

Maomao Li, Changsen Zhao, Qi Huang, Tianli Pan, Hervé Yesou, Françoise

Nerry, Zhao-Liang Li

<figure><img src="image_1.png" /></figure>

## To cite this version:

Maomao Li, Changsen Zhao, Qi Huang, Tianli Pan, Hervé Yesou, et al.. Combining Landsat 5 TM and UAV images to estimate river discharge with minimal groundbased flow velocity and water level observations. Remote Sensing of Environment, 2025, 318, pp.114610. 10.1016/j.rse.2025.114610. hal-05157714

HAL Id: hal-05157714

https://hal.science/hal-05157714v1

Submitted on 11 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

# Combining Landsat 5 TM and UAV images to estimate river discharge with minimal ground-based flow velocity and water level observations

2

3

Maomao Li<sup>a,♀</sup>, Changsen Zhao<sup>a,b,c,♀</sup>, Qi Huang<sup>d, e, f *</sup>, Tianli Pan<sup>a</sup>, Hervé Yesou<sup>b</sup>, Françoise Nerry<sup>b</sup>, Zhao-Liang Li<sup>b</sup>

4

<sup>a</sup> College of Water Sciences, Beijing Normal University, Beijing 100875, PR China

5

<sup>b</sup> ICube, UMR 7357, University of Strasbourg, UNISTRA, 300 Bld Sebastien Brant, CS 10413, 67412

6

Illkirch, France

7

<sup>c</sup> School of Environment & Sustainability, University of Saskatchewan, Saskatoon SK S7N 5C9 Canada

8

<sup>d</sup> Key Laboratory of Groundwater Conservation of Ministry of Water Resources, China University of Geosciences, Beijing 100083, China

9

10

<sup>e</sup> School of Water Resources and Environment, China University of Geosciences, Beijing 100083, China

11

<sup>f</sup> UK Centre for Polar Observation & Modelling, Centre for Excellence in Environmental Data Science, Lancaster University, Lancaster, LA1 4YW, UK.

12

13

♀: Contributing equally to this study and therefore first authorship shared;

14

*: Corresponding to: 2024010081@cugb.edu.cn (Qi Huang).

15

16

**Abstract:** River discharge plays an indispensable role in maintaining the stability of the hydrosphere system and eco-environment. Previous methods that utilize satellite imagery to estimate discharge over poorly gauged basins are generally tailored for large rivers and heavily reliant on ground-based measurements. Consequently, uncertainties often escalate when these methods are applied to medium-sized rivers. Based on Landsat 5 Thematic Mapper (TM) and unmanned aerial vehicle (UAV) images, this study proposed a framework for estimating the discharge of large and medium rivers in poorly gauged areas with minimal requirements on ground observations. It comprises (1) a modified C/M method, which considers the spatial heterogeneity of rivers using single-site observation data, and (2) a newly developed method for estimating river bathymetry with zero discharge measurements (RIBA-zero). Results show that, utilizing the modified C/M method, rivers wider than three times the satellite resolution (i.e., 90 m) exhibit a relative root mean square error (rRMSE) of 0.23 in the velocity estimation.

17

18

19

20

21

22

23

24

25

26
Narrower rivers display a slight increase in the rRMSE (0.41), which is still within an encouraging range.

For both types of river widths, the accuracy of flow velocity estimation is higher during high-flow periods compared with the low-flow counterparts. In terms of the flow area estimation, the RIBA-zero method is much more suited for parabola-shaped cross-sections (rRMSE=0.22) and flood seasons (rRMSE=0.35).

Additionally, when replacing 30-m Landsat 5 TM with 10 m-resolution Sentinel-2 imageries, the approaches make a significant improvement in velocity estimation for rivers narrower than 90 m across all periods, exhibiting great potential to estimate discharge in medium rivers with finer resolution satellite imageries. The framework requires a few ground observations for discharge estimates with the Nash–Sutcliffe efficiency coefficient (NSE) reaching ~0.9, thereby greatly facilitating hydrology-related studies in poorly gauged areas with profound implications for sustainable water resources management worldwide.

**Keywords:** poorly gauged rivers; river discharge; bathymetry; Landsat 5 TM; UAV

# 1. Introduction

River discharge plays an indispensable role in maintaining the stability of the hydrosphere system and eco-environment. However, traditional methods for measuring river discharge are often laborious, time-consuming, and face challenges related to rugged terrain, especially in poorly gauged areas. This can limit their widespread application, impacting environmental disaster response and emergency monitoring. Therefore, a more efficient approach is needed to derive accurate river discharge in near real-time for poorly gauged basins (Gleason and Durand 2020; Harlan et al. 2021; Ishitsuka et al. 2021; Muhebwa et al. 2024).

Multispectral satellite images with high spectral and spatial resolutions and wide coverage can be used for both large and medium rivers (Feng et al., 2019; Li et al., 2019; Wu et al., 2024). There are primarily two approaches to estimating river discharge using remote sensing techniques. The first one integrates remote sensing with hydrological models that are driven and/or calibrated by remotely sensed information (Huang et al., 2020; Sun et al. 2018). Although hydrological models are suitable for estimating long-term continuous river discharge, there still exist high uncertainties in model inputs and parameter calibration, and parameters are not easily transferable (Mersel et al., 2013; Huang et al.,
2022). The second approach is to establish empirical regression equations between remote sensing-observed river characteristic variables (such as width and water level) and discharge. It estimates river discharge through empirical regression and hydraulic formulas based on remotely sensed observations.

Depending on the number of hydraulic characteristics involved, this approach can be categorized into single-variable and multi-variable empirical methods. Single-variable methods focus on a single hydraulic variable (e.g., water level or width) to estimate discharge, inclusive of the water level-discharge rating curve method (estimating discharge based on the empirical relationship between water level and flow rate) (Papa et al. 2010; Papa et al. 2012; Kim et al. 2019), river width (water fraction)-discharge rating curve method (using river width/ water fraction as a proxy for flow rate) (Brombacher et al. 2020; Smith and Pavelsky 2008; Pavelsky 2014; Lin et al., 2024), and the calibration/measurement (C/M) signal method (calculating river discharge via the ratio of reflectance values during floods) (Tarpanelli et al., 2013; Filippucci et al., 2022). Multivariate methods incorporate multiple hydraulic variables (e.g., water level, width, and velocity) to more comprehensively capture the combined effects on river flow.

These algorithms rely on empirical relationships derived from multiple factors, thereby providing a more nuanced and accurate estimate of discharge (Brinkerhoff et al., 2020; Durand et al. 2016; Durand et al., 2023; Durand et al. 2024; Feng et al., 2019; Frasson et al., 2023; Hagemann et al., 2017; Huang et al., 2018; Lin et al., 2023; Zhao et al., 2019).

In the second approach mentioned above, the C/M signal method is straightforward and widely applied, relying on the spectral differences between water bodies and land objects in near-infrared (NIR) images (Lu et al., 2010; Li et al., 2019). Its basis is that in any given area, the increases in water level would cause a decrease in the NIR reflectance value. As a result, the ratio (C/M) between the surface reflectance of a dry pixel, named C (calibration), and that of a wet pixel, named M (measurement), changes accordingly. That is, due to the variations in water volume during flood events, the reflectance of the wet measurement pixels (M) typically decreases, while the reflectance of the dry calibration pixels (C) remains constant. Therefore, during flooding, the C/M ratio is highly sensitive to the increase in the number of wet pixels, which in turn reflects changes in river discharge. Moreover, because discharge and velocity change synchronously, the C/M ratio also correlates closely with flow velocity (Tarpanelli et al., 2013; Tarpanelli et al., 2015; Tarpanelli et al., 2019; Tarpanelli et al., 2020). Tarpanelli et al. (2013)
have effectively established the relationship between C/M-V based on hydrological data from several stations on the same river. Currently, this approach allows for large-scale, high-precision, and low-cost river discharge estimation. However, the spectral behavior (C/M) and river morphology at different locations often exhibit significant spatial heterogeneity that is closely related to river flow velocities. For example, flow velocity affects the morphology, suspended substance concentration, and transparency in the river, leading to differences in the reflectance of optical imagery (Kwon et al., 2023). Additionally, different river morphologies lead to distinct hydrodynamic characteristics that in turn influence the distribution of flow velocity. Therefore, establishing accurate velocity models typically requires calibration data from multiple hydrological stations. However, hydrological stations are sparsely distributed, and obtaining confidential hydrological data from multiple stations poses substantial challenges in many regions across the world. It is highly necessitated to develop an approach that considers the spatial heterogeneity of both spectral behavior and river morphology at a large scale, relying on minimal hydrological stations—ideally one station—to compute flow velocity for any segment of a river.

In addition, when the river flow velocity has been obtained, the accurate flow area is essential for the river discharge estimation. Satellites with stereo imaging modes are eligible for this, as they can acquire large-scale topography and river geometry. There have been explorations into estimating river depths based on remote sensing data, utilizing water surface width and elevation information acquired from satellites such as SRTM and ICESat-2 (Domeneghetti et al., 2016; Coppo Frias et al., 2023). River depth is estimated based on the relationship between water level and river width, while high-accuracy observations of river cross-sections are still challenging (Mersel et al., 2013; Schaperow et al., 2019). When the water-surface width is less than 10 m, a measurement error greater than 0.5 m in the horizontal or vertical direction will finally result in exponentially increased errors in the flow area and discharge estimation. Evidently, it is difficult for satellites, including commercial satellites, such as QuickBird-2 or Pleiades, with a superior spatial resolution (better than 1 m ground sampling distance resolution) and higher revisiting frequency (1-3 days; twice a day for the Pleiades Neo), to directly monitor underwater conditions (Pan et al., 2016) and generate the bathymetry of the waterbed.

Instead, low-altitude unmanned aerial vehicles (UAVs) demonstrate the advantages of high flexibility, convenience, and high resolution. They can be used to monitor the water surface slope, cross-sectional
profiles above the water surface, and other hydraulic variables, with a centimeter-level resolution, and data acquisition can be performed at flexible hours. Therefore, UAVs have become among the major low-altitude remote-sensing platforms available today (Colomina and Molina, 2014; Che et al., 2020; Jiang et al., 2022) and have been widely used in hydrological applications (Kim et al., 2019; Erena et al., 2019). Despite this, ground-measured data or mathematical modeling are still among the most efficient ways to derive underwater bathymetry (Leon et al., 2006; Bonnema et al., 2016; Legleiter, 2015; Moramarco et al., 2019). However, in a specific poorly gauged region, where ground data are limited, it is difficult to obtain detailed ground-measured data necessary for calculating river discharge, which constitutes a major barrier to river discharge estimation. Therefore, there is an urgent need to mathematically model river cross-section without any ground-based bathymetry measurements.

The objective of this study is to present a new framework for the estimation of river discharge in poorly gauged rivers based on satellite and UAV data. We propose a method to reconstruct river bathymetry to calculate flow area (A) without any in situ bathymetry data and a method to calculate flow velocity (V) that considers the spatial heterogeneity of river spectral information based on the C/M-V relationship derived from only one gauging station. We applied the developed approaches across different reaches along the river and the accuracy of the velocity, bathymetry, and discharge estimates was validated using the gauged counterparts. These methods can provide reliable discharge estimates of poorly gauged rivers with profound implications for global river monitoring.

## 2. Study Area

This study selected the largest tributary in the middle reach of the Yellow River in China, specifically the Wei River (widths range from 100 m upstream to 3,000 m downstream), as the test bed. In addition, the Wei River has a high sediment yield and concentration, with the sediment load during the flood season accounting for 75% to 94% of the annual total. The Wei River basin (WRB), located between 103.5-110.5°E and 33.5-37.5°N, with daily mean temperatures ranging from 6 to 14°C, is situated in a continental monsoon climatic zone. The average annual precipitation ranges from 450 mm to 700 mm, occurring mostly during the period from May to October, and the annual mean evaporation is between 1000 and 2000 mm (Zhang & Sun, 2015). The precipitation distribution in the WRB has been observed
to decrease from the southeast to the northwest, showing large spatial differences that induce corresponding influences on river discharge. The scarcity of rainfall during the non-flood season leads to small discharge, whereas intense flooding during the flood season results in a dramatic increase in discharge. The Wei River basin is the major region for agriculture, industry, and commerce in Northwest China. However, hydrological stations are scarce for such an important river in China, thereby making it difficult to implement effective measures to prevent floods or drought in advance for sustainable water resources management. Therefore, in this study, we selected six representative river cross-sections in the Wei River to propose a method for monitoring the discharge of poorly gauged rivers. Among the six sections, Tuoshi (TS) and Linjiacun (LJC) are in the upper streams, Weijiabu (WJB) and Xianyang (XY) are in the middle streams whereas Lintong (LT) and Huaxian (HX) are in the lower streams, as shown in Fig. 1.

<figure><img src="image_3.png" /><figcaption>Fig.1 Study area and the representative gauging stations.</figcaption></figure>

## 3. Data and Methods

### 3.1 Data

#### 3.1.1 Satellite images

Landsat 5 Thematic Mapper (TM) and Landsat 8 Operational Land Imager (OLI) images corresponding
to the ground-measurement periods were used (Table 1). The TM sensor has seven bands whereas the OLI sensor has nine bands. For velocity estimation, band 4 (0.76–0.96 µm) of TM and band 5 (0.851–0.879 µm) of OLI, both with a spatial resolution of 30 m, were utilized. Sentinel-2 images with a high spatial resolution (10 m) and high temporal resolution (5 days) are employed for further validation in this study. The Sentinel-2 sensor has 13 bands, and band 8 (10 m) is specifically chosen for velocity estimation. All images were acquired from Google Earth Engine (GEE), a free cloud computing platform established by Google in 2010, enabling users to perform geospatial analysis on a global scale using Google's infrastructure (Gorelick et al., 2017). The cloud cover percentage of the current image is calculated based on GEE, and images with cloud cover <10% are selected. The median filtering technique (dos Santos et al., 2022) is used to handle anomalies in NIR values to smooth the data, remove noise and improve data reliability. We selected 67 Landsat 5 images, 11 Landsat 8 images, and 128 Sentinel-2 images for the stations in the Wei River basin. Taking the Landsat 5 images for example, among all six stations, WJB is the only station covered by all 67 Landsat 5 images. Therefore, it can provide the most sufficient information for the calibration of C/M and velocity as a reference station. LJC and TS are covered by 33 images, and XY, LT, and HX are covered by 34 images. On the other hand, Sentinel-2 images were used for velocity estimation across the stations for further validation.

<table><caption>**Table 1** Satellite images used in this study.</caption><thead><tr><th>Satellite</th><th>Satellite Path Index</th><th>Satellite Row Index</th><th>Corresponding Hydrological Stations</th><th>Numbers of Images</th></tr></thead><tbody><tr><td rowspan="2">Landsat 5 (2007–2009)</td><td>127</td><td>36</td><td>WJB, XY, LT, HX</td><td>34 (17 in flood season, 17 in non-flood season)</td></tr><tr><td>128</td><td>36</td><td>WJB, LJC, TS</td><td>33 (11 in flood season, 22 in non-flood season)</td></tr><tr><td>Landsat 8 (2019–2020)</td><td>127</td><td>36</td><td>WJB, LT</td><td>11 (3 in flood season, 8 in non-flood season)</td></tr><tr><td rowspan="2">Sentinel-2 (2019–2020)</td><td>—</td><td>—</td><td>WJB</td><td>56 (18 in flood season, 38 in non-flood season)</td></tr><tr><td>—</td><td>—</td><td>LT, HX</td><td>29 (9 in flood season, 20 in non-flood season)</td></tr></tbody></table>
— — XY

18 (9 in flood season, 9 in non-flood season)

— — LJC, TS

25 (11 in flood season, 14 in non-flood season)

## 3.1.2 UAV Data

UAVs are suited for capturing aerial images, generating digital surface models (DSM), and extracting high-precision land and water surface elevation information in the Wei River basin, where the flight environment is favorable with stable river bathymetry. Therefore, the Phantom 3 Standard-DJI drone (Table 2) used in this study flew a total of 5 times at 5 sites at the end of June 2019. The river discharges at the five sites on the day of the flight were as follows: TS (30.76 m³/s), LJC (5.53 m³/s), XY (34.1 m³/s), LT (116.5 m³/s), HX (66.2 m³/s). We assume minimal changes in the river's cross-sectional features, i.e., considering the cross-sectional data above the water surface obtained from the UAV flights in 2019 as equivalent to those from 2007–2009 on the same dates. The flight control software used was Pix4D Capture from Pix4D S.A. (https://pix4d.com/). The UAV fly collected 200 to 300 images per river, with the photograph shooting overlap rate set to 90% to ensure the accuracy of the subsequent postprocessing (e.g., stereo image pair, point cloud, and digital surface model (DSM) datasets) using Pix4dMapper Software (Version 2.0.104) (https://www.pix4d.com, Zhao et al., 2022). The generated DSM had a spatial resolution of 2.22 cm when the flight altitude was set to 50 m, and the elevation precision was verified by the total station, which was implemented during the UAV flights.

<table><caption>Table 2 DJI Phantom-3-pro's basic parameters.</caption><thead><tr><th>UAV Product</th><th>Phantom-3-pro</th><th>UAV photo</th></tr></thead><tbody><tr><td>Camera</td><td>FC300X</td><td></td></tr><tr><td>Camera Sensor</td><td>Sony Exmor R CMOS</td><td></td></tr><tr><td>Max Photo Resolution</td><td>4000×3000</td><td></td></tr><tr><td>Max Aperture</td><td>f/2.8</td><td></td></tr><tr><td>Focal Length</td><td>20 mm</td><td></td></tr></tbody></table>

<figure><img src="image_4.png" /></figure>
Field of View (FOV) 94°

Max Flight Altitude 500 m

### 3.1.3 Ground Data

In the Wei River, daily river discharge, flow velocity, water level, and cross-sections of six typical stations were obtained and quality-controlled by the hydrological bureau of the Yellow River Conservancy Commission. The specific time range for data collection is shown in **Table 3**. Among these, the TS station was established in 2004 and the data collection has been continuous since then. The measurement criteria of the data for different years were consistent throughout the entire recorded period. Due to the incomplete continuity of gauged data in the hydrological bureau of the Yellow River Conservancy Commission, contemporaneous data for the satellite overpass dates may not be available. Therefore, we supplemented the missing data based on the gauged counterparts from the hydrological bureau using the cubic spline interpolation method (Table 3). Flow velocity at the reference station was used for calibration of the velocity estimation (referring to Section 3.2.1). In situ water levels were jointly used with mathematical modelling (referring to Section 3.2.2) to estimate cross-sections when there were no corresponding UAV measurements. In situ cross-section and discharge data across the other stations were used for validation purposes only.

<table><caption>**Table 3** Ground-measured data and their corresponding purposes for the six gauging stations (i.e., TS, LJC, WJB, XY, LT, and HX).</caption><thead><tr><th>Stations</th><th>Data</th></tr></thead><tbody><tr><td>TS</td><td>Discharge (2004-2011) (validation)<br>Velocity (2004-2011) (validation)<br>Water level (2007-2009) (calibration)<br>Cross-section (2007-2009) (validation)</td></tr><tr><td>LJC</td><td>Discharge (2001-2011) (validation)</td></tr><tr><td>XY</td><td>Velocity (2001-2011) (validation)</td></tr><tr><td>HX</td><td>Water level (2007-2009) (calibration)<br>Cross-section (2007-2009) (validation)</td></tr><tr><td>LT</td><td>Discharge (2001-2011, 2019-2020) (validation)<br>Velocity (2001-2011, 2019-2020) (validation)<br>Water level (2007-2009, 2019-2020) (calibration)</td></tr></tbody></table>
## 3.2 Methods

A modified C/M method is proposed to overcome the shortcomings of the spectral information-based method and calculate flow velocity. Then a novel algorithm for estimating River Benthymetry with zero ground measurements, termed the RIBA-zero method, is put forward to estimate flow area based on the UAV-monitored above-water cross-section. By multiplying the flow velocity and the flow area the discharge in the ungauged river section is finally derived.

### 3.2.1 Calculation of flow velocity by presenting a modified C/M method

The boundary of the largest water body (a static area with a size of 4 km × 4 km) in all analyzed Landsat /Sentinel-2 images is determined based on the Modified Normalized Difference Water Index (MNDWI) (Xu, 2005), and is termed the Maximum River boundary (MRB). The near-infrared reflectance of all pixels within this area of interest is extracted. The mean of near-infrared bands within MRB (the averaged reflectance value of all pixels within the MRB area) is considered as M, which changes with the flooded area of the river. An increase in river discharge within the MRB enlarges water proportion or wet area, thereby resulting in a decrease in MRB-averaged reflectance. Using MRB-averaged instead of individual wet pixel's reflectance effectively reduces the uncertainties and stabilizes the C/M-V relationship. All hardened urban areas on the images (e.g., road and roof) are extracted, of which the average value of the near-infrared bands on different images is counted as C.

Sediment affects water reflectance. Therefore, Filippucci et al. (2022) developed the application of the sediment correction (i.e., the CMW approach). Based on previous studies, we attempted to calibrate the M value of sediment-laden rivers based on changes in water proportion and water pixel reflectance to better eliminate the influence of sediment (Eq. 1). Subsequently, this study uses Landsat 5/8 images and in situ velocity data at the single hydrological station (i.e., reference station) to establish a regional C/M-V relationship in flood and non-flood seasons, respectively, for a poorly gauged basin (Eq. 2). We assume that the relationship of Eq. 2 is universal, that is, the coefficients a and b are the same at both
reference and ungauged stations. To make this assumption feasible, we propose Eqs. 3 & 4 to correct
M at ungauged stations by establishing the relationship between the M values of the reference and other
ungauged stations to make the coefficients a & b static in Eq. 2 when applied to other stations. Tarpanelli
et al., (2013) showed that stations along the same river have distinct C/M-V formulas, highlighting the
need for each station to gather a substantial number of gauged V observations before applying the C/M
method. The extensive data demand of the C/M method makes it challenging to apply to large-scale
rivers. To dramatically reduce the data requirement and therefore make the C/M method applicable to
large-scale data-scarce areas, we made improvements using Eqs. 2-4. In this study, Eq. 2 represents
the hypothesis and we ensure its validity by Eqs. 3 & 4 before applying Eq. 2 across the entire river.
Here, we plotted the data pairs of C/M and V at all stations along a river before (Fig.2a) and after (Fig.
2b) correction with Eqs. 3 & 4. It was found that R² increased from 0.41 (before correction) to 0.67 (after
correction). This justifies the assumption that after correction all stations along the same river adhere to
the same C/M-V relationship.

<figure><img src="image_5.png" /><figcaption>Fig. 2 C/M-V relationship diagrams before (a) and after (b) correction for the stations in the Wei River
basin. Red lines and colored dash lines stand for trends of the entire and individual stations that
determine the parameters in Eq. (2), respectively. Note that ground velocity measurements have been
normalized due to the confidentiality of the data using $V_{\text{nor}} = (V - V_{\text{min}})/(V_{\text{max}} - V_{\text{min}})$, where $V_{\text{nor}}$ is the
normalized velocity, $V$ is the original velocity (m/s), $V_{\text{min}}$ is the minimum velocity (m/s), and $V_{\text{max}}$ is the
max velocity (m/s).</figcaption></figure>
The MRB represents the maximum boundary where water can reach, generally with water in the center and soil on the sides, and the reflectance of soil within the MRB changes less than water. Therefore, the relationship between the reflectance of MRBs is mostly dominated by the reflectance of the water within the MRBs due to the change of water color (sediment or pollution). Because the central pixel of the MRB is heavily dominated by water, to establish a reliable relationship between the reference and ungauged MRB, we selected a single central pixel (water pixel) from each of them as a proxy for the relationship between the two. Moreover, after reflectance correction for all MRBs with Eqs. 3 & 4, the same relationship can be applied to all the MRBs in a hydrologically-connected river within a similar climatological setting.

Therefore, we calibrate c and d based on the values of the central pixels in the MRB of the reference station and ungauged station first, and then apply c and d to each pixel in the MRB to calculate the corrected pixel M value of the ungauged station based on Eq. 4. The reason is that there exists spatial-spectral heterogeneity in reflectance between ungauged and reference stations because water pollution and soil erosion induced by human activities varies among different reaches and regions in a basin. The skillful corrections by using Eqs. 3 & 4 ensure that Eq. 2 has uniform coefficients, thereby avoiding the recalibration of coefficients a and b for every ungauged station. In the end, we use Eq. 5 for flow velocity calculation at ungauged river sections (Fig. 3).

$$M = M' - (M'_{w} - M_{w}) \times W_{w} \qquad (1)$$

where M is the time series of the modified MRB-averaged NIR reflectance value of the river section, M' is the time series of the original MRB-averaged NIR reflectance of the river section, M'w is the time series of the NIR reflectance of the water pixels which are located at the inner part of the river, with water always present even during low flows, Mw is the minimum value of the time series of the NIR reflectance of the selected permanent water pixels of the river section (i.e., a single value), Ww is the time series of the water body proportion of the river section, calculated initially by extracting water bodies using the MNDWI index from remote sensing imagery. Specifically, the MNDWI index within MRB is calculated based on the green band and shortwave infrared band. Then the water body threshold is automatically extracted using the OTSU method (Dong et al., 2021; Otsu, 1979) to identify the water body portion.
$$274 \qquad V = a \cdot C/M + b \quad (\text{Universal formula based on the reference station}) \qquad (2)$$

where V is the flow velocity of the river section (m/s). Parameters a and b are calibrated based on the long-term time series of flow velocity and M value data at the reference station.

$$277 \qquad M_{\text{ref,cen}} = c \cdot M_{\text{u,cen}} + d \quad (\text{Correcting the M value at ungauged station to the reference station basis by using the link of central pixel in MRB between ungauged and reference stations}) \qquad (3)$$

where $M_{\text{ref,cen}}$ and $M_{\text{u,cen}}$ are the M values at central pixels in MRB on the reference and ungauged river stations. Parameters a and b are calibrated by using M values at the central pixels within MRB on the reference and ungauged stations.

$$283 \qquad Mu' = \frac{1}{N} \sum_{p=1}^{N} (c \cdot M_{u,p} + d) \quad (\text{Corrected M value within MRB at the ungauged station based on all pixels in MRB}) \qquad (4)$$

where $Mu'$ is the corrected M value of the ungauged river section, $M_{u,p}$ is the M value of a pixel in MRB on the ungauged station, p represents different pixels in MRB, and N is the number of pixels in MRB.

$$287 \qquad V_u = a \cdot C/Mu' + b \quad (\text{Formula for ungauged stations}) \qquad (5)$$

$V_u$ is the flow velocity of the ungauged river section (m/s).
<figure><img src="image_6.png" /><figcaption>Fig. 3 Flowchart of the modified C/M method.</figcaption></figure>

### 3.2.2 Calculation of flow area using the RIBA-zero method

#### (1) Estimation of the underwater river cross-section

To calculate flow area accurately a dataset of precise underwater cross-sections is essential. Despite the reported remote sensing techniques to monitor the cross-section under the water surface (e.g., through the sonar or green Light Detection And Ranging (LiDAR)), the drone used in the study cannot carry water penetration equipment for bathymetric measurement due to limitations in the payload capacity and power supply. Therefore, another option is to estimate the underwater cross-section by mathematical methods. Many studies introduce curve-fitting approaches for predicting unobserved symmetric cross-sections (Mersel et al., 2013; Schaperow et al., 2019). However, due to the natural curvature of rivers and the asymmetry of cross-sections, simple linear models or symmetric fittings may result in significant uncertainties in estimating flow area. Therefore, further mathematical exploration is needed for flow area estimation. From the viewpoint of mathematics, the cross-section of a river can be
expressed in the form of a power function, and the power value of the function is usually less than 2 (Dingman, 2007) (Fig. 4).

This section developed a new algorithm (termed the RIBA-zero method) for modelling underwater cross-sections. Because there are differences between the shapes of the two sides of a cross-section, the fitting formulas of the distance and water depth of the left and right sides of the cross-section are established respectively based on the UAV-derived DSM (Eq. 6).

$$ \begin{cases} y_1 = a_1(|x_1 - b_1|)^{t_1} + c_1 \\ y_r = a_r(|x_r - b_r|)^{t_r} + c_r \end{cases} \qquad (6) $$

where y represents the water depth of the river (m); x represents the distance from the starting point (m), or the value 0.0 in Fig. 4, and the starting point can also be any value instead of 0.0 if necessary; t represents the power of the cross-section shape; a represents the opening degree of the cross-section, i.e., a = ctg (α) in Fig. 4; b represents the central axis point of the fitting function (b = 0 in Fig. 4) and if necessary b can be other values than 0; the subscript l and r stand for the left and right side of the cross-section, respectively; c represents the minimum water depth (m). c = 0 in Fig. 4, and can be other values than 0 if necessary.

<figure><img src="image_7.png" /><figcaption>Fig. 4 Ideal cross-section or relationship between water depth and distance in Eq. 6 with $t \le 1$ (a) and $t \ge 1$ (b).</figcaption></figure>

In Eq. 6 and Fig. 4, the difference in power value represents the shape of the cross-section of the river. For instance, the curve with $t = 1$ represents a triangle-type cross-section, $t = 1.75$ the 'Lane Type B stable channel', and $t = 2$ a parabola shape cross-section (Dingman, 2007; Bjerklie et al., 2018). So far,
the cross-section formula is established. However, in reality, there will be slope break points in a cross-section that are not as smooth as in Fig. 4 due to the frequent eroding and sedimentation, as shown in Fig. 5.

<figure><img src="image_8.png" /><figcaption>Fig. 5 Example of the slope breakpoint on the cross-section. Thick black curve segment: the segmented section above but closest to the water surface, which is taken as a basis to fit underwater cross-section. To avoid overfitting, as long as possible the thick black curve segment should be selected.</figcaption></figure>

To model the actual across-section in Fig. 5, the slope (dy/dx) and the change rate of slope (d²y/dx²) of the above-water cross-section are calculated and the local maximum of d²y/dx² is used as the slope mutation point for extracting the segment above water cross-section that is closest to the water surface. The relationship between the slope of the cross-section and the starting distance can be formulated as Eq. 7.

$$ \begin{cases} \text{slope}_l = a_l \cdot t_l \cdot (x_r - b_l)^{t_l-1} \\ \text{slope}_r = a_r \cdot t_r \cdot (x_r - b_r)^{t_r-1} \end{cases} \qquad (7) $$

Depending on whether the power of the fitted relationship is less than 1, the segmented sections of the cross-section closest to the water surface are divided into the convex or concave cross-section, and the
underwater cross-section is determined according to three scenarios: (a) Scenario 1: $t_i < 1$ and $t_r < 1$ (Fig. 6a), i.e., both the left and right sides of the cross-section are convex. To fit the underwater cross-section with Eq. 6 using segments that are closest to the water surface, we use the points with $x_l = b_l$ and $x_r = b_r$ (i.e., slope equal to infinite) as the bottom points of the river on two sides; (b) Scenario 2: $t_i > 1$, $t_r < 1$, or $t_i < 1$, $t_r > 1$ (Fig. 6b), i.e., only one side of the cross-section is convex. We detect the bottom point of the cross-section by finding the point with a slope equal to infinite; (c) Scenario 3: $t_i > 1$ and $t_r > 1$ (Fig. 6c), the intersection of the formulas fitted by the left and right sides of the above water cross-section is used as river bottom point, thereby determining the underwater cross-section. Combined with Eq. 6 and Eq. 7, we can model the underwater cross-section using above water cross-section (including water level) produced by the UAV.

<figure><img src="image_9.png" /><figcaption>**Fig. 6** Example of three scenarios when determining the below water cross-section. (a) Scenario 1 & (b) Scenario 2: no intersection is found in fitting the formula of two sides of the cross-section; (c) Scenario 3 intersection is found in fitting the formula of two sides of the cross-section.</figcaption></figure>

## (2) Calculating flow area

After modelling the underwater cross-section with the RIBA-zero method (Eq. 6 and Eq. 7), we can then calculate the flow area (or underwater cross-section area) as shown in purple in Fig. 6. In Fig. 6a&b. The area for the purple polygon is calculated with Eq. 8. In Fig. 6c. The area for the purple polygon is calculated with Eq. 9.

$$ 
\left\{ 
\begin{array}{l}
A = A_{\text{all}} - (A_{\text{left}} + A_{\text{right}}) = |x_r - x_l| \cdot y_l - \left[ \int_{x_l}^{x_{l0}} (a_l(x - b_l)^{t_l} + c_l) + \int_{x_{r0}}^{x_r} (a_r(x - b_r)^{t_r} + c_r) \right] \\
y_l = \left| (a_l(x_l - b_l)^{t_l} + c_l) - (a_l(x_{l0} - b_l)^{t_l} + c_l) \right|
\end{array}
\right. 
\qquad (8)
$$
359

$$ egin{cases} A = A_{\text{all}} - (A_{\text{left}} + A_{\text{right}}) = |x_r - x_l| \cdot y_l - \left[ \left( \int_{x_l}^{x_o} (a_l(x-b_l)^{t_l} + c_l) \right) + \int_{x_o}^{x_r} (a_r(x-b_r)^{t_r} + c_r) \right] \\ y_l = \left| (a_l(x_l-b_l)^{t_l} + c_l) - (a_r(x_r-b_r)^{t_r} + c_r) \right| \end{cases} \quad (9) $$

360

361 where A represents the flow area (m²), Aall represents the whole area of the red rectangle (m²), Aleft and Aright represent the grey area under the curves (m²), xl and xr represent the distance of the left and right sides of the cross-section crossed by the water surface (m). In addition, xl0 and xr0 represent the lowest underwater points on the left and right sides (m); x0 represents the underwater intersection point (m), and yl represents the height difference from the water surface line to the underwater intersection point (m). xl and xr can be directly obtained based on UAV images. In cases where UAV measurements are missing, we calculate the average water-level difference between the ungauged station with UAV flights coupled with mathematical modelling and the reference station with ground water-level measurements at the corresponding date. Then, we can easily estimate the water level for ungauged stations based on the water level data from the reference station hereby to determine xl and xr based on the cross-sections fitted by the UAV. Unavoidably this may introduce uncertainty, which can be effectively reduced by increasing the number of UAV flights.

362

363

364

365

366

367

368

369

370

371

372

373 With the flow velocity and flow area determined above, we can finally calculate the river discharge with Eq. 10.

374

375

$$ Q = A \cdot V \quad (10) $$

376 where Q represents the river discharge (m³/s), A represents the flow area (m²), and V represents the flow velocity (m/s).

377

378

# 4. Results

379

## 4.1 Estimating flow velocity using Landsat 5 images

380

381

382

383

Before determining the CIM-V relationship, the corresponding areas of land and water were extracted to calculate C and M using methods specified in Section 3.2.1. For the calculation of C, a large number of urban pixels on the image were extracted from two Landsat 5 TM images (satellite path No. 127 and No. 128, Table 1). Because the cities in the study area developed alongside rivers, except for most of
384

the urban pixels near rivers we also used some urban pixels far from the river. Fig. 7 shows that the reflectance of the urban area in the near-infrared band is high.

385

<figure><img src="image_10.png" /></figure>

386

**Fig. 7** Extracted boundary for C value calculation from Landsat 5 TM. Sub-figures a & c: image No. 128; b & d: image No. 127. Images a & b are true color images, and c & d are near-infrared images. Red circles stand for urban areas selected with high reflectivity in the near-infrared band.

387

388

389

For the calculation of the *M* value, this study visually extracts the boundary of the water body area, or MRB recorded in the image during 2001–2011 based on Landsat 5 true-color images. Then the average reflectance of the pixel within the MRB was used as the *M* value. The maximum *M* value is 0.34 and the minimum one is 0.12. The MRB-based *M* values change significantly over time as water and land have different spectral characteristics in the NIR band over time. When the water level rises, the number of water pixels (land pixels) increases (decreases), and vice versa. Taking partial images of LT as an example (Fig. 8), during non-flood periods or low-flow periods, the river is clear, with a low reflectance (dark blue in Fig. 8c, d). During the flood period, especially during the peak flow, the river carries a large

390

391

392

393

394

395

396

397
amount of sediment, thereby resulting in reflectance values that are close to those of the adjacent shore

(Fig. 8a, b).

<figure><img src="image_11.png" /></figure>

400

Fig. 8 The difference in reflectance when the river carries sediment (a, b) and when it is clear (c, d). The red line is the maximum river boundary (MRB) of the LT (a downward hydrological station in Figs. 1 & 7) based on TM images. The yellow square represents the selected permanent water pixels. Blue and green colors show the value of the Landsat surface reflectance.

Therefore, we introduced Eq. 1 to correct the original *M* value (0.11–0.27), and the modified *M* value ranges from 0.09 to 0.24. After obtaining the values of *C* and modified *M* at each station, we established the relationship between the flood season and non-flood season by the least squares regression method based on 67 pairs of *C/M* and *V* of the reference station (WJB), as shown in Fig. 9.
<figure><img src="image_12.png" /><figcaption>**Fig. 9** The relationship between the C/M and velocity in different seasons for WJB. The R² value is used to evaluate the correlation between the flow velocity and the C/M. The red and blue areas represent the 95% confidence interval that is calculated using the least squares method.</figcaption></figure>

According to **Fig. 9**, the relationships between the C/M and V of the reference station (WJB) in flood and non-flood seasons are different. During the flood season, the slope of the C/M-V curve is 0.66 with an R² of 0.84 whereas during the non-flood season, the slope of the C/M-V curve is 0.57 with an R² of 0.65. According to the fitting formula, for all periods, the average relative error of the estimated flow velocity for the reference station (WJB) is 20.10%, which indicates a high accuracy. The average relative error in flood season and non-flood season is 20.85% and 19.56%, respectively.

The fitting relationship between M values at the reference station and ungauged river section is established according to Eq. 3 in the flood and non-flood seasons, as shown in **Table 4**. In flood seasons, the highest R² fitted by the TS station is 0.92, and the lowest R² fitted by the LJC is 0.60; in the non-flood season, the highest R² fitted by the TS station is 0.75, and the lowest R² fitted by the HX is 0.56. The reliability of the M value relationship among the five sites varies primarily due to the following reasons: variations in the water color and turbidity among different sites and vegetation growth when surfaces are exposed during the low-water period.

**Table 4** Fitting of M value relationship between ungauged river sections and the reference station (WJB).
<table><thead><tr><th>Sections</th><th>Average river width (m)</th><th>Flood Season Expressions</th><th>R<sup>2</sup></th><th>F-Sig</th><th>Non-flood season Expressions</th><th>R<sup>2</sup></th><th>F-Sig</th></tr></thead><tbody><tr><td>LT</td><td>158.31</td><td>M<sub>WJB</sub>=1.31M<sub>LT</sub>-0.0582</td><td>0.91</td><td>1.35E-7</td><td>M<sub>WJB</sub>=0.89M<sub>LT</sub>+0.0043</td><td>0.63</td><td>6.09E-3</td></tr><tr><td>HX</td><td>150.27</td><td>M<sub>WJB</sub>=1.40M<sub>HX</sub>-0.1087</td><td>0.87</td><td>1.09E-6</td><td>M<sub>WJB</sub>=0.51M<sub>HX</sub>+0.0558</td><td>0.56</td><td>7.65E-3</td></tr><tr><td>XY</td><td>63.92</td><td>M<sub>WJB</sub>=1.106M<sub>XY</sub>-0.0457</td><td>0.88</td><td>5.72E-6</td><td>M<sub>WJB</sub>=0.96M<sub>XY</sub>+0.0054</td><td>0.65</td><td>4.65 E-3</td></tr><tr><td>LJC</td><td>30.45</td><td>M<sub>WJB</sub>=0.74M<sub>LJC</sub>+0.028</td><td>0.60</td><td>5.15E-3</td><td>M<sub>WJB</sub>=0.90M<sub>LJC</sub>-0.0148</td><td>0.70</td><td>1.21E-6</td></tr><tr><td>TS</td><td>46.08</td><td>M<sub>WJB</sub>=0.90M<sub>TS</sub>-0.0169</td><td>0.92</td><td>1.18E-5</td><td>M<sub>WJB</sub>=0.77M<sub>TS</sub>+0.0183</td><td>0.75</td><td>1.96E-6</td></tr></tbody></table>

After obtaining the fitting relationship between the *M* value of each ungauged station and the reference station, we can calculate flow velocity with Eq. 4 for each ungauged river section. We used satellite images to extract the *M* value of the corresponding station at different times, calculate flow velocity for the ungauged stations and perform error analysis with the gauged velocity, as shown in **Fig. 10**.

<figure><img src="image_13.png" /><figcaption>**Fig. 10** The scatterplot between the estimated velocity and ground-measured counterpart. The grey dashed line represents the linear fit line between the ground-measured velocities (V<sub>G</sub>) and the estimated velocities (V<sub>E</sub>) at TS, LJC, and XY stations, where the average width of the river is less than 90 m. The blue dashed line represents the linear fit line between the V<sub>G</sub> and V<sub>E</sub> at LT and HX stations, where the average width of the river is larger than 90 m, and the black lines represent the 1:1 line. (a) flood season, (b) non-flood season. RW means the average river width. Note that the velocity has been normalized.</figcaption></figure>
The accuracy of flow velocity estimation differs among stations with different river widths. The five stations have been divided into two categories based on the multi-year average river widths derived from the hydrological bureau of the Yellow River Conservancy Commission. One category is for sites with river widths greater than 90 m (which is more than three times the image resolution), including the LT and HX stations. The other category is for sites with river widths less than 90 m (which is less than three times the image resolution), including the XY, LJC, and TS stations (Fig. 10). The average relative error of the estimated velocity for stations with river widths greater than 90 m (LT and HX) in the entire period is 21.35%, and the rRMSE is 0.23, indicating the relationship is reliable for large rivers. However, for stations with river widths less than 90 m, the average relative error increased by 17.05%, and the rRMSE increased by 0.18, indicating slightly worse estimation performance. In flood season, the average relative error for river width greater than 90 m/ less than 90 m is 16.54%/32.04%, and the rRMSE is 0.18/0.33, respectively. In non-flood season, the average relative error of stations with river width greater than 90 m/ less than 90 m is 26.31%/42.27%, and the rRMSE is 0.30/0.48 (Table 5). It can be concluded from the above results that our modified CIM method provides better estimation in flood season as opposed to non-flood season. Moreover, our method performs better on rivers wider than 90 m and exhibits an acceptable estimation accuracy for rivers narrower than 90 m.

<table><caption>**Table 5** The relative root-mean-square error (rRMSE) of the estimated flow velocity of all stations in the middle Yellow River.</caption><thead><tr><th rowspan="2">Period</th><th colspan="2">River width>90 m</th><th colspan="3">River width<90 m</th></tr><tr><th>LT</th><th>HX</th><th>XY</th><th>LJC</th><th>TS</th></tr></thead><tbody><tr><td>Flood season</td><td>0.19</td><td>0.17</td><td>0.23</td><td>0.36</td><td>0.40</td></tr><tr><td>Non-flood season</td><td>0.27</td><td>0.33</td><td>0.41</td><td>0.58</td><td>0.46</td></tr><tr><td>All periods</td><td>0.23</td><td>0.24</td><td>0.32</td><td>0.47</td><td>0.44</td></tr></tbody></table>

Among all stations, the result of LT and HX (river widths > 90 m) performed best both in flood and non-flood seasons, with rRMSE values of 0.19/0.17 and 0.27/0.33, respectively, whereas LJC and TS (RW < 90 m) performed worst with rRMSE values of 0.36/0.40 in the flood season and 0.58/0.46 in the non-flood season. For station XY, although the average width of the river is less than 90 m, it expands to over 90 m during the flood season, thereby showing better performance during this period. It can be concluded that when using the method in this study to estimate flow velocity based on Landsat 5 TM
462 (with a resolution of 30 m), the wider the river is, the higher the accuracy of flow velocity estimation.

463 **4.2 Estimating underwater cross-section area**

464 After validating the formulas for flow velocity calculation, we then evaluated the flow area estimation
465 model based on ground-measured cross-section data (2007-2009) from the hydrological bureau of the
466 Yellow River Conservancy Commission, which was conducted in the middle Yellow River, as shown in
467 **Fig. 11**. The measured flow area was used as a benchmark to analyze the accuracy of the modelled
468 flow area at each station.
<figure><img src="image_14.png" /><figcaption>**Fig. 11** Typical underwater cross-section estimations in 2007, 2008, and 2009 over LT(a-c), HX(d-f), XY (g-i), LJC (j-l), and TS (m-o).</figcaption></figure>

According to **Fig. 11**, the relative error of the estimated flow area of all stations in the middle Yellow River regardless of the season is 36.02%, the rRMSE is 0.37, and the relationship between the ground-measured and the estimated flow area was tested by using the F-test method with a significance of 3.08E-33 (less than 0.01). In the flood season, the relative error of the estimated flow area of all stations
in the middle Yellow River is 32.46%, the rRMSE is 0.35, and the relationship between the ground-measured and simulated cross-section is validated by the F-test with a significance of 6.02E-11 (less than 0.01). In the non-flood season, the relative error is 37.54%, the rRMSE is 0.36, and the relationship between the ground-measured and estimated flow area is validated by the F-test with a significance of 8.58E-24 (less than 0.01). The accuracy of estimating the flow area in the flood season is higher than that in the non-flood season. Among all the five stations, the estimated flow area error of XY (2009), where the shape of the cross-section (the underwater section) is most similar to a parabola, is the lowest, and the rRMSE is 0.22. In TS, where the shape of the cross-section is narrow and the riverbed is flat, the rRMSE is 0.56, which is the largest among the five stations. In other words, errors are positively correlated with the degree of bias (from an ideal parabola curve) of the actual underwater section of the cross-section. Based on the above analysis, it can be concluded that the RIBA-zero method is more suited for parabola-shaped cross-sections.

## 4.3 Estimating river discharge using single satellite sources

After deriving the flow area and flow velocity, we then used Eq. 10 to calculate the river discharge at each station. Results are shown in **Fig. 12** and **Table 6**. The average relative error of the estimated discharge for stations with river widths greater than 90 m in the entire period is ~ 30%, the Nash–Sutcliffe efficiency coefficient (NSE) is ~ 0.9, and the rRMSE is ~ 0.3, indicating that the relationship is reliable. For stations with river widths less than 90 m, the average relative error increased by 30%, the NSE decreased by 0.1, and the rRMSE increased by 0.2, indicating slightly worse estimation performance. Among them, the high average relative error is due to the LJC station's average river width being only 30.45 m. At certain times during the non-flood season, the river width even falls below 30m, leading to a significant bias in the estimation of low flows. When excluding the LJC station, that is, when only considering rivers with widths ranging from 40-90 m, NSE/rRMSE is 0.79/0.43. It can be concluded that the developed approaches are more suited for rivers wider than 90 m, but also perform well in rivers ranging from 40-100 m. In the flood season, the range of NSE is 0.91-0.96 and the range of rRMSE is 0.20-0.34 whereas in the non-flood season, the range of NSE and rRMSE is 0.42-0.69 and 0.31-0.44, respectively. It can be concluded that the discharge estimates are more accurate in the flood season than in the non-flood season.
Among all the stations, the flow estimation accuracy is better at the HX and LT stations, followed by XY, whereas LJC and TS have the lowest accuracy. This is possibly due to constraints imposed by the spatial resolution of the images in velocity estimation for narrow rivers and uncertainties caused by the irregular-shaped cross-sections. For example, the velocity estimation accuracy at the TS station is better than at the LJC station, but due to the irregular shape of the TS section, the errors in cross-section estimation are larger, thereby resulting in the greatest error in discharge estimation (Fig. 12). Additionally, the accuracy of LJC flow velocity estimation is low, but its flow estimation accuracy is high. This is because the overall velocity at the LJC station is overestimated whereas the flow area is underestimated. The joint error in discharge estimation combining velocity and flow area has been offset. Therefore, the performance of the estimates is still reliable.
<figure><img src="image_15.png" /><figcaption>**Fig. 12** Validation of the discharge estimates against the ground-measured discharge. (a) flood season, (b) non-flood season, the grey dashed line represents the linear fit line of ground-measured discharge ($D_G$) and estimated discharge ($D_E$) at TS, LJC, and XY stations where the average width of the river is less than 90 m. The blue dashed line represents the linear fit line between the $D_G$ and $D_E$ at LT and HX stations, where the average width of the river is larger than 90 m and the black lines represent the 1:1 line, (c-g) discharge hydrograph of LT (c), HX (d), XY (e), LJC (f), and TS (g). Sparse discharge estimates are due to cloud contamination during high flow periods, coupled with the long revisiting time of the satellite. Only part of the discharge estimates were shown with concurrent in situ data during 2007-2009. Note that ground discharge measurements have been normalized using $Q_{\text{nor}} = (Q - Q_{\text{min}})/(Q_{\text{max}} - Q_{\text{min}})$,</figcaption></figure>
where $Q_{nor}$ is the normalized discharge, $Q$ is the original discharge, $Q_{min}$ is the minimum discharge, and $Q_{max}$ is the max discharge.

<table><caption>**Table 6** Performance of discharge estimates at five sites in the Wei River basin (note that results of flood and non-flood seasons were aggregated due to limited data points).</caption><thead><tr><th>River widths</th><th>Site</th><th>No. of samples</th><th>rRMSE (-)</th><th>NSE (-)</th></tr></thead><tbody><tr><td rowspan="2">&gt; 90 m</td><td>LT</td><td>16</td><td>0.31</td><td>0.87</td></tr><tr><td>HX</td><td>15</td><td>0.26</td><td>0.95</td></tr><tr><td rowspan="3">&lt; 90 m</td><td>XY</td><td>16</td><td>0.37</td><td>0.90</td></tr><tr><td>LJC</td><td>11</td><td>0.41</td><td>0.95</td></tr><tr><td>TS</td><td>11</td><td>0.49</td><td>0.67</td></tr></tbody></table>

## 4.4 Densifying discharge estimates using Landsat and Sentinel-2 images

To complement the discharge estimation in Section 4.3 using a single satellite mission (i.e., Landsat 5), next we supplemented a scenario incorporating Landsat 8 with Sentinel-2 to densify the discharge estimates at the LT station during 2019-2020. To explore the added value of an additional satellite platform in discharge estimation, we plotted the estimated discharge against the gauged counterpart with a table showing the corresponding statistics for Landsat 8 and Sentinel-2, respectively.

Overall, it was found that with a temporal resolution of 5 days, the number of Sentinel-2 images was almost three times that of the Landsat 8 (29 vs 11), thereby proving more discharge estimates in the concurrent dates than its Landsat counterpart. The joint use of the two missions further contributed to the increasing number of discharge estimates, particularly in the non-flood season. The NSE was 0.87 for Landsat, 0.86 for Sentinel-2 and 0.87 for their combination (Fig 13 and Table 7).

However, given the limited useable images in the flood season due to the prevalent cloud contamination, it is not surprising to find that the joint use of Landsat and Sentinel-2 still missed many of the peaks during floods. This analysis confirms the added value of Sentinel-2 in discharge estimates but also necessitates the synergy of multisource remote sensing platforms in discharge estimation for further application (e.g., flood forecasting).
<figure><img src="image_16.png" /><figcaption>**Fig. 13** The joint use of Landsat 8 and Sentinel-2 in discharge estimation at the LT station during 2019-2020. Note that the discharge has been normalized due to the confidentiality of the in situ data.</figcaption></figure>

<table><caption>**Table 7** Performance of discharge estimates using Landsat 8 and Sentinel-2 at LT in the Wei River basin during 2019-2020 (note that results of flood and non-flood seasons were aggregated due to limited data points).</caption><thead><tr><th>Mission</th><th>No. of samples</th><th>rRMSE<br>(-)</th><th>NSE<br>(-)</th></tr></thead><tbody><tr><td>Landsat 8</td><td>11</td><td>0.31</td><td>0.87</td></tr><tr><td>Sentinel-2</td><td>29</td><td>0.55</td><td>0.86</td></tr><tr><td>Landsat 8 + Sentinel-2</td><td>40</td><td>0.53</td><td>0.87</td></tr></tbody></table>

# 5. Discussion

## 5.1 Comparison of flow velocity estimation approaches

The original C/M method requires multiple hydrological station data to determine a regional relationship
between C/M and V. Based on Landsat 5 TM data, this study proposes a modified C/M method that requires only in situ velocity and water level data from one station to enable the application to other images, thereby relieving the constraints on the spatial scale. Due to the differences in river morphology caused by human activities and other factors (Yu et al., 2016; Hoitink et al., 2017), the location of the M pixel determined based on the one-pixel method is difficult to extend to different locations of rivers in the study area (Tarpanelli et al., 2013; Tarpanelli et al., 2015). To overcome this, our method used the mean value of the MRB as opposed to one pixel, which is quite suited for estimating the future discharge of the reference station. The R² of C/M-V was 0.84 in flood seasons and 0.65 in non-flood seasons, respectively, when Landsat 5 TM images were used in this study, which is consistent with the finding reported by Tarpanelli et al. (2015). This justified the reliability of flow velocity estimation in this study when C/M methods were used. To improve the applicability of the study, we also attempted to change the reference stations to LT and HX (after the change, TS and LJC did not have concurrent images. Therefore, only XY, LT, HX, and WJB stations were discussed). We found that the range of rRMSE variation for flow velocity at XY, LT, HX, and WJB stations over the entire period ranged from 0.18 to 0.42 (not shown), demonstrating the stability of our method.

Because of different conditions, such as image lighting conditions, the relationship between C/M and V based on a single station is difficult to apply to different stations in the study area. In our modified method, the relationship is beyond the spatial limitations, thereby enabling the application across diverse stations. To explore the performance of the modified relationship, the accuracy of the original C/M method (not considering the impact of sediments) and the modified counterpart for estimating the flow velocity for ungauged sites was analyzed (Table 8).

**Table 8** Comparison of average relative errors of the modified C/M and original C/M in calculating the flow velocity.

<table><thead><tr><th>Method</th><th>XY</th><th>LT</th><th>HX</th><th>LJC</th><th>TS</th><th>Mean</th></tr></thead><tbody><tr><td>Modified C/M</td><td>28.48%</td><td>19.89%</td><td>23.12%</td><td>48.67%</td><td>38.23%</td><td>31.68%</td></tr><tr><td>Original C/M</td><td>41.05%</td><td>35.78%</td><td>55.95%</td><td>49.05%</td><td>58.53%</td><td>48.07%</td></tr></tbody></table>
According to **Table 8**, in the case of the original method, the average relative error (RE) of all stations is 48.07%, and in the station with the worst estimation (TS), the RE reaches 58.53%. In contrast, for the modified method, the accuracy of the flow velocity estimation is significantly improved with the average RE of all stations reaching 31.68%—decreased by more than 16% at all stations and by more than 30% at HX. Therefore, the modified method shows better performance for velocity estimation at ungauged river sections.

Our method can be well applied in Landsat 5 TM images for estimating flow velocity. To verify the reliability of the correlation between C/M and V, we took WJB as the reference station and used Landsat 5 images to estimate velocity in flood and non-flood seasons, respectively. Among the 28 sets of data in the flood season, 23 sets were selected to fit the relationship between C/M and V, showing that the average R² was 0.84. In the non-flood season, 30/39 sets of data were selected to fit the C/M-V relationship showing that the average R² was 0.65. This proved that the use of Landsat 5 images can also detect the C/M-V relationship, which is in accordance with the study of Li et al. (2019) that confirmed the superiority of the mean value of multiple pixels over a single pixel in generating stable discharge estimates.

Furthermore, results show that the wider the river is, the higher the accuracy of flow velocity estimation due to the impact of image resolution. To investigate the impact of spatial resolution on the results, Sentinel-2 images were also employed in this study. The flow velocity of each station was calculated using the same method (Table 9).

<table><caption>**Table 9** rRMSE (relative root-mean-square error) of the flow velocity of all stations in the middle Yellow River calculated based on Sentinel-2.</caption><thead><tr><th rowspan="2">Period</th><th colspan="2">River width>90 m</th><th colspan="3">River width<90 m</th></tr><tr><th>LT</th><th>HX</th><th>XY</th><th>LJC</th><th>TS</th></tr></thead><tbody><tr><td>Flood season</td><td>0.18</td><td>0.18</td><td>0.23</td><td>0.29</td><td>0.29</td></tr><tr><td>Non-flood season</td><td>0.19</td><td>0.17</td><td>0.34</td><td>0.34</td><td>0.35</td></tr><tr><td>All periods</td><td>0.19</td><td>0.18</td><td>0.27</td><td>0.32</td><td>0.32</td></tr></tbody></table>

Compared with Landsat 5 TM, estimating the flow velocity using Sentinel-2 images decreases the average relative error by 9.53%, and lowers the rRMSE by 0.08 on average (Table 9). When switching
to Sentinel-2 data sources, even the XY, LJC, and TS cross-sections within river widths ranging from 30 to 90 m show good performance, with the overall rRMSE ranging from 0.27 to 0.32 across all periods. It can be concluded that our modified C/M method provides better flow velocity estimation when using Sentinel-2 data with a resolution of 10 m compared to 30 m Landsat 5 TM images. Moreover, the developed approach performs better on rivers wider than three times the satellite resolution. Previous studies have found that river width extraction is more accurate at more than 3 times the satellite resolution (Allen and Pavelsky, 2018), which supports our findings. Tarpanelli et al. (2019) also believed that Sentinel-2 have the potential to estimate the velocity and discharge using the C/M method for narrow rivers (width <100 m).

## 5.2 Uncertainties in cross-section estimation

In Section 4.2, we postulated that the RIBA-zero method is much more suited for the parabola-shaped cross-sections while having difficulties in precisely determining the flow area of the trapezoidal cross-section. To further confirm our hypothesis, a parameter for describing the shape of the cross-section was designed to determine what factor would influence the performance of our method. In RIBA-zero, the slope of the above-water cross-section is the key information to estimate the underwater cross-section. Cross-sections that are approximately parabolic in shape can be represented by a parabola (Zhao et al., 2017). After segmenting the cross-section with the slope breakpoint, the segment of the cross-section closest to the water surface was chosen to estimate the underwater cross-section. Its real slope ($S_r$) and the simulated slope in a parabola fitted by the above-water cross-section ($S_f$) are compared. We define $|S_r - S_f|/S_f$ as the parameter describing the degree of the bias of the real cross-section to a parabola-shaped cross-section and explore how this degree influences the estimation of the underwater cross-section by calculating the elevation bias (RMSE) of the estimated underwater cross-section and in situ underwater cross-section at nine different water levels (not shown).

We found that with the increase in the degree of bias (from the parabola curve) of the segment of the cross-section closest to the water surface, the RMSE of the underwater elevation increases ($R^2$=0.72, RMSE=0.20 m), thereby further confirming that our proposed method is much more suited for the parabola-shaped cross-sections in determining the underwater cross-section. This finding is consistent
with Bjerklie et al. (2018), who assumed a parabolic cross-section based on the analysis of the in situ database to infer bed elevation from height-width measurements. This explains why our method performs better in flood seasons than in non-flood seasons. In this study, the average degree of bias during the non-flood season is 0.31 higher than during the flood season.

The trapezoidal cross-section is difficult to estimate using only the above-water cross-section without information on the river bottom point (Schaperow et al., 2019). As discussed in Tarpanelli et al. (2015), the ground elevation of the river bottom in entropy can improve the accuracy of constructing river cross-sections. Garambois & Monnier (2015) proposed a method combining remotely sensed surface elevation and slope, plus one in situ depth measurement to estimate river discharge in the Garonne River (France), and noted that one in situ depth measurement improved the performance of their method.
We also tested our method in determining trapezoidal cross-sections with riverbed elevation at nine different water levels, and the result is consistent with the studies of Garambois & Monnier (2015) and Tarpanelli et al. (2015). After adding underwater elevations to determine the underwater cross-section area, the rRMSE of the nine estimates is 0.30, which is 6.30% of the results (both flow area and discharge) derived without using underwater elevations.

## 5.3 Analysis of the river discharge errors and uncertainties

In this study, we established a relationship between velocity and C/M for flow velocity estimation, then combined it with the RIBA-zero method that estimates the cross-section to derive river discharge. We also attempted to establish a direct relationship between C/M and river discharge at the reference station, with R² reaching 0.72 in the flood season and 0.18 in the non-flood season, which is significantly lower than the relationship established between C/M and V (R² is 0.83 in the flood season and is 0.64 in the non-flood season). This finding is consistent with Tarpanelli et al. (2013) who showed that the correlation between the C/M and the velocity is higher than that between the C/M and discharge at the Italy Piacenza and Cremona gauging stations. Therefore, using V to build relationships is more reliable. Moreover, as discussed by Yu et al. (2015) and Yang et al. (2024), the spectral reflectance increases with the suspended sediment content in the river. In sediment-laden rivers, high concentrations of suspended sediment typically occur during the flood season, especially during high-flow periods. This
653 interferes with the original C/M-V and C/M-Q relationship. To demonstrate this point, we also compared
654 the relationship between C/M and river velocity before and after eliminating the influence of sediment at
655 reference stations based on Eqs. 1-5. The fitting results after eliminating the influence of sediment were
656 significantly better than those without considering sediment effects (R² is 0.10 in the flood season and
657 0.30 in the non-flood season). This indicates the reliability of our modified method in sediment-laden
658 rivers. Filippucci et al. (2022) also developed the CMW approach using the pixel within the water body
659 to eliminate the influence of sediment in Italy and found that the performance of the CMW method
660 improved significantly compared to the C/M method, demonstrating the importance of considering
661 sediment in the estimation of river discharge.

662 Apart from estimating river discharge through spectral information, there are also many remote sensing-
663 based discharge estimation methods based on river width and water level that can be readily observed
664 by satellites (Zhao et al., 2019; Huang et al., 2018). These methods are difficult to apply to rivers
665 narrower than the satellite resolution because the uncertainty of discharge estimation for small rivers
666 (widths less than 100m) may increase correspondingly. Gleason & Smith (2014) proposed the at-many-
667 stations hydraulic geometry (AMHG) method for river discharge estimation in ungauged rivers based on
668 Landsat 5 TM images and tested their method over multiple rivers in the United States, Canada, and
669 China. The rRMSE was within 0.20-0.30. However, as stated by Brinkerhoff et al. (2020) and Feng et al.
670 (2019), if the river width is less than the resolution of the satellite, the AMHG method will fail. Similarly,
671 the method proposed in this study is also constrained by resolution, making it difficult to estimate rivers
672 narrower than a single pixel. It is a pertinent topic for further investigation of how to decompose individual
673 water-land mixed pixels smaller than the satellite resolution. However, the method proposed in this study
674 demonstrates promising potential for estimating the discharge of medium-sized rivers (40-100 m), with
675 an average rRMSE of 0.42 and NSE of 0.84 (Table 6). Because the C/M method depends on the
676 reflectance of a river, its limitation due to image resolution is relatively smaller compared to methods
677 relying on river width calculation. Li et al. (2019) exploited spectral information by the C/M ratio in
678 estimating the discharge of rivers with widths less than 100 m in China, and the NSE value of their result
679 is 0.45, thereby indicating the potential of using the C/M-based method to monitor the discharge of
680 medium-sized rivers poorly gauged rivers. In addition, when floods exceed the bank range, the
uncertainty of discharge estimates may increase dramatically for many methods whereas the modified C/M method, based on the reflectance of water and soil, can still perform well. To keep the results more precise, the method is suggested to be used at rivers with similar climatic settings and flow regimes and it will benefit from the use of very high-resolution satellite images (e.g., PlanetScope (Flores et al. 2024; Harlan et al. 2023)).

## 6. Conclusions

The present study proposed a framework by combining satellites and UAVs to estimate the river discharge of poorly gauged rivers, catering to both sediment-laden and non-sediment-laden rivers. First, a modified C/M method for velocity estimation was proposed, and then a new method to estimate river bathymetry at poorly gauged rivers (i.e., the RIBA-zero method) was presented. Both methods were validated by ground observations with the following results.

(1) The modified C/M method, using MRB to calculate M values and considering the spatial heterogeneity of river spectral information, can establish a regional relationship between C/M and flow velocity based on only one hydrological station. The estimated velocity rRMSE (relative error) of stations with river width greater than 90 m in the Wei River is 0.23 (21.35%) in the entire year. However, for stations with river widths less than 90 m, the rRMSE increased by 0.18, and the average relative error increased by 17.05%. This method performs better on rivers wider than 90 m and exhibits higher accuracy during the flood season. An attempt to establish a direct relationship between C/M and river discharge shows that the R² is 0.72 in flood seasons and 0.18 in non-flood seasons, which is significantly lower than the relationship established between C/M and V (R² is 0.84 and 0.65 in different seasons). This suggests the relationship of C/M-V as opposed to C/M-Q is recommended when the modified C/M method is applied in practice. In addition, this method can provide better flow velocity estimation when replacing 30 m Landsat 5 TM with 10 m-resolution Sentinel-2 imageries, yielding an overall rRMSE ranging from 0.27 to 0.32 for rivers less than 90 m across all periods.

(2) The RIBA-zero method for determining the bathymetry is more suited for the parabola-shaped cross-sections (rRMSE is 0.22 and the relative error is 12.81%). The degree of bias (from the
parabola curve) of the segment of the cross-section that is the closest to the water surface is the key factor influencing the performance of RIBA-zero in estimating the underwater cross-section. This is the reason for the better performance of RIBA-zero in the flood season than in its non-flood counterpart. In the non-flood season, the average degree of bias (from the parabola curve) is 0.31 higher than that of the flood season. After adding the elevation of the river bottom to determine the underwater cross-section area, the rRMSE is 0.30, which is 6.30% of the counterpart derived without using the elevation of the river bottom.

(3) The modified C/M plus RIBA-zero method demonstrates higher accuracy for rivers wider than three times the satellite resolution, with an estimated discharge rRMSE of 0.29 and an NSE of 0.91 over the entire period. When applying our method to rivers smaller than three times the resolution, the performance metrics are also acceptable. Particularly, for medium-sized rivers (40-100 m), there is an encouraging potential for estimating river discharge with an average rRMSE of 0.42 and an NSE of 0.84. Since the C/M method takes the reflectance as its basis, it suffers from fewer constraints imposed by the image resolution compared to methods based on river widths.

Overall, it is feasible to apply the developed method to estimate river discharge with limited in-situ data (at least from one station). However, the method is currently not applicable to totally ungauged basins without any measured flow velocity. Future research is highly necessitated to explore the possibility of extending this method to similar basins with the joint efforts of multisource remote sensing observations.

## Data availability

Landsat 5/8 and Sentinel-2 images were acquired from https://earthengine.google.com. All ground data to reproduce our work, inclusive of water level, discharge, and UAV-derived cross-sectional data are available at https://zenodo.org/records/14119100.

## Acknowledgements

We acknowledge the reviewers and editors for their valuable suggestions and comments that helped improve the quality of this manuscript. This research was jointly supported by the National Key Project for R&D (Grant numbers: 2021YFC3201103), the National Natural Science Foundation Program of
734 China (Grant numbers: 52279004 & U1812401), the Beijing Advanced Innovation Program for Land

735 Surface Science and the 111 Project (B18006).
737 Allen, G.H. & Pavelsky, T.M. (2018). Global extent of rivers and streams. *Science*, **361**(6402),
738 585-587.

739 Angileri, S. E., Conoscenti, C., Hochschild, V., Märker, M., Rotigliano, E., & Agnesi, V. (2016).
740 Water erosion susceptibility mapping by applying stochastic gradient treeboost to the Imera
741 Meridionale river basin (Sicily, Italy). *Geomorphology*, **262**, 61-76.

742 Bjerklie, D. M., Birkett, C. M., Jones, J. W., Carabajal, C., Rover, J. A., Fulton, J. W., & Garambois,
743 P. A. (2018). Satellite remote sensing estimation of river discharge: Application to the Yukon
744 River Alaska. *Journal of hydrology*, **561**, 1000-1018.

745 Bonnema, M. G., Sikder, S., Hossain, F., Durand, M., Gleason, C. J., & Bjerklie, D. M. (2016).
746 Benchmarking wide swath altimetry-based river discharge estimation algorithms for the
747 Ganges river system. *Water Resources Research*, **52**(4), 2439-2461.

748 Brakenridge, R., & Anderson, E. (2006). MODIS-based flood detection, mapping and
749 measurement: the potential for operational hydrological applications. In *Transboundary
750 floods: reducing risks through flood management*, pp. 1-12, Springer, Houten, Dordrecht,
751 Netherlands.

752 Brinkerhoff, C. B., Gleason, C. J., Feng, D. & Lin, P. (2020). Constraining remote river discharge
753 estimation using reach-scale geomorphology. *Water Resources Research*, **56**(11),
754 e2020WR027949.

755 Brombacher, J., Reiche, J., Dijksma, R., & Teuling, A.J. (2020). Near-daily discharge estimation
756 in high latitudes from Sentinel-1 and 2: A case study for the Icelandic Þjórsá river. *Remote
757 Sensing of Environment*, **241**, 111684

758 Caballero, I., & Stumpf, R. P. (2019). Retrieval of nearshore bathymetry from Sentinel-2A and 2B
759 satellites in South Florida coastal waters. *Estuarine, Coastal and Shelf Science*, **226**, 106277.

760 Che, Y.P., Wang, Q., Xie, Z.W., Zhou, L., Li, S.W., Hui, F., Wang, X.Q., Li, B.G., & Ma, Y.T.
761 (2020). Estimation of maize plant height and leaf area index dynamics using an unmanned
762 aerial vehicle with oblique and nadir photography. *Annals of Botany*, **126**(4), 765-773.

763 Colomina, I., & Molina, P. (2014). Unmanned aerial systems for photogrammetry and remote
764 sensing: A review. *ISPRS Journal of Photogrammetry and Remote Sensing*, **92**, 79-97.

765 Coppo Frias, M., Liu, S.X., Mo, X.G., Nielsen, K., Ranndal, H., Jiang, L.G., Ma, J. & Bauer-
766 Gottwein, P. (2023). River hydraulic modeling with ICESat-2 land and water surface
767 elevation. *Hydrology and Earth System Sciences*, **27**(5), 1011-1032.

768 Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D. & Moore, R. (2017). Google Earth
769 Engine: Planetary-scale geospatial analysis for everyone. *Remote Sensing of Environment*,
770 **202**, 18-27.

771 Dingman, S. L. (2007). Analytical derivation of at-a-station hydraulic-geometry relations. *Journal
772 of Hydrology*, **334**(1-2), 17-27.

773 Domeneghetti, A. (2016). On the use of SRTM and altimetry data for flood modeling in data-sparse
774 regions. *Water Resources Research*, **52**(4), 2901-2918.

775 Dong, Z., Wang, G.J., Amankwah, S.O.Y., Wei, X.K., Hu, Y.F. and Feng, A.Q. 2021. Monitoring
776 the summer flooding in the Poyang Lake area of China in 2020 based on Sentinel-1 data and
777 multiple convolutional neural networks. *International Journal of Applied Earth Observation*
778 and Geoinformation 102.

779 dos Santos, E.P., da Silva, D.D., do Amaral, C.H., Fernandes-Filho, E.I., & Dias, R.L.S. (2022). A machine learning approach to reconstruct cloudy affected vegetation indices imagery via data fusion from Sentinel-1 and Landsat 8. *Computers and Electronics in Agriculture*, 194, 106753.

780 Durand, M., Gleason, C.J., Garambois, P.A., Bjerklie, D., Smith, L.C., Roux, H., Rodriguez, E., Bates, P.D., Pavelsky, T.M., Monnier, J., Chen, X., Di Baldassarre, G., Fiset, J.M., Flipo, N., Frasson, R.P.d.M., Fulton, J., Goutal, N., Hossain, F., Humphries, E., Minear, J.T., Mukolwe, M.M., Neal, J.C., Ricci, S., Sanders, B.F., Schumann, G., Schubert, J.E., & Vilmin, L. (2016). An intercomparison of remote sensing river discharge estimation algorithms from measurements of river height, width, and slope. *Water Resources Research*, 52, 4527-4549

781 Durand, M., Gleason, C.J., Pavelsky, T.M., Frasson, R.P.D., Turmon, M., David, C.H., Altenau, E.H., Tebaldi, N., Larnier, K., Monnier, J., Malaterre, P.O., Oubanas, H., Allen, G.H., Astifan, B., Brinkerhoff, C., Bates, P.D., Bjerklie, D., Coss, S., Dudley, R., Fenoglio, L., Garambois, P.A., Getirana, A., Lin, P.R., Margulis, S.A., Matte, P., Minear, J.T., Muhebwa, A., Pan, M., Peters, D., Riggs, R., Sikder, M.S., Simmons, T., Stuurman, C., Taneja, J., Tarpanelli, A., Schulze, K., Tourian, M.J. & Wang, J.D. (2023). A framework for estimating global river discharge from the surface water and ocean topography satellite mission. *Water Resources Research*, 59(4), e2021WR031614.

782 Durand, M., Dai, C., Moortgat, J., Yadav, B., de Moraes Frasson, R.P., Li, Z., Wadkwoski, K., Howat, I., & Pavelsky, T.M. (2024). Using river hypsometry to improve remote sensing of river discharge. *Remote Sensing of Environment*, 315, 114455

783 Erena, M., Atenza, J. F., García-Galiano, S., Domínguez, J. A., & Bernabé, J. M. (2019). Use of drones for the topo-bathymetric monitoring of the reservoirs of the Segura River Basin. *Water*, 11(3), 445.

784 Feng, D.M., Gleason, C.J., Yang, X. and Pavelsky, T.M. (2019). Comparing discharge estimates made via the BAM algorithm in high-order arctic rivers derived solely from optical cubeSat, Landsat, and Sentinel-2 Data. *Water Resources Research*, 55(9), 7753-7771.

785 Filippucci, P., Brocca, L., Bonafoni, S., Saltalippi, C., Wagner, W., & Tarpanelli, A. (2022). Sentinel-2 high-resolution data for river discharge monitoring. *Remote Sensing of Environment*, 281, 113255.

786 Frasson, R.P.D., Turmon, M.J., Durand, M.T. & David, C.H. (2023). Estimating the relative impact of measurement, parameter, and flow law errors on discharge from the surface water and ocean topography mission. *Journal of Hydrometeorology*, 24(3), 425-443.

787 Flores, J.A., Gleason, C.J., Brinkerhoff, C.B., Harlan, M.E., Lummus, M.M., Stearns, L.A., & Feng, D. (2024). Mapping proglacial headwater streams in High Mountain Asia using PlanetScope imagery. *Remote Sensing of Environment*, 306, 114124

788 Garambois, P. A., & Monnier, J. (2015). Inference of effective river properties from remotely sensed observations of water surface. *Advances in Water Resources*, 79, 103-120.

789 Garambois, P. A., Calmant, S., Roux, H., Paris, A., Monnier, J., Finaud-Guyot, P., & Santos da Silva, J. (2017). Hydraulic visibility: Using satellite altimetry to parameterize a hydraulic model of an ungauged reach of a braided river. *Hydrological processes*, 31(4), 756-767.

800 Gleason, C. J., & Smith, L. C. (2014). Toward global mapping of river discharge using satellite images and at-many-stations hydraulic geometry. *Proceedings of the National Academy of*
821 *Sciences*, **111**(13), 4788-4791.

822 Gleason, C.J., Smith, L.C., & Lee, J. (2014). Retrieval of river discharge solely from satellite imagery and at-many-stations hydraulic geometry: Sensitivity to river form and optimization parameters. *Water Resources Research*, **50**, 9604-9619

825 Gleason, C.J., & Durand, M.T. (2020). Remote Sensing of River Discharge: A Review and a Framing for the Discipline. *Remote Sensing*, **12**, 1107.

827 Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., & Moore, R. (2017). Google Earth Engine: Planetary-scale geospatial analysis for everyone. *Remote Sensing of Environment*, **202**, 18-27

830 Hagemann, M.W., Gleason, C.J., & Durand, M.T. (2017). BAM: Bayesian AMHG-Manning Inference of Discharge Using Remotely Sensed Stream Width, Slope, and Height. *Water Resources Research*, **53**, 9692-9707

833 Harlan, M., Gleason, C., Altenau, E., Butman, D., Carter, T., Chu, V., Cooley, S., Dolan, W., Durand, M., & Eidam, E. (2021). Discharge estimation from dense arrays of pressure transducers. *Water Resources Research*, **57**, e2020WR028714

836 Harlan, M.E., Gleason, C.J., Flores, J.A., Langhorst, T.M., & Roy, S. (2023). Mapping and characterizing Arctic beaded streams through high resolution satellite imagery. *Remote Sensing of Environment*, **285**, 113378

839 Hoitink, A. J. F., Wang, Z. B., Vermeulen, B., Huismans, Y., & Kästner, K. (2017). Tidal controls on river delta morphology. *Nature geoscience*, **10**(9), 637.

841 Huang, C., Chen, Y., Zhang, S., & Wu, J. (2018). Detecting, extracting, and monitoring surface water from space using optical sensors: A review. *Reviews of Geophysics*, **56**(2), 333-360.

843 Huang, Q., Long, D., Du, M., Zeng, C., Qiao, G., Li, X., Hou, A., & Hong, Y. (2018). Discharge estimation in high-mountain regions with improved methods using multisource remote sensing: A case study of the Upper Brahmaputra River. *Remote Sensing of Environment*, **219**, 115-134.

847 Huang, Q., Long, D., Du, M., Han, Z., & Han, P. (2020). Daily continuous river discharge estimation for ungauged basins using a hydrologic model calibrated by satellite altimetry: implications for the SWOT mission. *Water Resources Research*, **56**, e2020WR027309.

850 Huang, Q., Long, D., Han, Z., & Han, P. (2022). High-resolution satellite images combined with hydrological modeling derive river discharge for headwaters: A step toward discharge estimation in ungauged basins. *Remote Sensing of Environment*, **277**, 113030.

853 Ishitsuka, Y., Gleason, C.J., Hagemann, M.W., Beighley, E., Allen, G.H., Feng, D., Lin, P., Pan, M., Andreadis, K., & Pavelsky, T.M. (2021). Combining optical remote sensing, McFLI discharge estimation, global hydrologic modeling, and data assimilation to improve daily discharge estimates across an entire large watershed. *Water Resources Research*, **57**, e2020WR027794

858 Jiang, S., Jiang, W., & Wang, L. (2022). Unmanned aerial vehicle-based photogrammetric 3D mapping: A survey of techniques, applications, and challenges. *IEEE Geoscience and Remote Sensing Magazine*, **10**(2), 135-171.

861 Kim, D., Yu, H., Lee, H., Beighley, E., Durand, M., Alsdorf, D.E., & Hwang, E. (2019). Ensemble learning regression for estimating river discharges using satellite altimetry data: Central Congo River as a Test-bed. *Remote Sensing of Environment*, **221**, 741-755.
864 Kim, J. S., Baek, D., Seo, I. W., & Shin, J. (2019). Retrieving shallow stream bathymetry from UAV-assisted RGB imagery using a geospatial regression method. *Geomorphology*, **341**, 102-114.

867 Kim, S., & Sharma, A. (2019). The role of floodplain topography in deriving basin discharge using passive microwave remote sensing. *Water Resources Research*, **55(2)**, 1707-1716.

869 Kwon, S., Seo, I., & Lyu, S. (2023). Investigating mixing patterns of suspended sediment in a river confluence using high-resolution hyperspectral imagery. *Journal of Hydrology*, **620**, 129505.

871 Leon, J. G., Calmant, S., Seyler, F., Bonnet, M. P., Cauhopé, M., Frappart, F., Filizola, N., & Fraizy, P. (2006). Rating curves and estimation of average water depth at the upper Negro River based on satellite altimeter data and modeled discharges. *Journal of hydrology*, **328(3-4)**, 481-496.

875 Legleiter, C. J. (2015). Calibrating remotely sensed river bathymetry in the absence of field measurements: Flow resistance equation-based imaging of river depths (FREEBIRD). *Water Resources Research*, **51(4)**, 2865-2884.

878 Legleiter, C. J., & Harrison, L. R. (2019). Remote sensing of river bathymetry: Evaluating a range of sensors, platforms, and algorithms on the upper Sacramento River, California, USA. *Water Resources Research*, **55(3)**, 2142-2169.

881 Li, H., Li, H., Wang, J., & Hao, X. (2019). Extending the ability of near-infrared images to monitor small river discharge on the northeastern Tibetan Plateau. *Water Resources Research*, **55(11)**, 8404-8421.

884 Lin, H., Cheng, X., Liu, J., Shi, Q., Li, T., Zheng, L., Hou, X., & Du, J. (2024). Estimating river discharge across scales with a novel regional gauging method driven by Sentinel satellite data. *Remote Sensing of Environment*, **311**, 114266.

887 Lin, P. R., Feng, D. M., Gleason, C. J., Pan, M., Brinkerhoff, C. B., Yang, X., Beck, H. E., & Frasson, R. P. D. (2023). Inversion of river discharge from remotely sensed river widths: A critical assessment at three-thousand global river gauges. *Remote Sensing of Environment*, **287**, 113489.

891 Lu, S., Wu, B., Yan, N., Li, F., Wen, M., & Wang, J. (2010). Progress in river runoff monitoring by remote sensing. *Advance in Earth Sciences*, **25(8)**, 820-826.

893 Moramarco, T., Barbetta, S., Bjerklie, D. M., Fulton, J. W., & Tarpanelli, A. (2019). River bathymetry estimate and discharge assessment from remote sensing. *Water Resources Research*, **55(8)**, 6692-6711.

896 Mejia, A. I., & Reed, S. M. (2011). Role of channel and floodplain cross-section geometry in the basin response. *Water Resources Research*, **47(9)**.

898 Mersel, M.K., Smith, L.C., Andreadis, K.M. & Durand, M.T. (2013). Estimation of river depth from remotely sensed hydraulic relationships. *Water Resources Research*, **49(6)**, 3165-3179.

900 Muhebwa, A., Gleason, C.J., Feng, D., & Taneja, J. (2024). Improving discharge predictions in ungauged basins: Harnessing the power of disaggregated data modeling and machine learning. *Water Resources Research*, **60**, e2024WR037122.

903 Otsu, N. 1979. THRESHOLD SELECTION METHOD FROM GRAY-LEVEL HISTOGRAMS. *Ieee Transactions on Systems Man and Cybernetics* 9(1), 62-66.

905 Pan, F., Wang, C., & Xi, X. (2016). Constructing river stage-discharge rating curves using remotely sensed river cross-sectional inundation areas and river bathymetry. *Journal of Hydrology*,
907 540, 670-687.

908 Papa, F., Durand, F., Rossow, W.B., Rahman, A., & Bala, S.K. (2010). Satellite altimeter-derived
909 monthly discharge of the Ganga-Brahmaputra River and its seasonal to interannual variations
910 from 1993 to 2008. *Journal of Geophysical Research*, **115**, C12013.

911 Papa, F., Bala, S.K., Pandey, R.K., Durand, F., Gopalakrishna, V.V., Rahman, A., & Rossow,
912 W.B. (2012). Ganga-Brahmaputra river discharge from Jason-2 radar altimetry: An update to
913 the long-term satellite-derived estimates of continental freshwater forcing flux into the Bay of
914 Bengal. *Journal of Geophysical Research: Oceans*, **117**, C11021.

915 Pavelsky, T.M. (2014). Using width-based rating curves from spatially discontinuous satellite
916 imagery to monitor river discharge. *Hydrological Processes*, **28**, 3035-3040.

917 Smith, L.C., & Pavelsky, T.M. (2008). Estimation of river discharge, propagation speed, and
918 hydraulic geometry from space: Lena River, Siberia. *Water Resources Research*, **44**,
919 W03427.

920 Schaperow, J. R., Li, D., Margulis, S. A., & Lettenmaier, D. P. (2019). A curve-fitting method for
921 estimating bathymetry from water surface height and width. *Water Resources Research*,
922 55(5), 4288-4303.

923 Schramkowski, G. P., & De Swart, H. E. (2002). Morphodynamic equilibrium in straight tidal
924 channels: Combined effects of the Coriolis force and external overtides. *Journal of
925 Geophysical Research: Oceans*, 107(C12), 20-1-20-17.

926 Sun, W., Fan, J., Wang, G., Ishidaira, H., Bastola, S., Yu, J., Fu, Y.H., Kiem, A.S., Zuo, D., & Xu,
927 Z. (2018). Calibrating a hydrological model in a regional river of the Qinghai-Tibet plateau
928 using river water width determined from high spatial resolution satellite images. *Remote
929 Sensing of Environment*, 214, 100-114

930 Tarpanelli, A., Brocca, L., Lacava, T., Melone, F., Moramarco, T., Faruolo, M., & Tramutoli, V.
931 (2013). Toward the estimation of river discharge variations using MODIS data in ungauged
932 basins. *Remote Sensing of Environment*, **136**, 47-55.

933 Tarpanelli, A., Brocca, L., Barbetta, S., Faruolo, M., Lacava, T., & Moramarco, T. (2015). Coupling
934 MODIS and radar altimetry data for discharge estimation in poorly gauged river basins. *IEEE
935 Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, 8(1), 141-
936 148.

937 Tarpanelli, A., Santi, E., Tourian, M.J., Filippucci, P., Amarnath, G., & Brocca, L. (2019). Daily river
938 discharge estimates by merging satellite optical sensors and radar altimetry through artificial
939 neural network. *IEEE Transactions on Geoscience and Remote Sensing*, 57(1), 329-341.

940 Tarpanelli, A., Iodice, F., Brocca, L., Restano, M., & Benveniste, J. (2020). River flow monitoring
941 by Sentinel-3 OLCI and MODIS: Comparison and combination. *Remote Sensing*, 12, 3867.

942 Wu, T.X., Chen, R., Xu, Z.N., Yin, S.X., & Wang, S.D. 2024. Improving river medium-high flow
943 estimation by CM Hierarchical Classification (CMHC) method using Sentinel-2 imagery.
944 *Journal of Hydrology*, 633, 130990.

945 Xie, D., Gao, S., Wang, Z. B., Pan, C., Wu, X., & Wang, Q. (2017). Morphodynamic modeling of
946 a large inside sandbar and its dextral morphology in a convergent estuary: Qiantang
947 Estuary, China. *Journal of Geophysical Research: Earth Surface*, 122(8), 1553-1572.

948 Xu, H.-Q.: A Study on Information Extraction of Water Body with the Modified Normalized
949 Difference Water Index (MNDWI)[J], *Journal of Remote Sensing*, 5, 589-595, 2005.
950 Yang, H.Q., Mei, T.S., & Chen, X.Y. (2024). Variation of satellite-based suspended sediment concentration in the Ganges-Brahmaputra Estuary from 1990 to 2020. *Remote Sensing, 16*(2), 396.

951 Yu, G. A., Disse, M., Huang, H. Q., Yu, Y., & Li, Z. (2016). River network evolution and fluvial process responses to human activity in a hyper-arid environment-Case of the Tarim River in Northwest China. *Catena, 147*, 96-109.

956 Yu, Z.F., Zhang, J., Chen, Z.Y., Hu, Y.K., Shum, C.K., Ma, C.F., Song, Q.J., Yuan, X.H., Wang, B., & Zhou, B. (2023). Spatiotemporal evolutions of the suspended particulate matter in the Yellow River Estuary, Bohai Sea and characterized by Gaofen imagery. *Remote Sensing, 15*(19), 4769.

960 Zaji, A. H., Bonakdari, H., & Gharabaghi, B. (2018). Remote sensing satellite data preparation for simulating and forecasting river discharge. *IEEE Transactions on Geoscience and Remote Sensing, 56*(6), 3432-3441.

963 Zhao, C., Li, M., Wang, X., Liu, B., Pan, X., & Fang, H. (2022). Improving the accuracy of nonpoint-source pollution estimates in inland waters with coupled satellite-UAV data. *Water Research, 225*, 119208.

966 Zhao, C. S., Pan, T. L., Xia, J., Yang, S. T., Zhao, J., Gan, X. J., & Ding, S. Y. (2019). Streamflow calculation for medium-to-small rivers in data scarce inland areas. *Science of The Total Environment, 693*, 133571.

969 Zhao, C. S., Zhang, C. B., Yang, S. T., Liu, C. M., Xiang, H., Sun, Y., & Yu, Q. (2017). Calculating e-flow using UAV and ground monitoring. *Journal of Hydrology, 552*, 351-365.

971 Zhao, C., Pan, X., Yang, S., Liu, C., Chen, X., Zhang, H., & Pan, T. (2019). Measuring streamflow with low-altitude UAV imagery. *Acta Geographica Sinica, 74*(07), 1392-1408.

973 Zhang, R., & Sun, R. (2015). Evapotranspiration estimation and influence on water change in the Weihe River Basin, China. *Remote Sensing, 3*(4), 42-52.
976 **List of figure captions:**

977 **Figure 1** Study area and the representative gauging stations.

978 **Figure 2** C/M-V relationship diagrams before (a) and after (b) correction for the stations in the Wei River basin. Red lines and colored dash lines stand for trends of the entire and individual stations that determine the parameters in Eq. (2), respectively. Note that ground velocity measurements have been normalized due to the confidentiality of the data using $V_{\text{nor}} = (V - V_{\text{min}})/(V_{\text{max}} - V_{\text{min}})$, where $V_{\text{nor}}$ is the normalized velocity, $V$ is the original velocity, $V_{\text{min}}$ is the minimum velocity, and $V_{\text{max}}$ is the max velocity.

984 **Figure 3** Flowchart of the modified C/M method.

985 **Figure 4** Ideal cross-section or relationship between water depth and distance in Eq. 6 with $t \le 1$ (a) and $t \ge 1$ (b).

987 **Figure 5** Example of the slope breakpoint on the cross-section. Thick black curve segment: the segmented section above but closest to the water surface, which is taken as a basis to fit underwater cross-section. To avoid overfitting, as long as possible the thick black curve segment should be selected.

991 **Figure 6** Example of three scenarios when determining the below water cross-section. (a) Scenario 1 & (b) Scenario 2: no intersection is found in fitting the formula of two sides of the cross-section; (c) Scenario 3 intersection is found in fitting the formula of two sides of the cross-section.

995 **Figure 7** Extracted boundary for C value calculation from Landsat 5 TM. Sub-figures a & c: image No. 128; b & d: image No. 127. Images a & b are true color images, and c & d are near-infrared images. Red circles stand for urban areas selected with high reflectivity in the near-infrared band.

999 **Figure 8** The difference in reflectance when the river carries sediment (a, b) and when it is clear (c, d). The red line is the maximum river boundary (MRB) of the LT (a downward
1001 hydrological station in Figs. 1 & 6) based on TM images. The yellow square represents
1002 the selected permanent water pixels. Blue and green colors show the value of the
1003 Landsat surface reflectance.

1004 Figure 9 The relationship between the C/M and velocity in different seasons for WJB. The R² value is used to evaluate the correlation between the flow velocity and the C/M. The red and blue areas represent the 95% confidence interval that is calculated using the least squares method.

1008 Figure 10 The scatterplot between the estimated velocity and ground-measured counterpart. The
1009 grey dashed line represents the linear fit line between the ground-measured velocities
1010 (V<sub>G</sub>) and the estimated velocities (V<sub>E</sub>) at TS, LJC, and XY stations, where the average
1011 width of the river is less than 90 m. The blue dashed line represents the linear fit line
1012 between the V<sub>G</sub> and V<sub>E</sub> at LT and HX stations, where the average width of the river is
1013 larger than 90 m, and the black lines represent the 1:1 line. (a) flood season, (b) non-
1014 flood season. RW means the average river width. Note that the velocity has been
1015 normalized.

1016 Figure 11 Typical underwater cross-section estimations in 2007, 2008, and 2009 over LT(a-c),
1017 HX(d-f), XY (g-i), LJC (j-l), and TS (m-o).

1018 Figure 12 Validation of the discharge estimates against the ground-measured discharge. (a) flood
1019 season, (b) non-flood season, the grey dashed line represents the linear fit line of
1020 ground-measured discharge (D<sub>G</sub>) and estimated discharge (D<sub>E</sub>) at TS, LJC, and XY
1021 stations where the average width of the river is less than 90 m. The blue dashed line
1022 represents the linear fit line between the D<sub>G</sub> and D<sub>E</sub> at LT and HX stations, where the
1023 average width of the river is larger than 90 m and the black lines represent the 1:1 line,
1024 (c-g) discharge hydrograph of LT (c), HX (d), XY (e), LJC (f), and TS (g). Sparse
1025 discharge estimates are due to cloud contamination during high flow periods, coupled
1026 with the long revisiting time of the satellite. Only part of the discharge estimates were
1027 shown with concurrent in situ data during 2007-2009. Note that ground discharge
measurements have been normalized using $Q_{nor} = (Q - Q_{min})/(Q_{max} - Q_{min})$, where $Q_{nor}$ is the normalized discharge, $Q$ is the original discharge, $Q_{min}$ is the minimum discharge, and $Q_{max}$ is the max discharge.

Figure 13 The joint use of Landsat 8 and Sentinel-2 in discharge estimation at the LT station during 2019-2020. Note that the discharge has been normalized due to the confidentiality of the in situ data.
1035 **List of table captions:**

1036 Table 1 Satellite images used in this study

1037 Table 2 DJI Phantom-3-pro's basic parameters

1038 Table 3 Ground-measured data and their corresponding purposes for the six gauging stations
1039 (i.e., TS, LJC, WJB, XY, LT, and HX)

1040 Table 4 Fitting of *M* value relationship between ungauged river sections and the reference station
1041 (WJB).

1042 Table 5 The relative root-mean-square error (rRMSE) of the estimated flow velocity of all stations
1043 in the middle Yellow River.

1044 Table 6 Performance of discharge estimates at five sites in the Wei River basin (note that results
1045 of flood and non-flood seasons were aggregated due to limited data points).

1046 Table 7 Performance of discharge estimates using Landsat 8 and Sentinel-2 at LT in the Wei River
1047 basin during 2019-2020 (note that results of flood and non-flood seasons were aggregated
1048 due to limited data points).

1049 Table 8 Comparison of average relative errors of the modified *C/M* and original *C/M* in calculating
1050 the flow velocity.

1051 Table 9 rRMSE (relative root-mean-square error) of the flow velocity of all stations in the middle
1052 Yellow River calculated based on Sentinel-2.

1053

1054