Machine Learning Improvement of Streamﬂow Simulation by Utilizing Remote Sensing Data and Potential Application in Guiding Reservoir Operation

: Hydro-meteorological datasets are key components for understanding physical hydrological processes, but the scarcity of observational data hinders their potential application in poorly gauged regions. Satellite-retrieved and atmospheric reanalysis products exhibit considerable advantages in ﬁlling the spatial gaps in in-situ gauging networks and are thus forced to drive the physically lumped hydrological models for long-term streamﬂow simulation in data-sparse regions. As machine learning (ML)-based techniques can capture the relationship between different elements, they may have potential in further exploring meteorological predictors and hydrological responses. To examine the application prospects of a physically constrained ML algorithm using earth observation data, we used a short-series hydrological observation of the Hanjiang River basin in China as a case study. In this study, the prevalent mod è le du G é nie Rural à 9 param è tres Journalier (GR4J-9) hydrological model was used to initially simulate streamﬂow, and then, the simulated series and remote sensing data were used to train the long short-term memory (LSTM) method. The results demonstrated that the advanced GR4J9–LSTM model chain effectively improves the performance of the streamﬂow simulation by using more remote sensing data related to the hydrological response variables. Ad-ditionally, we derived a reservoir operation model by feeding the LSTM-based simulation outputs, which further revealed the potential application of our proposed technique.


Introduction
The availability of reliable hydro-meteorological material is an initial yet crucial part of water resource planning and management. Using hydrological simulation (or forecasting) as an example can have significant repercussions on socio-economic growth and development prospects from rainfall-runoff observations [1][2][3][4][5]. However, the scarcity of hydro-meteorological monitoring and inaccessibility issues pose obstacles to conducting effective integrated evaluation research, developing modeling frameworks, and recommending policies for resilience, especially for developing countries. Obtaining access to actual long-series hydro-meteorological processes is worthy of further investigation.
In recent decades, both satellite telemetry and data inversion techniques have been mined deeply, which compensate for the deficiencies of meteorological stations and provide

Study Area
The Hanjiang River, illustrated in Figure 1, was chosen as the case study. The river is the largest tributary of the Yangtze River. It lies between 30.28 • and 34.5 • N and 106.42 • and 114.55 • E, with a mainstream length of 1577 km and a total drainage area of 159,000 km 2 . It originates from the southern Qin Mountain, flows through the Shanxi and Hubei provinces and converges into the Yangtze River in Wuhan. Characterized by a subtropical monsoon climate and annual precipitation of between 700 and 1100 mm, this basin has abundant water resources; 75% of the total annual precipitation occurs in the flood season (June to September). During the flood season, the sudden rainstorms in early summer and persistent rainfalls in autumn typically induce large-scale flooding [25]. For water resource regulation, a key water conservancy project named the Danjiangkou Reservoir was built in the middle of the Hanjiang River Basin. It not only serves as a source for the Middle Route of the South-to-North Water Diversion Project (MSWDP) but also plays an important role in China's One Belt, One Road construction. Therefore, the sub-basin over Danjiangkou Reservoir is a useful candidate for conducting our proposed approach (in Section 3).

