Discharge Estimation Using Integrated Satellite Data and Hybrid Model in the Midstream Yangtze River

: Remotely sensing data have advantages in ﬁlling spatiotemporal gaps of in situ observation networks, showing potential application for monitoring ﬂoods in data-sparse regions. By using the water level retrievals of Jason-2/3 altimetry satellites, this study estimates discharge at a 10-day timescale for the virtual station (VS) 012 and 077 across the midstream Yangtze River Basin during 2009–2016 based on the developed Manning formula. Moreover, we calibrate a hybrid model combined with Gravity Recovery and Climate Experiment (GRACE) data, by coupling the GR6J hydrological model with a machine learning model to simulate discharge. To physically capture the ﬂood processes, the random forest (RF) model is employed to downscale the 10-day discharge into a daily scale. The results show that: (1) discharge estimates from the developed Manning formula show good accuracy for the VS012 and VS077 based on the improved Multi-subwaveform Multi-weight Threshold Retracker; (2) the combination of the GR6J and the LSTM models substantially improves the performance of the discharge estimates solely from either the GR6J or LSTM models; (3) RF-downscaled daily discharge demonstrates a general consistency with in situ data, where NSE/KGE between them are as high as 0.69/0.83. Our approach, based on multi-source remotely sensing data and machine learning techniques, may beneﬁt ﬂood monitoring in poorly gauged areas. underestimations during the ﬂood peaks exist. The KGE and NSE between downscaling results and in situ data are as high as 0.83 and 0.69, with NRMSE and BIAS of 0.15 and 21%, respectively. In addition, the scatter plots between discharge simulations and observations strong linear relationship (CC = 0.87) between them.


Introduction
A flood is the sudden rise of water volume caused by heavy rain, rapid melting of ice/snow and storm surges [1]. As one of the most destructive weather-related hazards, flood has caused huge adverse impacts on infrastructures and the ecosystem [2,3]. Under climate warming and anthropogenic activities, frequency and intensity of floods have substantially intensified in the recent decades, suggesting a pressing need to better monitor flood processes for sustainable development. Gauge stations can provide accurate measurements of river discharge at point scale, facilitating both flood warning and water resources monitoring, which contribute to a better understanding of the water cycle process [4,5]. However, the distribution of hydrological stations is sparse, and the number of observation stations worldwide has decreased continuously due to considerable financial and material costs [6,7]. Considering the unique advantages that covering inaccessible regions and ones not affected by extreme hydrological events, the remote sensing technologies provide a preferable tool to evaluate river discharge, particularly in poorly gauged areas [8]. In this context, multiple remote sensing platforms have been used to return adequate streamflow observations. This study classified methods to generate discharge data using remote sensing data into two main categories. The first category exploits statistical relation between discharge and some indices sensitive to discharge regardless of the hydraulic features of rivers. This kind of method was initially developed by Brakenridge et al. [9], who calculated the ratio of Advanced Microwave Scanning Radiometer (AMSR-E) band values between merely focused on single data sources, the combination of Jason-2/3 altimetry satellites and the GRACE mission in this study will greatly improve the discharge estimation.
However, as simple and practical estimating discharge from space is, the discontinuous observations and coarse-resolution of remote sensing data challenge their application in flood monitoring. Hence, combining remote sensing data with hydrological models is fundamental to produce continuous discharge data for poorly gauged catchments. Michailovsky et al. [41] assimilated radar altimetry data to a routing model and improved its performance. Hulsman et al. [42] used altimetry observations and GRACE data to select parameter sets of a hydrological model in a data-scarce catchment of Zambia. Assimilation of future pre-SWOT data clearly corrected hydrological models on a global and continental scale, by reducing discharge errors from~30% to~24% [43]. The remote sensing-enabled hydrological simulation has been conducted across many basins [44,45]. Furthermore, revolutionary advances in machine learning methods have demonstrated great potential in streamflow prediction. For example, Cheng et al. [46] established three machine learning models involving Artificial Neural Network (ANN), Support Vector Regression, and Long Short-Term Memory (LSTM) to predict discharge fluctuation in North China, and satisfactory accuracy was obtained. Gauch et al. [47] applied the multi-timescale LSTM model to 516 basins over the continental United States and achieved a significantly better runoff simulation performance than the US National Water Model. Machine learning models have been used to discharge modeling with remarkable success [48]. Several studies developed the so-called hybrid modeling approach by coupling physical-models with machine learning techniques to improve the modeling performance [49,50]. The so-called hybrid models exploit the predictive ability of physical models, then harvest available data to correct modeling bias [51]. Alternatively, the flood generation is a complex interacting phenomena, which is inherently characterized by peak discharge, flood volume and hydrograph shape, particularly in natural hazard risk assessment and water resources management [52]. The daily streamflow extremes, having a high correlation with peak discharge, represent a key readily available measure of extreme storms and are frequently used to assess flood hazard risk. However, most of the existing studies only derive the ten-day streamflow, which cannot be used for representing a flood event.
Although a few studies have estimated discharge using either remotely sensing data or a hybrid model across many basins, some shortcomings still exist, including: (1) their joint application has not been systematically validated; (2) the temporal resolution is too coarse for water management and water risk assessment. Therefore, the objectives of this study are to: (1) estimate river discharge using water level retrievals from satellite altimetry based on the improved Multi-subwaveform Multi-weight Threshold Retracker and the developed Manning formula; (2) develop a hybrid hydrological model using satellite-based discharge and other remotely sensing datasets (e.g., GRACE TWSA, rainfall) by coupling a GR6J model with a LSTM model; (3) statistically downscale the hydrological simulations into a daily scale using the Random Forest method.

