Discharge Estimates for Ungauged Rivers Flowing over Complex High-Mountainous Regions based Solely on Remote Sensing-Derived Datasets

: Reliable information about river discharge plays a key role in sustainably managing water resources and better understanding of hydrological systems. Therefore, river discharge estimation using remote sensing techniques is an ongoing research goal, especially in small, headwater catchments which are mostly ungauged due to environmental or ﬁnancial limitations. Here, a novel method for river discharge estimation based entirely on remote sensing-derived parameters is presented. The model inputs include average river width, estimated from Landsat imagery by using the modiﬁed normalized di ﬀ erence water index (MNDWI) approach; average depth and velocity, based on empirical equations with inputs from remote sensing; channel slope from a high resolution shuttle radar topography mission digital elevation model (SRTM DEM); and channel roughness coe ﬃ cient via further analysis and classiﬁcation of Landsat images with support of previously published values. The discharge of the Lhasa River was then estimated based on these derived parameters and by using either the Manning equation (Model 1) or Bjerklie equation (Model 2). In general, both of the two models tend to overestimate discharge at moderate and high ﬂows, and underestimate discharge at low ﬂows. The overall performances of both models at the Lhasa gauge were satisfactory: comparisons with the observations yielded Nash–Sutcli ﬀ e e ﬃ ciency coe ﬃ cient (NSE) and R 2 values ≥ 0.886. Both models also performed well at the upper gauge (Tanggya) of the Lhasa River (NSE ≥ 0.950) indicating the transferability of the methodology to river cross-sections with di ﬀ erent morphologies, thus demonstrating the potential to quantify streamﬂow entirely from remote sensing data in poorly-gauged or ungauged rivers on the Tibetan Plateau. (via Model 1 and Model 2) discharge at Lhasa gauge over the study period of 1988–2018. Both models tend to overestimate discharge at moderate and high ﬂows and underestimate discharge at low ﬂows, but in general both models consistently capture the temporal changes. calibration validation, M.G.K. and L.W.; formal analysis, M.G.K.; investigation, M.G.K. L.W.; resources, L.W. and Z.H.; data curation, M.G.K., L.W. and Z.H.; writing—original draft preparation, M.G.K.; writing—review and editing, L.W., K.Y., D.C. and X.L.; visualization, M.G.K.; supervision, L.W.; project administration, L.W.; funding acquisition,


Introduction
Water is a critical yet limited natural resource required for socioeconomic development and for all life on the Earth. Rivers are the main components of natural water courses playing a crucial role in the

Study Area
The Lhasa River basin (LRB) is found in the central and southern part of the Tibet Autonomous Region of China along the northern ranges of Himalayas, at 29 • -31.3 • N and 90 • -93.4 • E ( Figure 1). The Lhasa River being the main tributaries of the Yarlung Zangbo River is 551 km long. Its altitude ranges from 3570 m above sea level (a.s.l.) on the valleys to 7136 m a.s.l. in its headwaters, where the mean altitude exceeds 3600 m and has an average channel gradient of 0.29 % [63,64]. According to Wu et al. [64], the impoundment and operation of Zhikong Dam (built in 2006 and located few kilometers above the Tanggya gauge station) directly changed the hydrological regime of the downstream braided sections of the Lhasa River and highly affects its morphology. The slope of Lhasa River decreased from 0.26 to 0.10% along the channel from the upstream mountain ridges to the downstream valleys. The average valley width ranges from 1.5 km to 3 km. The mean annual suspended load and bed load are 0.720 and 0.151 million tons, respectively [64]. We considered an area of about 26,225 km 2 in the Lhasa River above Lhasa gauge. The majority of the basin's annual runoff (about 90%) is received in summer. The basin includes seasonal freezing ground, and 670 km 2 of glaciers (2 % of the basin area) which contribute significantly to the total runoff, but mostly during summer [65,66]. The mean annual discharge of the Lhasa River is 0.6 trillion m 3 (nearly 1/8th of Yellow River flow, China) [67], and is mostly sourced from groundwater, melt-water and precipitation [68].
All the surroundings of the LRB urban and rural areas community development activities are based on water resources of the Lhasa River. There are studies conducted in the LRB but mostly focused on climate change and extreme events monitoring [65,69], while a research approach based on space borne datasets considering the basin's in-situ data gap is still missing. Therefore, a study on the analysis and quantification of River discharge for this kind of widely exploited River using satellite technologies is timely to understand its dynamic change. Hence, by doing so, the present study aims to address this gap and to lay a foundation for the basin's water resources planning and management.