Study Area
The Hanjiang River, illustrated in Figure 1, was chosen as the case study. The river is the largest tributary of the Yangtze River. It lies between 30.28° and 34.5° N and 106.42° and 114.55° E, with a mainstream length of 1577 km and a total drainage area of 159,000 km 2 . It originates from the southern Qin Mountain, flows through the Shanxi and Hubei provinces and converges into the Yangtze River in Wuhan. Characterized by a subtropical monsoon climate and annual precipitation of between 700 and 1100 mm, this basin has abundant water resources; 75% of the total annual precipitation occurs in the flood season (June to September). During the flood season, the sudden rainstorms in early summer and persistent rainfalls in autumn typically induce large-scale flooding [25]. For water resource regulation, a key water conservancy project named the Danjiangkou Reservoir was built in the middle of the Hanjiang River Basin. It not only serves as a source for the Middle Route of the South-to-North Water Diversion Project (MSWDP) but also plays an important role in China's One Belt, One Road construction. Therefore, the sub-basin over Danjiangkou Reservoir is a useful candidate for conducting our proposed approach (in Section 3). The Danjiangkou Reservoir has eased chronic water shortages in several of China's provinces and urban cities, including the capital, Beijing. An official water diversion diagram developed by the Ministry of Water Resources of China is used for guidance of the water diversion. As shown in Figure 2, it defines a pre-set water diversion value. For example, if the reservoir water level at the initial time of the water diversion is in Region 3 (in Figure 2), the ideal water diversion flow should be 300 m 3 /s. Apart from water diversion, the Danjiangkou Reservoir also works as a hydropower source. Water supply and hydropower consist of the main positive purposes of the Danjiangkou Reservoir; these two objectives compete with each other, as part of the reservoir water release is redirected for water diversion instead of power generation. The basic reservoir parameters are listed in Table 1.  The Danjiangkou Reservoir has eased chronic water shortages in several of China's provinces and urban cities, including the capital, Beijing. An official water diversion diagram developed by the Ministry of Water Resources of China is used for guidance of the water diversion. As shown in Figure 2, it defines a pre-set water diversion value. For example, if the reservoir water level at the initial time of the water diversion is in Region 3 (in Figure 2), the ideal water diversion flow should be 300 m 3 /s. Apart from water diversion, the Danjiangkou Reservoir also works as a hydropower source. Water supply and hydropower consist of the main positive purposes of the Danjiangkou Reservoir; these two objectives compete with each other, as part of the reservoir water release is redirected for water diversion instead of power generation. The basic reservoir parameters are listed in Table 1

Data Collection
Three kinds of datasets (i.e., satellite-based observation, atmospheric reanalysis, and short-series streamflow) were collected and used in this study. The GPM Core Observatory is equipped with the first space-borne Ku/Ka-band dual-frequency radar and a multichannel microwave imager, which improves the monitoring ability of light and solid precipitation. Since the first release of Integrated multi-satellite retrievals for GPM (IMERG) products in 2015, it has undergone many improvements, and the latest version (V06B) has been retrospectively processed, including TRMM-era data since June 2000. Due to the infusion of the Global Precipitation Climatology Centre (GPCC) rain gauge data, the final operation of IMERG provides a more accurate estimation and was therefore adopted in this study.
ERA5 was used as another meteorological product, which is a global atmospheric reanalysis dataset developed by ECMWF. ERA5 data are generated by the combination of model simulations and observations using physics laws, which are based on data assimilation by the Integrated Forecasting System (IFS Cy31r2). This assimilation system includes a four-dimensional variational (4D-Var) analysis method and considers the exact time of observation and model evolution in the assimilation window to estimate the deviation between observations and select high-quality data from poor data. The hourly output resolution is 0.25° × 0.25°, which offers a more sophisticated simulation of weather

Data Collection
Three kinds of datasets (i.e., satellite-based observation, atmospheric reanalysis, and short-series streamflow) were collected and used in this study. The GPM Core Observatory is equipped with the first space-borne Ku/Ka-band dual-frequency radar and a multichannel microwave imager, which improves the monitoring ability of light and solid precipitation. Since the first release of Integrated multi-satellite retrievals for GPM (IMERG) products in 2015, it has undergone many improvements, and the latest version (V06B) has been retrospectively processed, including TRMM-era data since June 2000. Due to the infusion of the Global Precipitation Climatology Centre (GPCC) rain gauge data, the final operation of IMERG provides a more accurate estimation and was therefore adopted in this study.
ERA5 was used as another meteorological product, which is a global atmospheric reanalysis dataset developed by ECMWF. ERA5 data are generated by the combination of model simulations and observations using physics laws, which are based on data assimilation by the Integrated Forecasting System (IFS Cy31r2). This assimilation system includes a four-dimensional variational (4D-Var) analysis method and considers the exact time of observation and model evolution in the assimilation window to estimate the deviation between observations and select high-quality data from poor data. The hourly output resolution is 0.25 • × 0.25 • , which offers a more sophisticated simulation of weather processes. The hourly near-surface air temperature, dew point temperature (T dew ), and wind speed from the ERA5 dataset were considered. All the sub-daily satellite/reanalysis data covering 2002-2019 were aggregated into a daily scale.
As the boundary of the upper and mid-lower reaches of the Hanjiang River Basin, the short-series inflow of the Danjiangkou Reservoir was selected to calibrate the parameters of the hydrological model. Observed streamflow data spanning 2003-2007 were obtained from the Yangtze River Water Conservancy Commission.