Study Area
The Yangtze River Basin controls a drainage area of about 1,800,000 km 2 . Highmountain regions dominate the upper reach of the Yangtze River Basin and the very high-altitude results in the problem that the antenna needs to be very large in order to distinguish individual radar waveforms reflected by the water surface. In addition, the lower reach of the Yangtze River Basin characterized by a developed shipping and port industry [53] is likely to affect the reflection of echoes with high uncertainty [54]. Free from the impact of extremely high-altitude and intense human activities, this study focused on the middle reach of the Yangtze River Basin located between Zhicheng and Huangshigang (ZHY) because the river here is generally wider than 1000 m, and so it is feasible to extract water surface signals from radar waves (see Figure 1). Moreover, the ZHY is located downstream from the Three Gorges Dam with a total storage capacity of 39.3 km 3 , facing high flood risk in the context of increasing precipitation extremes. The flood intensity has changed from decadal to inter-decadal cycles due to the rapid expansion of urban populations, high shrinkage of river and lake systems, and fast construction of water facilities [55]. In addition, changes in flood frequency show an overall increasing trend in the ZHY over the past 50 years [56]). Hydrological simulation using remote sensing technologies is of great significance for the reduction of financial and life losses in flood disasters, particularly in the data-limited region [57]. This study divided the ZHY into three sub-basins, consisting of the poorly gauged Hanjiang River Basin (sub-basin 1), the Dongting Lake Basin (sub-basin 3) with sparse stations, and the remaining area within the ZHY (sub-basin 2). Moreover, two ground tracks of Jason-2/3 satellites were selected to derive the water level at the virtual station (VS) 012 and the VS077 along the Yangtze River.
The hydrological records at the Zhicheng station near the VS012 and the Huangshigang station near the VS077 were used for comparison (see Figure 1). because the river here is generally wider than 1000 m, and so it is feasible to extract water surface signals from radar waves (see Figure 1). Moreover, the ZHY is located downstream from the Three Gorges Dam with a total storage capacity of 39.3 km 3 , facing high flood risk in the context of increasing precipitation extremes. The flood intensity has changed from decadal to inter-decadal cycles due to the rapid expansion of urban populations, high shrinkage of river and lake systems, and fast construction of water facilities [55]. In addition, changes in flood frequency show an overall increasing trend in the ZHY over the past 50 years [56]). Hydrological simulation using remote sensing technologies is of great significance for the reduction of financial and life losses in flood disasters, particularly in the data-limited region [57]. This study divided the ZHY into three sub-basins, consisting of the poorly gauged Hanjiang River Basin (sub-basin 1), the Dongting Lake Basin (sub-basin 3) with sparse stations, and the remaining area within the ZHY (subbasin 2). Moreover, two ground tracks of Jason-2/3 satellites were selected to derive the water level at the virtual station (VS) 012 and the VS077 along the Yangtze River. The hydrological records at the Zhicheng station near the VS012 and the Huangshigang station near the VS077 were used for comparison (see Figure 1).