Data
SRTM DEM and Landsat 5, 7 and 8 satellite imageries with 30 m spatial resolution were acquired (https://earthexplorer.usgs.gov) and processed to estimate the channel slope and the reach-average effective width during the study period (1988-2018), respectively. We had manually inspected the quality of downloaded Landsat imageries, and those badly affected by cloud or snow were discarded from our image processing and analysis. We also used visual interpretation of the Landsat images,

Data
SRTM DEM and Landsat 5, 7 and 8 satellite imageries with 30 m spatial resolution were acquired (https://earthexplorer.usgs.gov) and processed to estimate the channel slope and the reach-average effective width during the study period (1988-2018), respectively. We had manually inspected the quality of downloaded Landsat imageries, and those badly affected by cloud or snow were discarded from our image processing and analysis. We also used visual interpretation of the Landsat images, in combination with previous studies, to identify the channel morphological conditions needed to assign the channel roughness coefficient.
Finally, observed discharge data at Lhasa gauge  and Tanggya gauge (1999-2018) on the Lhasa River were collected from the China Ministry of Water Resources (MWR). These data were used to validate the discharge estimated using the two models/equations. Generally, the Lhasa gauge has longer and better in-situ measured time series discharge data than the Tanggya gauge.
Remote Sens. 2020, 12, 1064 5 of 22 In fact, in our study particularly at Tanggya gauge site of Lhasa River, the number of Landsat images analyzed from 1999-2013 are much lower than the years from 2014-2018. Since the observed flow data at Tanggya gauge site are also limited and mainly available during the flooding season days, there is a mismatch between the satellite overpass date and in-situ measurement date which forces us to analyze the Landsat images only in those dates where in-situ discharge measurements are available. Moreover, particularly during the years of 2000, 2003 and 2004, we had the in-situ gauge data only in some limited monsoon flooding days and unfortunately we didn't get Landsat images matching with those in-situ measurement days to estimate Lhasa river discharge and due to this data gap we had excluded it from our analysis of discharge estimation.

Methodology
In this section, the methodology employed to extract the input surface hydraulic parameters required for the Model 1 and Model 2 discharge estimates are discussed in detail.

Estimating Water Surface Area (WSA) and Effective Width (W)
The resolution of the satellite images and heterogeneity of the land cover are essential considerations when classifying and identifying different features, including water bodies. The MNDWI has been tested and widely used during the last few decades to delineate surface water features, including rivers, using remotely sensed images at various scales [10,[49][50][51][52]56,70].
The accuracy of the MNDWI in distinguishing and delineating water features from other non-water features (all negative values) is reasonably high when compared to the other water body detection indices, particularly for rivers flowing over the TP mountainous regions [10,18,49,50,58,59]. The MNDWI can also be used to detect mixed pixels of small water bodies [50]. By using combined bands and manually adjusting threshold values based on the spectral signature and proportions of water and non-water features, the MNWDI can achieve a more accurate surface water retrieval [56]. The MNDWI is a more stable and reliable index than other indices to effectively distinguish water pixels from the other background features since it uses the shortwave infrared (SWIR) band, which can effectively represent the typical characteristics of water. The MNDWI can enhance and extract open water pixels by effectively suppressing the background noise from mixed pixels leading to bimodal (water and non-water) image histogram [71]. The spectral properties of water can be represented by the distribution of MNDWI values that varies with location and season. Based on our analysis of Landsat images over the LRB, MNDWI values range from −1.0 to 1.0. To overcome problems related to the misclassification of mixed pixels, we manually adjusted the MNDWI threshold value for each image processed over the study period in which its value ranges from −0.05 to 0.35 (the threshold value was lower in the dry season and higher in the wet monsoon season). The optimum threshold value was then chosen as 0.2 by trial and error. Terrain analysis via the SRTM DEM was used to avoid commission error caused by mountain shadows. Hence, the MNDWI is used to extract the water body by manually adjusting the optimum threshold values between water and non-water classes from the bimodal histogram inter-class variance following Otsu [72] (Figure 2): where R Green is the surface-reflectance of Band 2 for Landsat 5 and 7 images (Band 3 for Landsat 8 images) and R MIR is the surface-reflectance of Band 5 for Landsat 5 and 7 images (Band 6 for Landsat 8 images).
Several studies have shown that the in-bank river width can be successfully estimated using remote sensing images [4,42,73,74]. Here, we used Arc Map to digitize the WSA using Landsat images based on the MNDWI water body extraction approach in which the River width was calculated as the ratio of WSA and its appropriate river reach length (L) (Equation (2)) by selecting relatively stable and non-braided cross-sections to avoid errors in width measurements of braided river sections associated Remote Sens. 2020, 12, 1064 6 of 22 with contamination of mixed pixels [11,35]; this is the preferred method of measuring width for such types of river cross sections [10,30,36].
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 23 Several studies have shown that the in-bank river width can be successfully estimated using remote sensing images [4,42,73,74]. Here, we used Arc Map to digitize the WSA using Landsat images based on the MNDWI water body extraction approach in which the River width was calculated as the ratio of WSA and its appropriate river reach length (L) (Equation (2)) by selecting relatively stable and non-braided cross-sections to avoid errors in width measurements of braided river sections associated with contamination of mixed pixels [11,35]; this is the preferred method of measuring width for such types of river cross sections [10,30,36].

Estimating Channel Slope (S)
A viable channel slope estimation using the SRTM DEM can be derived by considering an appropriate reach length [11,62,75]. Using this approach, we computed the Lhasa River slope from the SRTM DEM by extending the reach lengths to a distance sufficient to accommodate height variations and minimize the associated errors. Specifically, the altitude of the river's centerline was obtained from the SRTM DEM and hence the slope was estimated in meters per kilometer. In order to do this, a total Lhasa River reach length of more than 250 km from upstream of the Zhikong dam to its downstream end was considered, and the resulting mean gradient was 1.98 m km -1 (average channel slope of 0.002) with a standard deviation of 3.67 m km -1 .

Estimating Channel Slope (S)
A viable channel slope estimation using the SRTM DEM can be derived by considering an appropriate reach length [11,62,75]. Using this approach, we computed the Lhasa River slope from the SRTM DEM by extending the reach lengths to a distance sufficient to accommodate height variations and minimize the associated errors. Specifically, the altitude of the river's centerline was obtained from the SRTM DEM and hence the slope was estimated in meters per kilometer. In order to do this, a total Lhasa River reach length of more than 250 km from upstream of the Zhikong dam to its downstream end was considered, and the resulting mean gradient was 1.98 m km −1 (average channel slope of 0.002) with a standard deviation of 3.67 m km −1 .

Estimating Channel Roughness Coefficient (n)
The channel roughness or resistance is strongly affected by vegetation cover, river bank and bed composition, obstruction, channel sinuosity (degree of meander) and channel irregularity [76]. Roughness can be estimated using several approaches based on visual interpretation and the assignment of tabulated values of Chow [77]. Ground surveys, in particular for narrow channels [78], and the use of empirical equations/formulas [79] are well-established methods. In this study, we used visual interpretation by inspecting Landsat images along the LRB. The Lhasa River bed, specifically the 150 km Remote Sens. 2020, 12, 1064 7 of 22 downstream reach is composed of gravel [64]. Sichangi et al. [62] retrieved the roughness coefficient from Landsat images when estimating discharge of the Yangtze River, China. Their methodology was applied in the present study, but using the Lhasa River's morphological characteristics. Classification of Landsat images and band combination approach was used to identify the prevailing channel conditions and other features, from which the required table values were derived as set out and discussed below.
Generally, the morphology of middle reach of Lhasa River is more stable than the lower reach [64]. The banks of Lhasa River around the Tanggya gauge station reach are stable, the channel is nearly straight without sharp bends, the size and shape of the river cross-sections are not changing significantly, and the degree of meandering (i.e., the ratio of the total length of the meandering channel in the reach being considered to the straight length of the channel reach) is 1.05. While, the Lhasa gauge station reach is not stable with changing crossing sections and higher meandering ratio of 1.33, and this reach being located in the downstream end of Lhasa River has three tributaries transporting sediments that can affect the channel roughness and also the dynamics of the river flow. According to Chow [77], meanders can increase the n values by as much as 30 percent where flow is confined within a stream channel.
The lookup table values to compute n based on the channel morphological conditions including types of channel materials involved, degree of irregularity, variation of channel cross-section, relative effect of obstruction, vegetation and degree of meandering can be found from previous publications [18,62,77]. These channel morphological conditions can be also analyzed using satellite imageries including Landsat and with an integration of published literatures [18,62]. Therefore, we calculated n based on detailed classification and inspection of Landsat images around the gauging sites following the approaches of Chow's [77] channel roughness lookup table values via Equation (3). Sichangi et al. [62] and Genanu et al. [18] had also used this method to successfully derive a reasonable channel roughness values for different Global Rivers and the same approach was also applied here. First, the basic channel condition (n 0 ) is assigned a value of 0.025 [80]. Other factors, which are included additively, are: a minor to moderate degree of irregularity (n 1 = 0.006); occasionally alternating river cross-section (n 2 = 0.005); minor obstructions which affect the channel flow pattern (n 3 = 0.006); low vegetation (n 4 = 0.004); and appreciable channel meandering with a ratio of 1.33 between channel length and straight-line distance of a reach (n 5 = 1.15). These factors were combined in Equation (3), yielding n = 0.053 for the Lhasa River around the Lhasa gauging station. Similarly, the value of n at Tanggya gauging station was estimated to be 0.043 as discussed in Genanu et al. [18]: n = (n 0 + n 1 + n 2 + n 3 + n 4 )n 5 (3)

Estimating River Depth and Velocity
Bjerklie [81] developed a general regression equation to estimate bank-full depth as a function of observed width and slope, in the absence of direct depth measurements, but noted the very high standard error. Various remote sensing methods have been developed and applied to estimate river depth [82]. Water level or discharge can be monitored via satellite altimetry, mainly applied to wide rivers (e.g., width > 800 m) over relatively longer periods (> 1 week) due to its low spatio-temporal resolution [11,83]. Satellite altimetry based water level or discharge estimated errors in a narrow high mountainous rivers are very high as compared to the wider rivers [10], and very few attempts have been made to apply this technique to establish a database of small to medium-sized rivers and lakes [83].
Likewise, there is lack of altimetry data to estimate water level or depth for small-to medium-size rivers such as the Lhasa River. Therefore, we have estimated the depth of the river empirically as a function of hydraulic variables derived from remote sensing datasets (Equation (4)). The method is easy and effective for estimating average depths at the sub-reach scale.
The well-known Manning velocity equation [84] can be simplified and equated to the traditional discharge formula (i.e., discharge equal to the product of mean cross-sectional flow width (W), depth Due to a lack of either direct in-situ or indirect depth observations from satellite altimetry, we estimated the reach-average depth of the Lhasa River as a function of average velocity, channel roughness coefficient and slope as indicated by Equation (4) and velocity was estimated using Equation (5).
Bjerklie et al. [35] derived a statistically-based discharge formula, from which Tourian et al. [85] developed a model that can estimate average velocity (V) with reasonable accuracy as a function of average river width and slope as shown in Equation (5).

Estimating River Discharge
Direct measurement of river discharge is very difficult; therefore, discharge is indirectly measured using surface hydraulic flow parameters which are themselves derived from either in-situ measurements or from remote sensing datasets via mathematical equations or calibrated relationships (rating curves). River discharge estimates could potentially be achieved via remote sensing techniques, given accurate estimates of the input surface hydraulic variables like depth, width, velocity, slope and channel resistance [1,11,15,62,75,82] are available. A simplified hydraulic equation (Equation (6)) that can estimate river discharge based on Manning's velocity equation for relatively wide, open channel flow (i.e., when the width/depth ratio is at least 10) has been applied in many rivers globally [10,11,75]. The width to depth ratio of the Lhasa River is consistent with this approach, so Equation (6) could be used to estimate the discharge.
Discharge estimates based on statistical multi-variate equations with reasonably computed input parameters are showing great potential at present. Bjerklie et al. [15] derived a hydraulic formula (Equation (7)) to compute discharge of a wide range of rivers and the same approach using remote sensing derived parameters was also used here.
In summary, we extracted the slope from a 30 m resolution SRTM DEM, width from Landsat images, used reasonable estimates of channel resistance coefficient, and estimated velocity and depth using empirical equations. These values were substituted into Equation (6), herein after Model 1, and into Equation (7), herein after Model 2, to estimate discharge of the Lhasa River reach around the Lhasa gauge site ( Figure 3). Finally, the performances of discharge estimates based on these two models were compared and evaluated using several statistical metrics.
Generally, this study was applied in non-braided sections of the Lhasa River and all the methodologies may not be directly used in another multithread braided rivers due to the dynamic variability of the channel geometry and morphology, and it needs further validation over different hydro-climatological channel environments. Optical images including Landsat are also contaminated by cloud and vegetation cover along the river banks and this will result in misclassification of water pixels affecting the river width estimation. Currently, there is also a tradeoff between spatial and temporal resolution of satellites and due to these limitations Landsat satellite cannot continuously monitor important hydrological events such as peak floods during the monsoon season.
variability of the channel geometry and morphology, and it needs further validation over different hydro-climatological channel environments. Optical images including Landsat are also contaminated by cloud and vegetation cover along the river banks and this will result in misclassification of water pixels affecting the river width estimation. Currently, there is also a tradeoff between spatial and temporal resolution of satellites and due to these limitations Landsat satellite cannot continuously monitor important hydrological events such as peak floods during the monsoon season.

Evaluation of Models Performance
The Lhasa River discharge estimates from the two models were evaluated based on the following performance metrics.

Evaluation of Models Performance
The Lhasa River discharge estimates from the two models were evaluated based on the following performance metrics. Relative where Q o is the observed discharge, Q e is the estimated discharge, Q o is the average observed discharge and k is the total number of observations.

Water Surface Area (WSA) and Width
Satellite datasets provide frequent spatiotemporal coverage of many features on the Earth, and allow, for example, surface water dynamics at various scales to be monitored. The WSA of the Lhasa River reach around Lhasa gauge was mapped using the MNDWI approach from Landsat 5, 7 and 8 images acquired during the study period . The non-braided river reach was selected, and then the reach-average effective width was computed as ratio of the WSA and the corresponding reach length.
The number of available Landsat images during the high flow monsoon season is much lower than that during the moderate and low flow seasons, resulting in denser sampling during the dry periods. The results of time series plots of both WSA and width show that Landsat images have the potential to capture seasonal variations in surface water ( Figure 4). The estimated WSA and width were higher during the wet flood season and lower during the dry season. The maximum and minimum effective widths of the Lhasa River reach around Lhasa Gauge site for all Landsat images analyzed in the study period were 198 m and 32 m, respectively ( Figure 4). Relative where Q o is the observed discharge, Q e is the estimated discharge, Q o is the average observed discharge and k is the total number of observations.

Water Surface Area (WSA) and Width
Satellite datasets provide frequent spatiotemporal coverage of many features on the Earth, and allow, for example, surface water dynamics at various scales to be monitored. The WSA of the Lhasa River reach around Lhasa gauge was mapped using the MNDWI approach from Landsat 5, 7 and 8 images acquired during the study period . The non-braided river reach was selected, and then the reach-average effective width was computed as ratio of the WSA and the corresponding reach length.
The number of available Landsat images during the high flow monsoon season is much lower than that during the moderate and low flow seasons, resulting in denser sampling during the dry periods. The results of time series plots of both WSA and width show that Landsat images have the potential to capture seasonal variations in surface water ( Figure 4). The estimated WSA and width were higher during the wet flood season and lower during the dry season. The maximum and minimum effective widths of the Lhasa River reach around Lhasa Gauge site for all Landsat images analyzed in the study period were 198 m and 32 m, respectively ( Figure 4). Generally, the WSA and effective width extracted from Landsat images around Lhasa gauge site show a good positive linear relation as best fit of the scatter plot ( Figure 5) and also some plot points are seen deviating from the trend line due to the outliers. The high correlation (R 2 = 0.765) between these two parameters demonstrates their consistency and also, to some extent, that dynamic changes in characteristics and trends depend on the hydrological wet and dry seasons. Higher values of both WSA and effective width were attained during the summer monsoon seasons than during the winter Generally, the WSA and effective width extracted from Landsat images around Lhasa gauge site show a good positive linear relation as best fit of the scatter plot ( Figure 5) and also some plot points are seen deviating from the trend line due to the outliers. The high correlation (R 2 = 0.765) between these two parameters demonstrates their consistency and also, to some extent, that dynamic changes in characteristics and trends depend on the hydrological wet and dry seasons. Higher values of both WSA and effective width were attained during the summer monsoon seasons than during the winter dry seasons. To minimize the errors associated with estimates of WSA and effective width, we selected a non-braided reach with stable banks around the Lhasa gauge site.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 23 dry seasons. To minimize the errors associated with estimates of WSA and effective width, we selected a non-braided reach with stable banks around the Lhasa gauge site.

River Depth and Velocity
Reasonable estimates of river depth and velocity from remotely sensed hydraulic parameters have been reported by many studies [10,11,35,62,81,83,85]. In the present study, due to the lack of insitu and remote-sensing datasets with which to derive either depth or velocity of the Lhasa River reach around the Lhasa gauge location, we empirically estimated the reach average depth and velocity. This approach yielded reasonable accuracy, and supplements estimates of river discharge in ungauged or poorly-gauged river basins like Lhasa River basin. The time series of empiricallyderived reach-averaged depth and velocity, based on remote sensing datasets, demonstrate acceptable accuracy in the estimates of both parameters. Generally, temporal patterns of both parameters follow the hydrological wet and dry seasons ( Figure 6). Increasing flow during the wet season results in rising water level and faster velocity, and vice-versa. The long-term average depth and velocity of the Lhasa river reach around the Lhasa gauge station were estimated as 1.44 m and 1.07 m/s, respectively.

River Depth and Velocity
Reasonable estimates of river depth and velocity from remotely sensed hydraulic parameters have been reported by many studies [10,11,35,62,81,83,85]. In the present study, due to the lack of in-situ and remote-sensing datasets with which to derive either depth or velocity of the Lhasa River reach around the Lhasa gauge location, we empirically estimated the reach average depth and velocity. This approach yielded reasonable accuracy, and supplements estimates of river discharge in ungauged or poorly-gauged river basins like Lhasa River basin. The time series of empirically-derived reach-averaged depth and velocity, based on remote sensing datasets, demonstrate acceptable accuracy in the estimates of both parameters. Generally, temporal patterns of both parameters follow the hydrological wet and dry seasons ( Figure 6). Increasing flow during the wet season results in rising water level and faster velocity, and vice-versa. The long-term average depth and velocity of the Lhasa river reach around the Lhasa gauge station were estimated as 1.44 m and 1.07 m/s, respectively.
The time series mean velocity was estimated from the reach-averaged width derived from Landsat imageries, the mean slope derived SRTM DEM, and the average depth was estimated as a function of velocity, slope and channel roughness in the corresponding empirical equations. The scatter plot clearly shows the non-linear relationship of the average depth, velocity and effective width. The hydraulic relationship between average depth, velocity and effective width shows a power function relationship with strong correlation (R 2 = 0.998), which can potentially be used to predict the unknown parameter with reasonable accuracy if the others are known (Figure 7). The time series mean velocity was estimated from the reach-averaged width derived from Landsat imageries, the mean slope derived SRTM DEM, and the average depth was estimated as a function of velocity, slope and channel roughness in the corresponding empirical equations. The scatter plot clearly shows the non-linear relationship of the average depth, velocity and effective width. The hydraulic relationship between average depth, velocity and effective width shows a power function relationship with strong correlation (R 2 = 0.998), which can potentially be used to predict the unknown parameter with reasonable accuracy if the others are known (Figure 7). The power function trend lines between depth and velocity, and between depth and effective width, show two distinct relationships (one from low to moderate values and another one from  The time series mean velocity was estimated from the reach-averaged width derived from Landsat imageries, the mean slope derived SRTM DEM, and the average depth was estimated as a function of velocity, slope and channel roughness in the corresponding empirical equations. The scatter plot clearly shows the non-linear relationship of the average depth, velocity and effective width. The hydraulic relationship between average depth, velocity and effective width shows a power function relationship with strong correlation (R 2 = 0.998), which can potentially be used to predict the unknown parameter with reasonable accuracy if the others are known (Figure 7).  (4)) versus average velocity (estimated using Equation (5)); (b) scatterplots of average depth versus effective width (estimated using Equation (2)) at Lhasa gauging station during the study period 1988-2018. Both plots show strong relationships between depth and average effective width or velocity, with two distinct curves for lower values and moderate to high values.
The power function trend lines between depth and velocity, and between depth and effective width, show two distinct relationships (one from low to moderate values and another one from  (4)) versus average velocity (estimated using Equation (5)); (b) scatterplots of average depth versus effective width (estimated using Equation (2)) at Lhasa gauging station during the study period 1988-2018. Both plots show strong relationships between depth and average effective width or velocity, with two distinct curves for lower values and moderate to high values.
The power function trend lines between depth and velocity, and between depth and effective width, show two distinct relationships (one from low to moderate values and another one from moderate to high values), with the trend lines breaking at 2.4 m and 1.6 m/s, or 2.4 m and 110 m, respectively. The accuracy of our average depth estimate was evaluated though statistical correlation with the observed stage data, yielding a strong linear relationship between the observed stage and estimated average depth (R 2 = 0.676) (Figure 8). Generally, as compared with the observed stage, the estimated depth using Equation (4) shows underestimation trend except its overestimation during a historical peak discharge event of 1706 m 3 s −1 at Lhasa gauge recorded on September 1, 2018 in which, the in-situ stage reading was 4.0 m while the estimated depth was 4.6 m (Figure 8).
respectively. The accuracy of our average depth estimate was evaluated though statistical correlation with the observed stage data, yielding a strong linear relationship between the observed stage and estimated average depth (R 2 = 0.676) (Figure 8). Generally, as compared with the observed stage, the estimated depth using Equation (4) shows underestimation trend except its overestimation during a historical peak discharge event of 1706 m 3 s -1 at Lhasa gauge recorded on September 1, 2018 in which, the in-situ stage reading was 4.0 m while the estimated depth was 4.6 m (Figure 8).

River Discharge
All the datasets required to estimate river discharge can be derived from space-borne platforms. After deriving the input surface hydraulic parameters entirely from high-resolution satellite datasets, we used Equation (6) (Model 1) and Equation (7) (Model 2) to estimate discharge of the Lhasa river in a non-braided reach around the Lhasa gauge site. Results show reasonable accuracy, demonstrating the potential of remote sensing-based discharge estimates to complement existing flow data and to fill gaps, especially in poorly-gauged or ungauged basins. Model 1 discharge estimates were always lower than Model 2 estimates (Figure 9) because Model 1 includes the channel resistance coefficient (indirectly related to discharge, as evident in Equation (6)), whereas Model 2 does not include it. However, the channel resistance coefficient was also used as an input in the calculation of river depth, which was a parameter required to estimate river discharge in both models.
Temporal variations in discharge estimated via Model 1 and Model 2 match up with the observed discharge, capturing the low and peak flows during the dry and wet seasons, respectively. In general, Model 1 and Model 2 overestimate discharge during moderate and high flows but underestimate discharge during low flows (Figure 9). Meanwhile, Model 2 discharge estimate was generally higher than that of Model 1, especially during high flows. Generally, the tendencies to overestimate discharge in Model 1 and Model 2 were slightly enhanced during the high flow seasons of discharge higher than 1000 m 3 s -1 . Perhaps, depth estimated based on roughness coefficient, channel slope and velocity (all estimated values); river width extracted from Landsat images that has not been validated by observed values; the assumption of temporally constant channel slope and especially roughness coefficient plays an important role in overestimation of discharge. Moreover, Figure 8. Observed stage and estimated average depth (using Equation (4)) at Lhasa gauging station. The strong correlation demonstrates the accuracy of our average depth estimation methodology.

River Discharge
All the datasets required to estimate river discharge can be derived from space-borne platforms. After deriving the input surface hydraulic parameters entirely from high-resolution satellite datasets, we used Equation (6) (Model 1) and Equation (7) (Model 2) to estimate discharge of the Lhasa river in a non-braided reach around the Lhasa gauge site. Results show reasonable accuracy, demonstrating the potential of remote sensing-based discharge estimates to complement existing flow data and to fill gaps, especially in poorly-gauged or ungauged basins. Model 1 discharge estimates were always lower than Model 2 estimates (Figure 9) because Model 1 includes the channel resistance coefficient (indirectly related to discharge, as evident in Equation (6)), whereas Model 2 does not include it. However, the channel resistance coefficient was also used as an input in the calculation of river depth, which was a parameter required to estimate river discharge in both models.
Temporal variations in discharge estimated via Model 1 and Model 2 match up with the observed discharge, capturing the low and peak flows during the dry and wet seasons, respectively. In general, Model 1 and Model 2 overestimate discharge during moderate and high flows but underestimate discharge during low flows (Figure 9). Meanwhile, Model 2 discharge estimate was generally higher than that of Model 1, especially during high flows. Generally, the tendencies to overestimate discharge in Model 1 and Model 2 were slightly enhanced during the high flow seasons of discharge higher than 1000 m 3 s −1 . Perhaps, depth estimated based on roughness coefficient, channel slope and velocity (all estimated values); river width extracted from Landsat images that has not been validated by observed values; the assumption of temporally constant channel slope and especially roughness coefficient plays an important role in overestimation of discharge. Moreover, errors in the discharge measurements, and other controlling factors, are also additional sources of uncertainties in discharge estimates.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 23 errors in the discharge measurements, and other controlling factors, are also additional sources of uncertainties in discharge estimates.

Further Validation of the Methodology
Having obtained promising results for the Lhasa River as described above, and to further validate the transferability of our methodology to different river morphologies, we applied this methodology to the Tanggya gauge in the Lhasa River (located about 85 km upstream of the Lhasa gauge). We used the same methodologies described in Sections 2.3.1-2.3.5 to derive the input surface hydraulic variables for Model 1 and Model 2 discharge estimations, then we also followed the same approaches mentioned in Section 2.3.6 to evaluate the two models performances. The time series plots of the estimated model input parameters including width, depth and velocity at Tanggya gauge are also presented ( Figure 10) demonstrating the temporal variation of these parameters. We did not compare the estimated width and velocity data with the corresponding in-situ measurements due to lack of observed data, but for depth we had compared the estimated depth with an in-situ stage data ( Figure 11). Unlike the Lhasa gauge, Tanggya gauge has very limited in-situ discharge or stage observations, requiring us to process and analyze only a limited number of Landsat images matching those days with in-situ observations for validation of the two models. Analysis of all the Landsat images for the study period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)

Further Validation of the Methodology
Having obtained promising results for the Lhasa River as described above, and to further validate the transferability of our methodology to different river morphologies, we applied this methodology to the Tanggya gauge in the Lhasa River (located about 85 km upstream of the Lhasa gauge). We used the same methodologies described in Section 2.3.1-Section 2.3.5 to derive the input surface hydraulic variables for Model 1 and Model 2 discharge estimations, then we also followed the same approaches mentioned in Section 2.3.6 to evaluate the two models performances. The time series plots of the estimated model input parameters including width, depth and velocity at Tanggya gauge are also presented ( Figure 10) demonstrating the temporal variation of these parameters. We did not compare the estimated width and velocity data with the corresponding in-situ measurements due to lack of observed data, but for depth we had compared the estimated depth with an in-situ stage data ( Figure 11). Unlike the Lhasa gauge, Tanggya gauge has very limited in-situ discharge or stage observations, requiring us to process and analyze only a limited number of Landsat images matching those days with in-situ observations for validation of the two models. Analysis of all the Landsat images for the study period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)  The estimated average depth and in-situ gauge stage data at Tanggya gauge shows a very strong correlation (R 2 = 0.992), providing strong support for the accuracy of our average depth estimation methodology ( Figure 11).
Both Model 1 and Model 2 discharge estimates at Tanggya gauge tended to underestimate discharge at low flows and over-estimate discharge at high flows, consistent with the trend for Lhasa gauge ( Figure 12). However, in contrast to the Lhasa gauge, discharge estimated by Model 1 at Tanggya gauge was slightly higher than that estimated by Model 2, because the river bank of the reach at Tanggya gauge is relatively stable and the estimated channel roughness coefficient at this gauge (n = 0.043) was lower than that of the Lhasa gauge reach. Evaluation of the Model 1 and The estimated average depth and in-situ gauge stage data at Tanggya gauge shows a very strong correlation (R 2 = 0.992), providing strong support for the accuracy of our average depth estimation methodology ( Figure 11). Figure 11. The scatterplots of observed stage and estimated average depth (using Equation (4)) at Tanggya gauging station for further validation of the methodology shows a very good correlation demonstrating the accuracy of our average depth estimation methodology especially in data scarce river basins.
Both Model 1 and Model 2 discharge estimates at Tanggya gauge tended to underestimate discharge at low flows and over-estimate discharge at high flows, consistent with the trend for Lhasa gauge ( Figure 12). However, in contrast to the Lhasa gauge, discharge estimated by Model 1 at Tanggya gauge was slightly higher than that estimated by Model 2, because the river bank of the reach at Tanggya gauge is relatively stable and the estimated channel roughness coefficient at this The estimated average depth and in-situ gauge stage data at Tanggya gauge shows a very strong correlation (R 2 = 0.992), providing strong support for the accuracy of our average depth estimation methodology ( Figure 11). Figure 11. The scatterplots of observed stage and estimated average depth (using Equation (4)) at Tanggya gauging station for further validation of the methodology shows a very good correlation demonstrating the accuracy of our average depth estimation methodology especially in data scarce river basins.
Both Model 1 and Model 2 discharge estimates at Tanggya gauge tended to underestimate discharge at low flows and over-estimate discharge at high flows, consistent with the trend for Lhasa gauge ( Figure 12). However, in contrast to the Lhasa gauge, discharge estimated by Model 1 at Tanggya gauge was slightly higher than that estimated by Model 2, because the river bank of the reach at Tanggya gauge is relatively stable and the estimated channel roughness coefficient at this Figure 11. The scatterplots of observed stage and estimated average depth (using Equation (4)) at Tanggya gauging station for further validation of the methodology shows a very good correlation demonstrating the accuracy of our average depth estimation methodology especially in data scarce river basins.
2 discharge estimates by comparison with in-situ discharge data (1999-2018) generally showed strong correlation, with NSEs of 0.97 and 0.95, RRMSEs of 16.5 % and 23.8 %, MBEs of -25.9 m 3 s -1 and -59.7 m 3 s -1 , RMSEs of 48.8 m 3 s -1 and 70.3 m 3 s -1 and REs of -8.8% and -20.2% for Model 1 and Model 2, respectively. Therefore, the additional validation of Model 1 and Model 2 at Tanggya gauge on the Lhasa River is consistent with the promising results presented for the Lhasa gauge and further supports the transferability of the proposed methodology to rivers with different morphological characteristics.

Discussions
The potential of river discharge estimation entirely from remotely sensed hydraulic variables is clearly verified in this study and the applicability of the methodology for ungauged and poorly gauged basins with different channel morphological settings is also promising. By combining and equating all the input parameters derived from satellite datasets via rigorous analysis and computations, i.e., W estimated using Equation (2), n using Equation (3), D using Equation (4), V using Equation (5) and S from SRTM DEM, the Lhasa River discharge was estimated using Equation (6) or Model 1 and Equation (7) or Model 2 and finally the models performance was evaluated via various metrics yielding in a good agreement with the gauge's discharge observation (Figure 9 and Figure 12). The effective width of Lhasa River was extracted from Landsat imageries with a 30 m spatial resolution via MNDWI approach and currently, there is always a tradeoff between satellite sensors spatial and temporal resolutions. Landsat's temporal resolution (16 days) has a short come to continously monitor important hydrological events such as peak floods and also it is contaminated by cloud cover. Therefore, an integrative utilization of multiple satellite products with better spatial and temporal resolution such as CubeSat, Sentinel, the upcoming Surface Water and Ocean Topography (SWOT) and others by accounting the size of the targeted water body can greatly help to tackle this problem. Additionally, the utilization of panchromatic bands can also improve results as they have very high signal compared to multispectral bands.
River discharge estimated from multiple hydraulic parameters (width, depth and slope) is considered to achieve higher accuracy than using a single parameter [10,11,15], as supported by our

Discussions
The potential of river discharge estimation entirely from remotely sensed hydraulic variables is clearly verified in this study and the applicability of the methodology for ungauged and poorly gauged basins with different channel morphological settings is also promising. By combining and equating all the input parameters derived from satellite datasets via rigorous analysis and computations, i.e., W estimated using Equation (2), n using Equation (3), D using Equation (4), V using Equation (5) and S from SRTM DEM, the Lhasa River discharge was estimated using Equation (6) or Model 1 and Equation (7) or Model 2 and finally the models performance was evaluated via various metrics yielding in a good agreement with the gauge's discharge observation (Figures 9 and 12). The effective width of Lhasa River was extracted from Landsat imageries with a 30 m spatial resolution via MNDWI approach and currently, there is always a tradeoff between satellite sensors spatial and temporal resolutions. Landsat's temporal resolution (16 days) has a short come to continously monitor important hydrological events such as peak floods and also it is contaminated by cloud cover. Therefore, an integrative utilization of multiple satellite products with better spatial and temporal resolution such as CubeSat, Sentinel, the upcoming Surface Water and Ocean Topography (SWOT) and others by accounting the size of the targeted water body can greatly help to tackle this problem. Additionally, the utilization of panchromatic bands can also improve results as they have very high signal compared to multispectral bands.
River discharge estimated from multiple hydraulic parameters (width, depth and slope) is considered to achieve higher accuracy than using a single parameter [10,11,15], as supported by our results. Our approach is unique in deriving average depth and velocity from simplified empirical equations, yet with reasonable accuracy; these parameters can then be used to estimate river discharge, especially in poorly gauged and ungauged basins.  4.2% and 5.3%, respectively. These results demonstrate the reasonable accuracy and great potential of reliably estimating river discharge entirely from remote sensing-derived surface hydraulic parameters, with either Model 1 or Model 2. This is encouraging given the lack of observed flow data in remote and inaccessible river basins. The errors in Model 1 and Model 2 discharge estimates are largely associated with the uncertainties in model input parameters estimates (including river width, depth, slope and channel resistance coefficient). Errors in the discharge measurements, and other controlling factors, are also additional sources of uncertainties in discharge estimates [86][87][88][89][90]. In general, the uncertainties in the two models tested here are within the ranges reported in other studies.
Although the Model 1 and Model 2 discharge estimates achieve varying accuracy in different rivers [91], both models performed well in rivers with discharge less than 20,000 m 3 s −1 , in which they showed strong correlation with relatively high NSE and R 2 values [11,92]. Evaluation of Model 1 for the Yangtze River at two virtual stations yielded NSEs of 0.50 and 0.76, and RRMSEs of 26.34 and 34.63% [62]. According to Sichangi et al. [11], the NSE and RRMSE were in the ranges 0. Moreover, the transferability of the proposed methodology was also tested in another section of Lhasa River (i.e., at Tanggya gauge, which is located about 85 km far from the Lhasa gauge site in the upstream direction) with different channel morphological characteristics and the results of the two models discharge estimation were also reliable and promising with reasonable accuracy indicating the potential capability of the methodology to estimate river discharge with input parameters entirely derived from remote sensing datasets which can potentially complement the current gaps and challenges of streamflow data especially over poorly gauged and ungauged basins.

Conclusions
Satellite remote sensing and space-based technologies provide the most effective means to monitor dynamic changes in surface water resources, including river discharge, by continuously observing the spatiotemporal changes in multiple surface hydraulic parameters at a range of scales. In this study, a discharge estimation methodology based entirely on remote sensing-derived surface hydraulic parameters (including Surface Water and Ocean Topography (SWOT) observable quantities of width, depth and slope) were evaluated in a non-braided segment of Lhasa River around the Lhasa gauge site. The effective width was extracted from Landsat 5, 7 and 8 imageries covering the study period from 1988-2018, using the MNDWI approach after adjusting the threshold values to distinguish between water and non-water features. The channel resistance coefficient was assigned following Chow [77] and by visually inspecting and analyzing the Landsat images. The slope of the Lhasa River was estimated from a 30 m resolution SRTM DEM by considering an appropriate reach length. Average depth and velocity were then estimated using simplified empirical equations. Finally, river discharge was estimated using Model 1 and Model 2 by combining all the remotely sensed surface hydraulic parameters.
In spite of a lack of both observed cross sectional and satellite altimetry depth, the method developed in this study is notable in its success at estimating river depth with reasonable accuracy. This approach uses simplified empirical equations which include velocity, slope and roughness coefficient, based on Manning's equation and the general open-channel flow equations. We have also shown that the time series average velocity estimated from an empirical equation can be potentially used to estimate river discharge with reasonable accuracy. The estimated time series average width, depth and velocity of Lasa River at the Lhasa gauge site were 71.2 m, 1.44 m and 1.07 m s −1 , respectively. The discharge of the Lhasa River estimated by Model 1 and Model 2 generally corresponded very closely with the observed discharge at the Lhasa gauge site over the study period 1988-2018, in particular capturing the temporal changes between hydrological dry and wet seasons.
Both models showed a general tendency to overestimate discharge in moderate and high flows and to underestimate discharge in low flows. Model 2 discharge was always higher than that of Model 1, especially during high flows (discharge > 1000 m 3 s −1 ), it is to some extent magnified. Model 2 also minimizes the variance related with the channel roughness coefficient in its parametrization. The errors in Model 1 and Model 2 discharge estimates are largely associated with the uncertainties in model input parameters estimates (including river width, depth, slope and channel roughness coefficient). Errors in the discharge measurements, and other controlling factors, are also additional sources of uncertainties in discharge estimates. The overall performances of both Model 1 and Model 2 during the study period (1988-2018) were evaluated using different statistical metrics, yielding NSEs of 0.956 and 0.886, RRMSEs of 25.7 % and 41.4 %, MBEs of −7.5 m 3 s −1 and 9.5 m 3 s −1 , RMSEs of 46.2 m 3 s −1 and 74.5 m 3 s −1 and REs of −4.2% and 5.3% for Model 1 and Model 2, respectively. In general, river discharge estimated using Model 1 and Model 2 entirely from remote sensing-derived datasets for the Lhasa River, where it flows through a high mountainous region of the TP, was in close agreement with observed discharge, especially for low to moderate flows. Both models tended to overestimate discharge at high flows. This method presents a novel approach to estimating river discharge, without relying on observed hydraulic parameters, and could play a key role in river discharge estimates and understanding the hydrological processes: especially in poorly-gauged and ungauged basins. We further validated the transferability and performance of the methodology by estimating discharge in another reach of the Lhasa River, with different channel morphology, at Tanggya gauge (located about 85 km upstream of Lhasa gauge); in this case, both Model 1 and Model 2 showed similarly good performances with NSE values ≥ 0.950.