Methodology
A flowchart of the framework module is presented in Figure 3, which is further elaborated in the following sections. It is worth mentioning that our proposed framework can be used in hydrological simulation rather than forecasting. processes. The hourly near-surface air temperature, dew point temperature (Tdew), and wind speed from the ERA5 dataset were considered. All the sub-daily satellite/reanalysis data covering 2002-2019 were aggregated into a daily scale.
As the boundary of the upper and mid-lower reaches of the Hanjiang River Basin, the short-series inflow of the Danjiangkou Reservoir was selected to calibrate the parameters of the hydrological model. Observed streamflow data spanning 2003-2007 were obtained from the Yangtze River Water Conservancy Commission.

Methodology
A flowchart of the framework module is presented in Figure 3, which is further elaborated in the following sections. It is worth mentioning that our proposed framework can be used in hydrological simulation rather than forecasting.

Hydrological Model
The modèle du Génie Rural à 9 paramètres Journalier (GR4J-9) hydrological model was used to initially simulate the hydrology of the upper Hanjiang River watershed. The GR4J-9 model is a daily lumped nine-parameter rainfall-runoff model that integrates a traditional GR4J (five-parameter version) hydrological model with the CemaNeige snowfall accumulation and snowmelt module (which occupies four parameters). The GR4J-9 model belongs to the family of soil moisture accounting models that route runoff through two interconnected reservoirs (i.e., production and routing reservoirs) and two-unit hydrographs. It has several main parameters: the maximum capacity of the production reservoir, the groundwater exchange coefficient, the 1-day maximum retention capacity of

Hydrological Model
The modèle du Génie Rural à 9 paramètres Journalier (GR4J-9) hydrological model was used to initially simulate the hydrology of the upper Hanjiang River watershed. The GR4J-9 model is a daily lumped nine-parameter rainfall-runoff model that integrates a traditional GR4J (five-parameter version) hydrological model with the CemaNeige snowfall accumulation and snowmelt module (which occupies four parameters). The GR4J-9 model belongs to the family of soil moisture accounting models that route runoff through two interconnected reservoirs (i.e., production and routing reservoirs) and two-unit hydrographs. It has several main parameters: the maximum capacity of the production reservoir, the groundwater exchange coefficient, the 1-day maximum retention capacity of the routing reservoir, and the time base of the unit hydrograph. This model has been tested in a large sample of catchments and has shown competitive performance compared to more complex models with more parameters [26,27]. Yang et al. [28] showed that the performance of GR4J is more stable than other models (i.e., WASMOD, HBV, and XAJ) in a changing climate. They also set the fixed coefficient of percolation leakage as a free parameter to better fit the study area and calibrated it for the objective watershed. The potential evaporation in the GR4J-9 model is obtained from the temperature-based Oudin method [29].
The observed daily precipitation (P), maximum and minimum temperature (T max and T min, respectively) and flow discharge were fed to calibrate and validate the GR4J-9 model for the experimental watershed. We optimized the parameters of the hydrological model using the Shuffled Complex Evolution (SCE-UA) method developed at the University of Arizona [30]. The SCE-UA method integrates the advantages of several effective global optimization concepts and employs both deterministic search strategies and random schemes to achieve an effective search ability.