Methods
The overall implementation details of our experiments are demonstrated in Figure 2. We aim to apply a satellite-derived water level to estimate river discharge and further calibrate the GR6J hydrological model to produce continuous discharge; then, improve the modeling results using the Long Short-Term Memory model combined with GRACE TWSA before downscaling the results into a daily time scale based on the random forest method. Different from the conventional offline downloading method, Google Earth Engine (GEE) was used to extract regional precipitation, potential evaporation, and air temperature in the ZHY from GPM and ERA-interim datasets. GEE is a freely available cloud computing platform for scientific analysis and visualization of geospatial datasets. Specifically, the open-source QGIS software is used for spatial analysis and map preparation and the GEE Code Editor requires the JavaScript Application Programming Interface (API) to communicate with the data catalogue. The powerful cloud servers and petabytes, and publicly available remote sensing archived datasets of GEE, enable a planetary-scale

Methods
The overall implementation details of our experiments are demonstrated in Figure 2. We aim to apply a satellite-derived water level to estimate river discharge and further calibrate the GR6J hydrological model to produce continuous discharge; then, improve the modeling results using the Long Short-Term Memory model combined with GRACE TWSA before downscaling the results into a daily time scale based on the random forest method. Different from the conventional offline downloading method, Google Earth Engine (GEE) was used to extract regional precipitation, potential evaporation, and air temperature in the ZHY from GPM and ERA-interim datasets. GEE is a freely available cloud computing platform for scientific analysis and visualization of geospatial datasets.
Specifically, the open-source QGIS software is used for spatial analysis and map preparation and the GEE Code Editor requires the JavaScript Application Programming Interface (API) to communicate with the data catalogue. The powerful cloud servers and petabytes, and publicly available remote sensing archived datasets of GEE, enable a planetary-scale environmental data analysis, greatly reducing the computation cost of academic users [58].
is the elevation of the river bottom, and denotes the water level retrievals from satellite radar altimetry. Parameters a and ℎ are calibrated using the least square method combined with a few in situ discharges from between 2008 and 2010.
The trained Manning formula was employed to estimate discharge using water levels from satellite altimetry. It is noted that the developed Equation (5) is only feasible for poorly gauged areas with prior information in bed level and cross-sections of the river but not for ungauged sites.

Hybrid Model
The GR6J hydrological model is a lumped, conceptual, and parsimonious model with six parameters to identify, which is simple and feasible that only requires precipitation, air temperature and potential evapotranspiration as inputs, and has been widely used due

Improved Multi-Subwaveform Multi-Weight Threshold Retracker
To avoid land contaminations as much as possible, this study exploited the highresolution cloud-free Landsat-8 images of the VS012 and VS077 between November and April to exclude footprints outside the river bank, which is called footprint selection in Figure 3. Moreover, the raw waveform selection was conducted using the pulse peakiness (PP) index proposed by Strawbridge and Laxon [59] to select the waveforms that show dominant peaks and consistent peak locations. The PP index is defined as follows: where P max denotes the maximum power of the waveform, and P(i) represents the return power of the i th (i = 1, 2, . . . , n) sampling gate; n is the number of sampling gates (104 for Jason-2/3). The waveforms with a PP index lower than 1.2 are discarded due to the difficulty of extracting water surface height information from them [60]. Yuan et al. [60] proposed the Multi-subwaveform Multi-weight Threshold Retracker (MSMWTR) method to improve the quality of satellite altimetry data. The MSMWTR method first identifies subwave forms, then derives a retracking gate using the threshold algorithm for each one [61] and obtains the subwaveform retracked ranges for water level extraction. The major novelty of the MSMWTR method is calculating the weight for each subwaveform as Equation (2); hence, the abnormal outliers can be removed by inspecting the relative weight of each subwaveform artificially.
where m is the number of subwaveforms; p j k represents the return power at the sampling gate j of the k th subwaveform. where and are the simulated and observed discharge, respectively. Similarly, / and s ( )/ ( ) represent the mean value and the standard deviation of / , respectively. is the number of discharge simulations. The NRMSE and BIAS were used to quantify the relative deviations of the estimates from the observations considering the distance between the VS and hydrological stations. The CC ranging from -1 to 1, NSE ranging from -∞ to 1 and KGE changes between 0 and 1 were applied for evaluation of the model performance. The closer CC, NSE and KGE were to 1, the better the models perform.  Although the MSMWTR method performed well in the middle reach of the Yangtze River [60], the approach relied heavily on manual intervention to remove the outlies. Hence, this study can improve the original MSMWTR method as follows: (1) Set an empirical threshold (0.5) after repeated experiments.
(2) Implement the MSMWTR algorithm and exclude the subwaveform retracked ranges whose weight is lower than the threshold for each cycle. The improved MSMWTR method relies less on manual intervention and has the potential to be used in poorly gauged regions, as any prior knowledge about the river is not required. After resolving the satellite altimetry range, water level retrievals can be calculated as follows: where H represents the derived water level (m, asl), H sat means the altitude of the satellite, and R is the retracked range. Cor dry , Cor wet , Cor iono , Cor set , and Cor pt are the corrections of dry tropospheric, wet tropospheric, ionospheric, solid Earth tide, and polar tide, respectively. EGM96 denotes the EGM96 geoid. All correction items, the altitude of the satellite, and the EGM96 geoid are included in the SGDRs product of Jason-2 and Jason-3 satellites.

Developed Manning Formula
The Manning formula can be expressed as follows: where Q is the discharge (m 3 /s) and n represents the Manning roughness factor; A denotes the cross-section area of the flow (m 2 ), R is the hydraulic radius (m) and S denotes the slope of the energy gradient of the river channel. The Manning formula considers the hydraulic geometry of the river and estimates river discharge simply and feasibly. Figure 2 presents the measured cross-sections of the Zhicheng station (near the VS012) and the Huangshigang station (near the VS077). We assume that the cross-sections between the two locations keep a similar width-to-depth radio and approximately belong to the trapezoidal shape, and the bed level in the ZHY is known. Moreover, based on the fact that an accurate estimation of discharge can be achieved without inputting the Manning roughness factor and the difficulty of measuring the changing slope of the energy gradient along the river channel [18], the slope and roughness coefficient in the Manning formula were considered as constants. Hence, the developed Manning formula by Huang et al. [63] for trapezoidal cross-section was applied as follows: where a represents the influence of water surface slope and the roughness coefficient, h is the elevation of the river bottom, and H denotes the water level retrievals from satellite radar altimetry. Parameters a and h are calibrated using the least square method combined with a few in situ discharges from between 2008 and 2010. The trained Manning formula was employed to estimate discharge using water levels from satellite altimetry. It is noted that the developed Equation (5) is only feasible for poorly gauged areas with prior information in bed level and cross-sections of the river but not for ungauged sites.

Hybrid Model
The GR6J hydrological model is a lumped, conceptual, and parsimonious model with six parameters to identify, which is simple and feasible that only requires precipitation, air temperature and potential evapotranspiration as inputs, and has been widely used due to the comparable performance to other more complex models over various catchments worldwide [64,65].
The LSTM model is capable of overcoming the weakness of the traditional recurrent neural network (e.g., vanishing and exploding gradients) to learn long-term dependencies. As a special kind of recurrent neural network structure, the LSTM model is one of the most popular neural network models in nonlinear time series analyses because it can store and relate previous information in a sequence, enabling it to predict time series in near real-time. LSTM has been widely used in hydrological simulation and prediction, extreme flood monitoring, and precise water resource management [66,67].
Combining hydrological simulation with machine learning technologies has been validated effective to improve the accuracy of stream flow simulations [68]. The discharge modeling based on the hybrid model was carried out from 2009 to 2016 due to the data gap and the decomposition of the GRACE mission. As illustrated in Figure 2, the rainfall-runoff model GR6J [69] was firstly applied for regional stream flow modeling for the whole of ZHY. The input variables were: the satellite-based discharge at the VS012 that flows into the sub-basin 2 using the developed Manning formula, IMERG Final precipitation, air temperature and potential evaporation from ERA5. The discharge estimates based on the developed Manning formula at the VS077 was used to calibrate the GR6J model. The datasets from 2009 to 2012 (~%50 of samples) were used for training the GR6J hydrological model and those from 2013 to 2016 (~%50 of samples) were used for testing the model's performance. Subsequently, the Long Short-Term Memory model (LSTM) method was employed to improve the discharge modeling accuracy of the GR6J hydrological model. The inputs of the LSTM model compromise GRACE TWSA for the sub-basin 1 and subbasin 3 representing the water volume flowing into the mainstream of the Yangtze River, GR6J-modeled discharge at the VS077, and remote sensing meteorological forcing data (i.e., GPM precipitation, ERA air temperature, and ERA potential evapotranspiration) for the whole of the ZHY. Satellite-based discharge at the VS077 was still used to train the LSTM model for direct comparison with the simulations not coupled with the LSTM model.
Similarly, the datasets between 2009 and 2013 were spilt as training (2009 to 2012,~%50 of samples) and testing (2013 to 2016,~%50 of samples) periods, respectively. In addition, we designed another experiment solely using the LSTM model to evaluate the performance improvement of the hybrid model compared with the GR6J and LSTM model, respectively. We note that the datasets and experiment setting of the single LSTM model were the same as the GR6J hydrological model. Moreover, for the single GR6J and LSTM model, as well as the hybrid model, the flow discharge data at the Huangshigang gauge station (near the VS077) were used for comparison.