LSTM Model
An LSTM model is a specific kind of RNN designed to overcome the drawbacks caused by a vanishing gradient or exploding in the process of training the RNN using backpropagation through time (BPTT) [31,32]. It sets up a dedicated memory cell that stores information over long periods, potentially making it an ideal candidate for modeling dynamic systems such as watersheds. An unfolded computational graph, as depicted in Figure 4, can reveal the working principle of the LSTM method. One LSTM unit is composed of an input gate, a forget gate, a memory cell and an output gate. The input gate decides which new value regulated by the memory cell will be updated in the cell state, and the forget gate controls the information to remove or retain in the cell state. A general memory block of an LSTM structure can be described by the following equations: where C(t + 1) and C(t) are the cell state at time t + 1 and t, respectively; X(t + 1) and H(t + 1) are the network input and the recurrent input at time t + 1, respectively. At the initial time step, both the cell and hidden states are initialized as a vector of zeros. w and W are the weights of the link between gates and layers, respectively to more complex models with more parameters [26,27]. Yang et al. [28] showed that the performance of GR4J is more stable than other models (i.e., WASMOD, HBV, and XAJ) in a changing climate. They also set the fixed coefficient of percolation leakage as a free parameter to better fit the study area and calibrated it for the objective watershed. The potential evaporation in the GR4J-9 model is obtained from the temperature-based Oudin method [29]. The observed daily precipitation (P), maximum and minimum temperature (Tmax and Tmin, respectively) and flow discharge were fed to calibrate and validate the GR4J-9 model for the experimental watershed. We optimized the parameters of the hydrological model using the Shuffled Complex Evolution (SCE-UA) method developed at the University of Arizona [30]. The SCE-UA method integrates the advantages of several effective global optimization concepts and employs both deterministic search strategies and random schemes to achieve an effective search ability.