Random Forest Model
The 10-day discharge time series from the hybrid model is generally too infrequent for flood simulation. Therefore, this study applied the random forest (RF) model to downscale simulated discharge data into the daily scale based on the strong statistical correlation between river height from spaceborne radar altimeters and the difference between daytime and nighttime land surface temperatures [70]. Random forest (RF) regression is an ensemble learning method based on classification and regression trees (CART). Each tree computes bootstrap samples chosen randomly from the predictor dataset at each node so the RF model can produce reasonable predictions by averaging the results of CARTs [71,72]. As the RF models do not require specific distributional assumptions of the original dataset and have an easy implementation, they have been widely used in the statistical downscaling field [73]. This study established the RF model between simulated discharge and three predictors, including the maximum skin temperature (i.e., at 14 p.m.), the minimum skin temperature (i.e., at 4 a.m.), together with the difference between them (see Figure 3).

Performance Metrics
Discharge estimation performance metrics were accessed by the correlation coefficient (CC), Nash-Sutcliffe efficiency coefficient (NSE) [74], the Kling-Gupta efficiency (KGE) [75], the normalized root mean square error (NRMSE), and the relative bias (BIAS): where Q s and Q o are the simulated and observed discharge, respectively. Similarly, Q s /Q o and std(Q s )/ std(Q o ) represent the mean value and the standard deviation of Q s /Q o , respectively. n is the number of discharge simulations. The NRMSE and BIAS were used to quantify the relative deviations of the estimates from the observations considering the distance between the VS and hydrological stations. The CC ranging from -1 to 1, NSE ranging from -∞ to 1 and KGE changes between 0 and 1 were applied for evaluation of the model performance. The closer CC, NSE and KGE were to 1, the better the models perform.

Jason-2/3 Missions
As the successor of Jason-1, Jason-2 was launched in 2008 with the Poseidon-3 radar altimeter on board. The sensor geophysical data records (SGDRs) of the Jason-2 satellite from the CNES Archiving, Validation and Interpretation of Satellite Oceanographic (AVISO+) team, including ranges, waveforms, geo-physical, and atmospheric corrections during 2008-2016 were obtained to retrieve water level at a high sampling rate of 20 Hz. The waveforms from the SGDRs were retracked using the improved MSMWTR algorithm proposed in this study [60] (refer to Section 2.2.1 for details) to improve the accuracy of water level retrievals at the VS012 and VS077 (see Figure 1). Jason-3 inherited the orbit of Jason-2 after its launch in 2016, with Jason-2 being shifted into an interleaved orbit and continuing functioning until 2019; this, increased the spatial coverage of the Jason altimetry series to some degree, with Poseidon-3B being the main instrument [76]. The water level retrievals of water level at the VS012 and VS077 between 2016 and 2019 were obtained from Jason-3 SGDRs data, and their waveforms were retracked using the improved MSMWTR algorithm for higher accuracy. The revisit period of the Jason-2/3 satellites is 10 days; hence, the sub-monthly knowledge of water level can be derived.

ITSG-Grace 2018 Product
This study used daily GRACE TWSA to improve hydrological simulations (see Figure 3). Solved by Graz University of Technology, the ITSG-Grace2018 gravity field product provides unrestricted daily GRACE solutions processed from GRACE-measured Level-1B data. As the GRACE data coverage within one day is not sufficient to resolve the global gravity field, the daily ITSG-Grace2018 solution was stabilized using an autoregressive model of order three in order to recover sub-monthly gravity variations, which was derived by fitting coefficients to the de-aliasing error estimates contained in the ESA's Earth system model and the de-trended residuals of land surface discharge model [77]. For each day of the observation period, a set of spherical harmonic coefficients for degrees ranging from 2 to 40 were estimated. The adjustment delivered daily solutions, even if there were no GRACE data available for a specific day. The daily solution shares the same signal definition as the official monthly GRACE products. They contain the full hydrological, cryospheric and GIA signal. The atmosphere and ocean signals are relative to AOD1B RL06 [78,79]. The daily TWSA data were resampled with 10-day intervals to be consistent with Jason-2/3 data.

Forcing Data for the Hydrologic Model
The daily remote sensing precipitation dataset used for hydrological simulation is from the Global Precipitation Measurement (GPM) mission, which is an international network of satellites that provide next-generation global observations of rain and snow. The Core Observatory satellite of the GPM, launched in February 2014, carried out the Dual-frequency Precipitation Radar and GPM Microwave Imager instruments to generate a global precipitation dataset. The Integrated Multi-Satellite Retrievals (IMERG) algorithm was employed to merge microwave-calibrated infrared satellite estimates, satellite microwave precipitation estimates, and in situ precipitation data with fine accuracy and consistency all over the globe [80]. The GPM IMERG Final version 06 product was applied with a 3-month latency. Moreover, the daily precipitation dataset from the GPM mission was with a high spatial resolution of 0.1 • × 0.1 • , which showed good applicability and high accuracy in the Yangtze River Basin [81]. Alternatively, daily potential evaporation (mm) and air temperature (K) from ERA-interim reanalysis datasets were used in this study for hydrological modeling, which contains physical data for atmospheric and surface analyses since 1979 based on the European Centre for Medium-Range Weather Forecasts' Integrated Forecast System Release Cy31r2 [82]. The daily precipitation, potential evaporation, and air temperature data were resampled with 10-day intervals, with radar altimetry observations. We also calculated the maximum and the minimum air temperature of the day from hourly ERA-interim product to downscale 10-day period discharge into a daily timescale using the random forest method [70].