LSTM Model
An LSTM model is a specific kind of RNN designed to overcome the drawbacks caused by a vanishing gradient or exploding in the process of training the RNN using backpropagation through time (BPTT) [31,32]. It sets up a dedicated memory cell that stores information over long periods, potentially making it an ideal candidate for modeling dynamic systems such as watersheds. An unfolded computational graph, as depicted in Figure 4, can reveal the working principle of the LSTM method. One LSTM unit is composed of an input gate, a forget gate, a memory cell and an output gate. The input gate decides which new value regulated by the memory cell will be updated in the cell state, and the forget gate controls the information to remove or retain in the cell state. A general memory block of an LSTM structure can be described by the following equations: where C(t + 1) and C(t) are the cell state at time t + 1 and t, respectively; X(t + 1) and H(t + 1) are the network input and the recurrent input at time t + 1, respectively. At the initial time step, both the cell and hidden states are initialized as a vector of zeros. w and W are the weights of the link between gates and layers, respectively; bi, bf, bc, and bo are learnable bias parameters for each gate; [ ] σ ⋅ is the sigmoid function and tanh[ ] ⋅ is the hyperbolic tangent function; both are activation functions with objective values ranging from 0 to 1. In this study, the two hyperparameters (the initial learning rate and the number of hidden nodes) diversely affecting LSTM model performance needed to be determined. A larger number of hidden neurons lead to a fully trained model, which may overfit the  In this study, the two hyperparameters (the initial learning rate and the number of hidden nodes) diversely affecting LSTM model performance needed to be determined. A larger number of hidden neurons lead to a fully trained model, which may overfit the data; conversely, a small number of hidden neurons may cause randomness with high bias. A more detailed process of hyper-parameter tuning is described in Section 4.2. For simplicity, we chose a three-layer LSTM network as the fully connected structure, which consists of one input layer, one hidden layer, and one output layer. Preliminary investigations found that LSTMs with one single hidden layer are capable of simulating streamflow in Hanjiang River Basin [33]. The BPTT algorithm [34] was used to train the LSTM model and the adaptive moment estimation (ADAM) algorithm [35] was employed as the learning rate method. Finally, the mean square error was treated as the loss index. The general model implementation can be accessed from the Statistics and Machine Learning Toolbox of the MATLAB software (website: https://ww2.mathworks.cn/products/statistics.html# machine-learning, accessed on 22 March 2021) [36].

Input Variable Selection (IVS)
The appropriate determination of inputs substantially influences the development of the LSTM model [37]. For streamflow simulation, a dataset of candidate inputs typically covers observed predictors (e.g., initial basin conditions or climate elements) as well as lagged streamflow observations. A wide array of potentially hydrological components can be fed into the model; however, many may only add redundancy or a high level of noise into the model. Furthermore, some candidate inputs (e.g., long-lagged temporal data) may add little or no value to the rainfall-runoff model. To this end, the idea of input variable selection (IVS) [38] is introduced.
A tree-based IVS method developed by Galelli and Castelletti [39] was implemented to identify the optimal input combination. The extremely random tree (extra-tree) method is a non-parametric tree regression approach that partitions the input space into mutually exclusive regions on a predefined principle of splitting nodes [40]. In this particular structure, the extra-tree can rank the importance of the input variables by scoring each input variable by evaluation of the relative variance reduction. It adopts a goodness-of-fit criterion of the coefficient of determination (R 2 ) to systematically select the most significant and non-redundant input space, which was found to consistently indicate the optimum LSTM structure in water resource modeling applications [41].
We used basin-averaged daily mean air temperature, precipitation, wind speed, relative humidity (RH), and the simulated daily discharge (corresponding to the observed reservoir inflow) from the simulations of the GR4J-9 model as the candidate inputs of the LSTM, and used observed daily discharge as the response output. Considering a typical e-folding time scale (recession time) of streamflow, we set the time lag at 4 days.

Simulation Performance Assessment
To ensure that the trained simulation model did not contain known or detectable defects and could be used on any unseen data, a comprehensive assessment was performed considering different aspects of the modeled simulation flow. Specifically, the Kling-Gupta efficiency (KGE) index was selected to describe the statistical accuracy of the hydrological models, and the objective function was to maximize the KGE value during the calibration period.
where r refers to Pearson's linear correlation coefficient between the observation and the simulations, and α (β) indicates the ratio of standard deviations (mean value) of the observed and simulated streamflow. KGE varies (−∞, 1]; a value closer to 1 represents a better simulation. As shown by previous studies [42][43][44], another metric (mean relative absolute error (MRAE)) can be coupled with KGE to evaluate the overall deterministic performance.

Operation Model
As mentioned in Section 2.1, water supply and power generation are the two main yet conflicting objectives of the Danjiangkou Reservoir. Their mathematical formulations can be expressed by Equations (5) and (6).
where W and E are water supply yield (m 3 ) and power generation (kW·h) per year, respectively; N t is the power output at time t (kW); Q d t and Q P t are water diversion flow and release discharge for power generation at time t (m 3 /s), respectively; k is hydropower generation efficiency; H t is the average water head at time t (m); ∆t is the time step (s); and T is the total number of operational periods. The reservoir operation model obeys some physical constraints, which were outlined by He et al. [45]. The mathematical equations of these constraints are omitted for the sake of brevity.

Operating Strategy of Reservoir Release
The optimal reservoir operation determines the reservoir release sequence Q out t during the whole operating period for the maximization of W and E. As the optimization strategy of reservoir release involves a high-dimensional and non-linear property, Gaussian radial bias functions (RBFs) are taken as the operating policy, since they are flexible to make decisions with strong universal approximation [46,47]. In the RBFs method, Q out t can be expressed in Equations (7) and (8).
where U is the total number of RBFs ϕ(·); ω u is the weight of the uth RBF, the sum of all weights is 1, e.g., M is the number of input variables of X t ; and c m,u and b u are the mth-dimensional center and radius of the uth RBF, respectively. For an individual reservoir, X t usually consists of three variables, namely time at time t, current reservoir storage (V t ) and reservoir inflow information (Q in t ) [48]; thus, M was set to 3. Moreover, each RBF can be regarded as one pattern of decision-making in reservoir operation based on X t and, ultimately, decision-making is determined by the combination of four patterns (i.e., U is 4) as suggested by Yang et al. [48].
Consequently, there were 20 parameters to be calibrated for the sum of the RBFs. We optimized the parameter combination based on the parameterization-simulationoptimization (PSO) framework using the non-dominated sorting genetic algorithm II (NSGA-II). To converge the Pareto front, the evolutionary NSGA-II algorithm adopts the fast non-dominated sorting and crowding distance strategies. The experimental setup of NSGA-II was: population size = 100, generation number = 1000, crossover probability = 0.9, and mutation probability = 0.1.

Initial Simulation of the GR4J-9 Model
We first calibrated and validated the GR4J-9 model using observed reservoir inflow during 2003-2007. With the first year serving as the warm-up period, the KGE values of the calibration (2004-2006, three years) and validation periods (2007, one year) were 0.76 and 0.54, respectively; the MRAE values were 0.56 and 0.73, respectively. As recommended by previous studies, the model performance is judged to be satisfactory for flow simulations if the daily KGE is greater than 0.5 and the MRAE is less than 0.85 for watershed-scale models [2,49]. Figure 5 depicts the simulation results for reservoir inflow. It shows that the GR4J-9 model could fit the inflow hydrograph of the calibration well, especially for a low inflow regime. However, it achieved a relatively low KGE value in the validation period, with a serious underestimation. In detail, the GR4J-9 model cannot capture the peak discharge of 28,900 m 3 /s, and instead, gives a lower value of 12,461 m 3 /s. Similar results were also found at other different times, which were caused by several aspects. From the view of model inputs, the basin-averaged daily precipitation still has a relatively low resolution compared to actual precipitation observations. From the view of model structure, GR4J-9 has a simple structure with nine model parameters. Although it has superior performance compared to distributed models in the ungauged basin (the latter need more sub-basin observations to calibrate parameters), it inevitably leads to a model error where the simple structure assumption fails to cater to the actual hydrological condition. Therefore, it is necessary to develop a state-of-the-art technique to further improve streamflow accuracy.

Initial Simulation of the GR4J-9 Model
We first calibrated and validated the GR4J-9 model using observed reservoir inflow during 2003-2007. With the first year serving as the warm-up period, the KGE values of the calibration (2004-2006, three years) and validation periods (2007, one year) were 0.76 and 0.54, respectively; the MRAE values were 0.56 and 0.73, respectively. As recommended by previous studies, the model performance is judged to be satisfactory for flow simulations if the daily KGE is greater than 0.5 and the MRAE is less than 0.85 for watershed-scale models [2,49]. Figure 5 depicts the simulation results for reservoir inflow. It shows that the GR4J-9 model could fit the inflow hydrograph of the calibration well, especially for a low inflow regime. However, it achieved a relatively low KGE value in the validation period, with a serious underestimation. In detail, the GR4J-9 model cannot capture the peak discharge of 28,900 m 3 /s, and instead, gives a lower value of 12,461 m 3 /s. Similar results were also found at other different times, which were caused by several aspects. From the view of model inputs, the basin-averaged daily precipitation still has a relatively low resolution compared to actual precipitation observations. From the view of model structure, GR4J-9 has a simple structure with nine model parameters. Although it has superior performance compared to distributed models in the ungauged basin (the latter need more sub-basin observations to calibrate parameters), it inevitably leads to a model error where the simple structure assumption fails to cater to the actual hydrological condition. Therefore, it is necessary to develop a state-of-the-art technique to further improve streamflow accuracy.

LSTM Performance
As stated in Section 3.2.2, we fed the GR4J-9 model output into the advanced LSTM model, which also included wind speed and RH. To derive RH data for LSTM inputs, we used the daily dew point temperature (Tdew) and daily mean temperature (Tmean) from ERA5. The actual vapor pressure (e) and saturated vapor pressure ( ) were derived using the Clausius-Clapeyron equation. The determined input variables were normalized to eliminate the influence of magnitude, thereby improving the accuracy and efficiency of network learning.

LSTM Performance
As stated in Section 3.2.2, we fed the GR4J-9 model output into the advanced LSTM model, which also included wind speed and RH. To derive RH data for LSTM inputs, we used the daily dew point temperature (T dew ) and daily mean temperature (T mean ) from ERA5. The actual vapor pressure (e) and saturated vapor pressure (e sa ) were derived using the Clausius-Clapeyron equation. The determined input variables were normalized to eliminate the influence of magnitude, thereby improving the accuracy and efficiency of network learning.
As anticipated, the tree-based IVS algorithm could score and rank input variables in terms of their relevance to the output. To ensure a reliable ranking result, the experiment was cross-validated multiple times with different shuffled datasets, and the inputs were sorted in decreasing order. The score results of the IVS run are presented in Table 2, which illustrates the importance of each selected variable. This tree-based method can make sense for the ex-post physical interpretation of the cause-effect relationships captured by the model. For the upper Hanjiang River Basin, the simulated flow at time t by the GR4J-9 model (Q sim t ), precipitation (P t ), antecedent simulated flow and precipitation with 1-or 2-time lag (Q sim t−1 , Q sim t−2 and P t−1 , respectively) were the top five most important variables (about 68% of the ensemble total score), followed by relative humidity and wind speed at time t (i.e., RH t and WD t , respectively). Except for the simulated streamflow element, which was highly related to the observed inflow, we found that P t−1 and P t ranked in the top positions, with relative scores of 13% and 6%, respectively. This may be due to the hydraulic characteristics of this large catchment, which is drained by base flow with a long period of concentration. The basin-averaged RH and wind speed are less important, but not negligible. Finally, antecedent simulated flow with 1-4 time lags (Q sim t ,Q sim t−1 ,Q sim t−2 ,Q sim t−3 , and Q sim t−4 , respectively), antecedent simulated flow with 1-3 time lags (P t ,P t−1 ,P t−2 and P t−3 , respectively), RH t and WD t consisted of 11 inputs for the LSTM machine learning model. With a hidden layer of 64 neurons and an initial learning rate of 0.1 (identified by the trialand-error method), the LSTM-based model could substantially improve the accuracy of the streamflow simulation compared to the GR4J-9 benchmark model. The daily streamflow trajectories are shown in Figure 6a. The model achieved a high KGE value of 0.87 and MRAE of 0.56 for the daily discharge in the calibration period; additionally, the KGE was 0.68 and the MRAE value was 0.71 in the validation period. The results demonstrated that this method can competently ameliorate the hydrological data scarcity. Compared to the benchmark GR4J-9 model, the LSTM model uses related hydrological variables (i.e., RH and wind speed) as driving inputs to improve streamflow accuracy. Compared to physically distributed hydrological models, which have limited applications in basins with shortseries datasets due to the complex characteristics [50], the LSTM model can sufficiently use remote sensing data and has the features of easy-to-use and highly efficient. However, the LSTM model displayed serious overestimation behavior in the validation period. This is due to the LSTM model being overly reliant on the calibrated data without exception. Figure 6a,b show that the streamflow simulations corrected by the LSTM model are closer to the peak discharge in the calibration. This overfitting performance was unfortunately carried into the validation period, which caused the LSTM model to provide a relatively high value for streamflow discharge under low flow regimes. In general, the LSTM model has room for improvement, but this does not hinder its application value.

Potential Application in Reservoir Management
We acquired a long-series daily streamflow simulation from the period of 2008-2019 by feeding the remote sensing data into the calibrated LSTM model. We could formulate more scientific strategies for basins with hydrological data scarcity. Using the mediumand long-term management of the Danjiangkou Reservoir operation as an example, we aimed to improve hydropower benefits and water supply yield and balance them as much as possible. Before performing the NSGA-II reservoir optimization method, daily scale simulated streamflow was converted into 10-day average runoff.
The optimization results of W and E by NSGA-II are presented in Figure 7, and the operation results merely based on short-series observation (2003)(2004)(2005)(2006)(2007) are also provided for further comparison. Both of the Pareto fronts under different scenarios are widely and evenly distributed between the two conflicting objectives. Considering the pursuit of maximum economic profit (official electricity price: 0.21 RMB/kWh; water price for the MSWTP project: 0.13 RMB/m 3 ), the final two optimal operating rules were chosen for

Potential Application in Reservoir Management
We acquired a long-series daily streamflow simulation from the period of 2008-2019 by feeding the remote sensing data into the calibrated LSTM model. We could formulate more scientific strategies for basins with hydrological data scarcity. Using the mediumand long-term management of the Danjiangkou Reservoir operation as an example, we aimed to improve hydropower benefits and water supply yield and balance them as much as possible. Before performing the NSGA-II reservoir optimization method, daily scale simulated streamflow was converted into 10-day average runoff.
The optimization results of W and E by NSGA-II are presented in Figure 7, and the operation results merely based on short-series observation (2003)(2004)(2005)(2006)(2007)   With these two optimal operating rules guiding reservoir operation during the period of 2003-2019, their objective results are summarized in Table 3. It can be inferred from water supply, hydropower, and economic profit (W, E, and H, respectively) that Solution I prefers a higher value of W than Solution II, no matter whether in the observation period or in the whole period when decision-makers want to pursue more economic profit. As shown in Table 3, while the performance of Solution I (obtained through LSTM-based streamflow simulation) is similar to that of Solution II, it provides decision-makers with another viable operating way that requires more real observational data to verify.

Conclusions
Hydro-meteorological data scarcity impairs hydrological simulation, manifesting a pressing need to develop an alternative scheme in this field. Satellite-based and atmospheric reanalysis estimation may provide feasible access to reproduce the hydrological recycle process and may have potential value. To this end, this paper proposed a novel method to integrate open-source remote sensing data and ML-based inversion techniques for hydrology. Furthermore, we applied it in reservoir management. According to the results, we reached the following conclusions: (1) Driven by the synthetic data generated by a lumped hydrological GR4J model, satellite-based data and ERA5 could overcome the limitation of historical observation With these two optimal operating rules guiding reservoir operation during the period of 2003-2019, their objective results are summarized in Table 3. It can be inferred from water supply, hydropower, and economic profit (W, E, and H, respectively) that Solution I prefers a higher value of W than Solution II, no matter whether in the observation period or in the whole period when decision-makers want to pursue more economic profit. As shown in Table 3, while the performance of Solution I (obtained through LSTM-based streamflow simulation) is similar to that of Solution II, it provides decision-makers with another viable operating way that requires more real observational data to verify.

Conclusions
Hydro-meteorological data scarcity impairs hydrological simulation, manifesting a pressing need to develop an alternative scheme in this field. Satellite-based and atmospheric reanalysis estimation may provide feasible access to reproduce the hydrological recycle process and may have potential value. To this end, this paper proposed a novel method to integrate open-source remote sensing data and ML-based inversion techniques for hydrology. Furthermore, we applied it in reservoir management. According to the results, we reached the following conclusions: (1) Driven by the synthetic data generated by a lumped hydrological GR4J model, satellite-based data and ERA5 could overcome the limitation of historical observation scarcity. With a KGE value of 0.54 in the validation period for the upper Hanjiang River Basin, the traditional lumped hydrological model can be applied to the ungauged basins but there is still room for improvement. (2) Compared to the traditional GR4J model, the ML-based data-driven model showed its superior performance in capturing the long-series time sequence using a sophisticated network structure. Along with inheriting the simulation output of the traditional hydrological model, the LSTM model can further mine the value of remote sensing data related to hydrological variables. Compared to traditional distributed hydrological models, its model structure is simple and can be highly efficient. (3) The LSTM-based streamflow simulation scenario can provide the basis for another scientific operation way for reservoir managers, which requires future validation of the potential value of the LSTM-based method.
Despite the outstanding performance of the developed methodology, some work remains for further exploration. First, hydrological models have different behaviors depending on the flow regimes. However, in this study, the hydro-meteorological data with a one-day timescale were fed to drive one set of hydrological models, yet the separation of flood seasons and non-flood seasons was neglected. Besides, a simpler LSTM model should be taken into consideration and compared with our proposed hybrid model for hydrological performance. Secondly, the methodology was merely applied for hydrological simulation rather than hydrological forecasts. Some products such as the Global Ensemble Forecast System (GEFS) Reforecast can be included to improve this methodology. In the future, we will explore the ML-based method with remote sensing data to verify its generalizability in more ungauged basins.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.