In Situ Measurement
We collected in situ observations of water level and river discharge at the Zhicheng station (near the VS012, see Figure 1) and the Huangshigang station (near the VS077, see Figure 1). The in situ data follow a strict quality control policy. The measurements were used mainly to validate water level retrievals from Jason-2/3 data and corresponding discharge estimations. In addition, a 3-year discharge time series was applied to fit hydrodynamic parameters in the developed Manning formula (refer to Section 2.2.2 for details). Figure 4 presents satellite-derived water level time series at the VS012 and the VS077, and water level measurements at the Zhicheng station (near the VS012) and the Huangshigang station (near the VS077) were used for comparison. The results show that water level retrievals of Jason-2/3 satellites are generally consistent with observed series, particularly in the wet season between May and October as the waveforms are less noisy with obvious dominant peaks and consistent peak locations [60]. Contrarily, the obvious underestimation during the dry season between November and April suggests more noise and various peak locations among waveforms. We furthermore established the developed Manning formula using the satellite-based water level and a few in situ observed discharges between 2008 and 2010; simulated discharge between 2011 and 2019 were used for testing. Figure 5 shows the calibration and validation results of the developed Manning formula and relative performance metrics are summarized in Table 1.    At the VS012, the trained model overestimates discharge during the dry season with a CC of 0.88 between 2008 and 2010, and the NSE and KGE reaches 0.74 and 0.83, respectively (see Table 1). The poorer model performance during the test period is because of the underestimation of streamflow peaks during the wet season (CC/NSE/KGE = 0.75/0.66/0.70). However, the developed Manning formula presents better performance at the VS077 with CC, NSE and KGE as high as 0.97, 0.93 and 0.96 during the calibration period. The trained model illustrates continuous satisfactory performance during the test period according to Table 1, although some underestimation during the dry season still occurs because of the water level bias. Figure 6  relationship (CC = 0.96) between estimates of discharge and observations is discovered at the VS077, consistent with the calibration and validation results for the VS012 and VS077.   Table 1). The poorer model performance during the test period is because of the underestimation of streamflow peaks during the wet season (CC/NSE/KGE = 0.75/0.66/0.70). However, the developed Manning formula presents better performance at the VS077 with CC, NSE and KGE as high as 0.97, 0.93 and 0.96 during the calibration period. The trained model illustrates continuous satisfactory performance during the test period according to Table 1, although some underestimation during the dry season still occurs because of the water level bias. Figure 6 presents the scatter plot between the satellite-estimated discharge and the in situ measurements during the whole period during 2009-2019. The results indicate that the fitted linear regression line at the VS012 is with a slope of 0.88 and the corresponding correlation coefficient of 0.91, and a stronger linear  relationship (CC = 0.96) between estimates of discharge and observations is discovered at the VS077, consistent with the calibration and validation results for the VS012 and VS077.  Figure 7 presents discharge estimates from the GR6J, LSTM and hybrid models as well as observations for the ZHY from 2009 to 2016. The comparison indicates that the LSTM model outperforms the GR6J model, especially during the dry season. In addition, discharge data from the hybrid model are generally in better agreement with in situ data than that solely from either the GR6J or LSTM model, which is consistent with the scatter plot and Tayler diagram between simulations and observations in Figure 8. A stronger correlation was discovered between in situ discharge and hybrid model results (CC = 0.96) than the GR6J results (CC = 0.92) and LSTM data (CC = 0.92). Table 2 Figure 7 presents discharge estimates from the GR6J, LSTM and hybrid models as well as observations for the ZHY from 2009 to 2016. The comparison indicates that the LSTM model outperforms the GR6J model, especially during the dry season. In addition, discharge data from the hybrid model are generally in better agreement with in situ data than that solely from either the GR6J or LSTM model, which is consistent with the scatter plot and Tayler diagram between simulations and observations in Figure 8. A stronger correlation was discovered between in situ discharge and hybrid model results (CC = 0.96) than the GR6J results (CC = 0.92) and LSTM data (CC = 0.92). Table 2 Figure 9a demonstrates the comparison between downscaled discharge simulated by the hybrid model at the VS077 and in situ data at the Huangshigang station from 2009 to 2016. The daily streamflow is generally consistent with in situ results although some underestimations during the flood peaks exist. The KGE and NSE between downscaling results and in situ data are as high as 0.83 and 0.69, with NRMSE and BIAS of 0.15 and 21%, respectively. In addition, the scatter plots between discharge simulations and observations in Figure 9b indicate a strong linear relationship (CC = 0.87) between them.

Discussion
A few previous studies have simulated discharge using either remote sensing data or hydrological models in the Yangtze River Basin. Xu et al. [83] used very high-resolution QuickBird-2 data to establish a width-discharge curve and estimate flow discharge in the Yangtze River with high accuracy. However, the method heavily relies on the spatial resolution and data quality of satellite images. Based on the accurate water level retrievals of

Discussion
A few previous studies have simulated discharge using either remote sensing data or hydrological models in the Yangtze River Basin. Xu et al. [83] used very high-resolution QuickBird-2 data to establish a width-discharge curve and estimate flow discharge in the Yangtze River with high accuracy. However, the method heavily relies on the spatial resolution and data quality of satellite images. Based on the accurate water level retrievals of satellite altimetry data over the Yangtze River [84], Sichangi et al. [85] initially combined satellite altimetry-derived water level and MODIS-derived river width to estimate flow discharge at the Datong gauge station in the Yangtze River between 2003 and 2008. However, the accuracy of discharge estimations is relatively lower (NSE =~0.5) than this study, with NSE ranging from 0.66 to 0.93. Chen et al. [86] also successfully estimated total basin discharge solely based on remote sensing data, including TRMM precipitation, GRACE TWSA, and MODIS images of the Yangtze River Basin during 2003-2012. Although achieving a reasonable accuracy with correlation coefficient of 0.89 and RMSE of~5900 m 3 /s, the monthly discharge estimations are too infrequent to support flood simulation. Alternatively, the discharge prediction ability of hydrological models has been comprehensively validated across many sub-basins of the Yangtze River Basin [87,88]. However, these works mainly applied observed discharge to calibrate the models; hence, experiments in some poorly gauged regions (e.g., upper reach of the Yangtze River Basin) are rarely carried out previously. In addition, Xu et al. [89] applied the LSTM model to model river flow in the upper reach of the Yangtze River and achieved better performance compared with several hydrological models. Liu et al. [90] tested several deep learning neural networks, including the LSTM and the Empirical Mode Decomposition model in the ZHY and achieved good forecasting accuracy in terms of NSE and RMSE because of the powerful capacity in learning non-linear and complex processes of machine learning technologies. Nevertheless, the hydrological models and machine learning methods are mostly employed for comparison rather than combination.
Even though our developed approach successfully generated daily outlet discharge for the ZHY between 2009 and 2016 using remote sensing data and the hybrid model, some limitations still exist as follows: (1) Satellite radar altimetry observations are difficult to extract the water height signals of the rivers located on high-mountain regions due to the very high-altitude and complicated topography, which challenge the applications of our algorithm in such kind of basins. Huang et al. [63] developed a combined waveform retracker to retrieve water levels in the upstream Brahmaputra River Basin with high mountains, while the method showed poor performance in the Yangtze River due to intense human activities. Moreover, the coarse spatial resolution of GRACE TWSA restricts the discharge estimation ability of our approach in medium-sized basins [91]. Advanced post-processing methods, including mass concentration blocks and constrained forward modeling, can enable the small-scale hydrological applications such as: flood and drought monitoring, groundwater depletion detection, and mass balance of Greenland and Antarctica [92]. It is still interesting to establish the skill of the adopted algorithm in smaller rivers that has higher confidence in discharge measurements, which will be further investigated and fully evaluated in the near future after solving the above-mentioned problems. (2) Multistep post-processing dynamics of our developed approach generated uncertainty in data, parameters, and structures of the developed Manning formula, GR6J model, and machine learning methods. Quantifying uncertainty contribution of each model and dataset is fundamental for our better understanding of hydrological forecasting mechanism, which will be further investigated and fully evaluated in the near future using uncertainty analysis tools such as analysis of variance [93]. (3) Due to the relatively short operation life of the Jason-2/3 satellite from 2009 up to present day, and the GRACE mission between April 2002 and June 2017, the length of the discharge series in this study is limited. Combining multi-source satellite data, including Envisat, TOPEX/Poseidon, and SARAL missions may be beneficial to overcome this shortage. Moreover, the SWOT satellite that is expected to launch in late 2021 will produce high-resolution surface water knowledge, including water level and slope, providing a great innovation for estimating discharge [36]. In addition, the in-orbit GRACE Follow-On satellite has provided accurate monthly TWSA data around the world since June 2018. The reconstruction of the historical satellite data and accumulation of the future satellite data will greatly extend our study time span [94].

Conclusions
This study estimated discharge using the water level retrievals from satellite radar altimetry based on the developed Manning formula at the VS012 and VS077 across the Yangtze River between 2008 and 2019, and further used the satellite-based discharge to calibrate the GR6J model. Subsequently, the LSTM method was employed to improve the hydrological modeling results combined with GRACE TWSA. Finally, the 10-day discharge time series at the outlet of the ZHY was statistically downscaled into a daily scale using the RF approach. The main conclusions of this study are summarized as follows: (1) The developed Manning formula shows an effective discharge estimation ability at the VS012 and the VS077. The fitted linear regression line between discharge estimates and in situ observations at the VS012 is with a slope of 0.88 and CC of 0.91 between them. A stronger linear relationship between estimated and observed flow discharges was discovered at the VS077 (k/CC = 0.95/0.96). (2) The hybrid model clearly shows the improvements of discharge estimation in comparison with either the sole GR6J or LSTM models. The CC, KGE and NSE of the hybrid model increased by 3.4/3.4%, 1.1/2.1% and 8.4/3.4% more than the GR6J/LSTM models, respectively. In addition, a noticeable decrease in NRMSE of 14.3/25.0% and BIAS of 13.3/23.5% was detected by comparing with the GR6J/LSTM models, respectively.