Mekong Delta Runoff Prediction Using Standardized Remotely-Sensed Water Balance Variables

A suitable routing model for predicting future monthly water discharge (WD) is essential for operational hydrology, including water supply, and hydrological extreme management, to mention but a few. This is particularly important for a remote area without a sufficient number of in-situ data, promoting the usage of remotely sensed surface variables. Direct correlation analysis between ground-observed WD and localized passive remotely-sensed surface variables (e.g., indices and geometric variables) has been studied extensively over the past two decades. Most of these related studies focused on the usage of constructed correlative relationships for estimating WD at ungauged locations. Nevertheless, temporal prediction performance of monthly runoff (R) (being an average representation of WD of a catchment) at the river delta reconstructed from the basin’s upstream remotely-sensed water balance variables via a standardization approach has not been explored. This study examined the standardization approach via linear regression using the remotely-sensed water balance variables from upstream of the Mekong Basin to reconstruct and predict monthly R time series at the Mekong Delta. This was subsequently compared to that based on artificial intelligence (AI) models. Accounting for less than 1% improvement via the AI-based models over that of a direct linear regression, our results showed that both the reconstructed and predicted Rs based on the proposed approach yielded a 2–6% further improvement, in particular the reduction of discrepancy in the peak and trough of WD, over those reconstructed and predicted from the remotely-sensed water balance variables without standardization. This further indicated the advantage of the proposed standardization approach to mitigate potential environmental influences. The best R, predicted from standardized water storage over the whole upstream area, attained the highest Pearson correlation coefficient of 0.978 and Nash–Sutcliffe efficiency of 0.947, and the lowest normalized root-mean-square error of 0.072.


Introduction
Owing to a lack of funding for facility maintenance [1], there has been a continuous decline in the number of ground-observed hydrological sites for observing water discharge (WD) worldwide since the 1970s [2]. Therefore, remote sensing (RS) has been advocated because of its global coverage. Normalized difference vegetation index, river width, floodplain inundation, and land surface temperature are common land surface variables derived from passive RS [3][4][5][6]. Regardless of their relationships with ground-observed WD, the passive RS surface variables are directly correlated with water level or WD as long as they are well correlated with each other in the RS field of study. However, this is not a (i.e., MD) that causes agricultural and economic losses (e.g., [37][38][39]). This yields the concept of using upstream WBVs for downstream WD prediction. While the downstream R has been reconstructed and estimated by using upstream RS-derived surface variables (e.g., [40]), satellite radar altimetric water level (e.g., [8,14,41]), and gravimetrically-determined TWS [16], the temporal prediction of the estuarine R at the downstream using the upstream RS WBVs has not been explored.
the concept of using upstream WBVs for downstream WD prediction. While the downstream R has been reconstructed and estimated by using upstream RS-derived surface variables (e.g., [40]), satellite radar altimetric water level (e.g., [8,14,41]), and gravimetrically-determined TWS [16], the temporal prediction of the estuarine R at the downstream using the upstream RS WBVs has not been explored. In addition, many dams at the upstream MB were constructed during the 1990-2007 period [42]. This significantly alters the WD at the downstream MB in different seasons compared to that before dam construction [43,44]. However, no significant annual change of WD in the MD has been found [42,44]. Since the WD has already been altered in different seasons at the downstream MD, the dam operation effect would accumulate for each season. In case it is accumulated, in the form of a partial bias, this effect can be partly mitigated by the data standardization approach in this study.
Given the above justification, this study aims to examine temporal prediction performance of monthly R at the MD reconstructed from the upstream standardized RS WBVs employing the data standardization approach. The R time series is predicted based on the reconstructed relationship, followed by comparing to the data time period independent of the period used for R reconstruction. In addition, many dams at the upstream MB were constructed during the 1990-2007 period [42]. This significantly alters the WD at the downstream MB in different seasons compared to that before dam construction [43,44]. However, no significant annual change of WD in the MD has been found [42,44]. Since the WD has already been altered in different seasons at the downstream MD, the dam operation effect would accumulate for each season. In case it is accumulated, in the form of a partial bias, this effect can be partly mitigated by the data standardization approach in this study.
Given the above justification, this study aims to examine temporal prediction performance of monthly R at the MD reconstructed from the upstream standardized RS WBVs employing the data standardization approach. The R time series is predicted based on the reconstructed relationship, followed by comparing to the data time period independent of the period used for R reconstruction. To show the effectiveness of our proposed approach, it is subsequently compared with those predictions from AI-based models.

Ground-Observed Discharge and Passive RS Surface Variables
Within the MD, Tan Chau, Chau Doc, Can Tho, and My Thuan are the four ground-observed discharge stations (Figure 1). Note that the Tan Chau and Chau Doc station pair (hereafter denoted as TC-CD pair), and Can Tho and My Thuan station pair (hereafter denoted as MT-CT pair) arẽ 220 km and 70-100 km away from the estuary mouth, respectively. The latter pair is subject to a more significant ocean tidal backwater effect than that of the former pair. To reduce semi-diurnal and diurnal ocean tidal effects from the South China Sea, the WD time series of TC-CD pair were summed up. The same procedure was applied to the MT-CT pair. Two pairs of stations were used to assess the remaining ocean tidal backwater effect.
The above WD station time series were available on request at http://www.mrcmekong.org. Since most dams were constructed before 2009, the ground-observed WD time series during 2009-2012 and 2013-2014 were chosen for R reconstruction and prediction assessment, respectively. Note that the WD (in m 3 /s) were daily sampled. This had to be converted into the R at a monthly rate (in mm/mon), which was calculated by adding up the daily WD per month and subsequently dividing by the approximate MB surface area (i.e., 795,000 km 2 ). To show the effectiveness of our proposed approach, it is subsequently compared with those predictions from AI-based models.

Ground-Observed Discharge and Passive RS Surface Variables
Within the MD, Tan Chau, Chau Doc, Can Tho, and My Thuan are the four ground-observed discharge stations (Figure 1). Note that the Tan Chau and Chau Doc station pair (hereafter denoted as TC-CD pair), and Can Tho and My Thuan station pair (hereafter denoted as MT-CT pair) are ~220 km and 70-100 km away from the estuary mouth, respectively. The latter pair is subject to a more significant ocean tidal backwater effect than that of the former pair. To reduce semi-diurnal and diurnal ocean tidal effects from the South China Sea, the WD time series of TC-CD pair were summed up. The same procedure was applied to the MT-CT pair. Two pairs of stations were used to assess the remaining ocean tidal backwater effect.
The above WD station time series were available on request at http://www.mrcmekong.org. Since most dams were constructed before 2009, the ground-observed WD time series during 2009-2012 and 2013-2014 were chosen for R reconstruction and prediction assessment, respectively. Note that the WD (in m 3 /s) were daily sampled. This had to be converted into the R at a monthly rate (in mm/mon), which was calculated by adding up the daily WD per month and subsequently dividing by the approximate MB surface area (i.e., 795,000 km 2 ).  The land surface temperature (LST) (i.e., MOD11C3) and normalized difference vegetation index (NDVI) (i.e., MOD13C2), which can be downloaded at https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table, were two of the passive RS surface data products obtained from MODIS. These data directly served as baseline results for the performance assessment against the standardized RS WBVs.

Water Balance Variables and Their Standardization Based on RS Data
The three RS WBVs (S, ET, and P) were measured or inferred from GRACE, MODIS, and TRMM, respectively. To avoid the dams constructed before 2009, the monthly RS WBVs during The land surface temperature (LST) (i.e., MOD11C3) and normalized difference vegetation index (NDVI) (i.e., MOD13C2), which can be downloaded at https://lpdaac.usgs.gov/dataset_discovery/ modis/modis_products_table, were two of the passive RS surface data products obtained from MODIS. These data directly served as baseline results for the performance assessment against the standardized RS WBVs.

Water Balance Variables and Their Standardization Based on RS Data
The three RS WBVs (S, ET, and P) were measured or inferred from GRACE, MODIS, and TRMM, respectively. To avoid the dams constructed before 2009, the monthly RS WBVs during 2009-2012 and 2013-2014 were extracted for R reconstruction and prediction assessment, respectively. All passive RS surface variables and RS WBVs are visualized in Figure 3.

GRACE Terrestrial Water Storage (S) and its Standardization
Given time-variable gravity observations from GRACE, globally gridded S time series can be inferred by Equation (14) in [45] divided by the averaged density of water. The Center for Space Research (CSR) Release (RL) 06 solution (hereinafter called CSR RL06) were utilized. The monthly RL06 data were in the form of Stokes coefficients (SCs) that represent the relative mass changes expanded up to degree 60. These data were available at http://icgem.gfz-potsdam.de/series. Post-processing steps were required before usage, which included the addition of the degree-1 and

GRACE Terrestrial Water Storage (S) and Its Standardization
Given time-variable gravity observations from GRACE, globally gridded S time series can be inferred by Equation (14) in [45] divided by the averaged density of water. The Center for Space Research (CSR) Release (RL) 06 solution (hereinafter called CSR RL06) were utilized. The monthly RL06 data were in the form of Stokes coefficients (SCs) that represent the relative mass changes expanded up to degree 60. These data were available at http://icgem.gfz-potsdam.de/series. Post-processing steps were required before usage, which included the addition of the degree-1 and the replacement of the C 20 term in SCs for restoring the geocenter motion and geoid components [46,47], de-striping process, and Gaussian filtering with a 350-km radius for minimizing correlated errors in space [48,49].
The inferred monthly S time series were then utilized to correlate with the ground-observed R and to calculate its standardization (hereinafter abbreviated as SI), followed by correlating with the standardized R. The SI can be calculated as: where S j,k , median(S k ), and s k are the observed S, and the median and sampled standard deviation of GRACE S for year j and month k, respectively. The same procedure was applied to ET when calculating its corresponding standardization.

MODIS Evapotranspiration (ET) and Its Standardization
Given the observed environmental variables from MODIS, global inferred ET can be inferred by an improved algorithm in [50]. The ET data (MOD16A2), with 0.5 • gridded resolution, were utilized and available from NASA Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/ dataset_discovery/modis/modis_products_table).
Similar to the above, ET values can be used to calculate its standardized form (hereinafter abbreviated as SETI) using Equation (1). Despite long-term continuous ET time series beingdesirable for climatic study, the ET time series spanning from 2000 to 2014 were sufficient in the context of water balance. Consequently, SETI during 2000-2014 were calculated, whereas the calculated standardized values between 2009and 2014 were utilized in this study.

TRMM Precipitation (P) and Its Standardization
Given the measured P from TRMM [51], the monthly P data (TRMM 3B43 version 7), with a gridded resolution of 0.25 • , were utilized and available from the NASA Goddard Space Flight Center at https://disc.gsfc.nasa.gov/datasets/TRMM_3B43_V7/summary. Notice that these data products were dependent on the Global Precipitation Climatology Centre (GPCP) rain gauge data, as assimilation and calibration processes were conducted (e.g., [52]).
Since P does not occur every day, its standardization (hereinafter denoted as standardized precipitation index (SPI) [53]), which is different from the above Equation (1), involves calculations of the standardization of cumulative probability distribution [54]. Note that equations are not shown as it is a well-known meteorological drought index. For a similar reason to that of ET, SPI during 1998-2014 were calculated, while the standardized values between 2009 and 2014 were extracted for usage in this study.

Linear Regression Analysis and Data Standardization
The processed RS WBVs and their standardizations were averaged over Yunnan Province, China and the whole upstream area (from Tibetan Plateau to southern boundary of Cambodia) ( Figure 1). Before linear regression analysis, a five-month moving averaging process was conducted for further smoothing. Note that one-month (between R and NDVI and LST) and two-month backward time shift (between R and TRMM P and MODIS ET) were applied in order to maximize the correlation due to inherent hysteretic properties of a river basin within a hydrological cycle [21]. Figure 4 displays an initial insight into the time series pattern of the RS WBVs against the ground-observed R during 2009-2012, showing that the general temporal patterns of all RS WBVs were largely consistent with the ground-observed R. Note that the passive RS surface variables (i.e., LST and NDVI) are not displayed because of the presence of large offsets (Table 1). A variable time lag within a year was also noticed for all RS WBVs at the upstream MB against the ground-observed R at the MD. We speculate that different precipitation patterns between the upstream MB and the MD every year might be caused by El Niño-Southern Oscillation (ENSO) variability [55], and hence the adaptive changes in evapotranspiration and terrestrial water storage. After all, the MB spans across 25 • of latitude.
Water 2020, 12, x 7 of 17 MD every year might be caused by El Niño-Southern Oscillation (ENSO) variability [55], and hence the adaptive changes in evapotranspiration and terrestrial water storage. After all, the MB spans across 25° of latitude. For a direct correlation between RS-observed surface variables and ground-observed R in a traditional practice, a linear regression analysis is conducted to fit with two parameters, which is expressed as: where t y and δ The two estimated parameters and their corresponding standard deviation of Equation (2) for the RS surface variables are shown in Tables 1 and 2. Though a large offset and its standard deviation were presented for the ground-observed R fitted with LST, the fitted relationship should have been as well as other RS surface variables because the relative errors among them were almost the same, whether they were averaged over Yunnan Province or the whole upstream area (Table 1 and 2). In addition, all estimated parameters were significantly different from zero in terms of statistics. These two determined parameters were then directly employed to reconstruct R time series using the same values of RS surface variables for data fitting during 2009-2012, while R time  For a direct correlation between RS-observed surface variables and ground-observed R in a traditional practice, a linear regression analysis is conducted to fit with two parameters, which is expressed as: where y t and x t−δ are the ground-observed R (either at the MT-CT or TC-CD station pair) at a monthly time t, and each RS WBVs at a monthly time t − δ; δ is the time shift applied to RS-observed surface variables (i.e., LST, NDVI, TRMM P, and MODIS ET), as mentioned; a and b are the two parameters to be estimated via the linear regression analysis using the data time period between 2009 and 2012. The two estimated parameters and their corresponding standard deviation of Equation (2) for the RS surface variables are shown in Tables 1 and 2. Though a large offset and its standard deviation were presented for the ground-observed R fitted with LST, the fitted relationship should have been as well as other RS surface variables because the relative errors among them were almost the same, whether they were averaged over Yunnan Province or the whole upstream area (Tables 1 and 2). In addition, all estimated parameters were significantly different from zero in terms of statistics. These two determined parameters were then directly employed to reconstruct R time series using the same values of RS surface variables for data fitting during 2009-2012, while R time series were predicted by the two determined parameters using the values of RS surface variables during 2013-2014. For our data standardization approach, the ground-observed R time series were standardized.The standardized ground-observed R, sR j,k , was achieved by: where R j,k , med(R k ), and s k are the ground-observed R, the median, and sampled standard deviation of ground-observed R for year j and month k, respectively. The standardized ground-observed R was then correlated with the standardized P (i.e., SPI), ET (i.e., SETI), and S (i.e., SI), as mentioned in Section 2.2, via Equation (2), to determine the corresponding a and b. In this case, these two determined parameters were then employed to calculate the corresponding standardized R using the same values of RS WBVs for data fitting during 2009-2012. The calculated corresponding standardized R generated from RS WBVs, together with the previously calculated med(R k ) and s k ,were then used to generate the reconstructed R time series via Equation (3), while R time series were predicted in the same manner using RS WBVs during 2013-2014.
After the above procedures, the reconstructed and predicted R time series were then compared against the respective ground-observed R to assess the reconstruction and prediction performance during different periods, as mentioned.

AI-Based Models
Various neural network models have been used in the hydrology field, including streamflow forecast, flood prediction, groundwater quality estimation, water resource management, and so on. Artificial neural network (ANN) and long-short term memory (LSTM) models were employed in this study. The ANN model includes an input layer, hidden layer, and output layer. In this study, the predictors were NDVI, LST, P, ET, and S. R was placed at the output layer as an output. The hidden layers consisted of two layers, and each layer contained four neurons. A scaled exponential linear unit (SELU) served as an activation function of each layer. For the hyperparameter settings, the loss function was mean absolute error. The Adam's optimizer was employed with the epoch set to 100 times. The callback function "ReduceLROnPlateau" within Keras was applied to adjust the learning rate automatically. A tenth of the training data was divided into validation data sets.
LSTM, proposed based on recurrent neural networks (RNNs), is one of the most powerful type of neural networks and is capable of processing sequences of arbitrary input patterns to overcome the gradient problem existing in traditional RNNs. In this study, the construction of the neural network formed by LSTM cells was similar to ANN; two hidden layers and the same activation function were employed. Due to non-unique results during the training stage., the two models were run ten times each. The means of the valuation parameters obtained from each training were the final valuation parameters.

Assessment Indicators
To examine the reconstructed and predicted R time series against the ground-observed R 0 , the Pearson correlation coefficient (PCC), normalized root-mean-square error (NRMSE), and Nash-Sutcliffe efficiency (NSE) model coefficient were used.
PCC measures the collinearity between two datasets. It ranges from −1 to 1, which is calculated as NRMSE measures the RMSE in terms of a fraction of error (e.g., [56]), where RMSE does not indicate the amount of relative error with respect to the observed quantities. NRMSE is defined as Since the 1970s, the NSE model coefficient [57] has become a conventional assessment indicator for evaluating the efficient gain in the performance of the estimated R when compared to the ground-observed R (e.g., [58][59][60]). It ranges from −∞ to 1. The best performance of the estimated R is attained when the calculated NSE model coefficient tends to one. It is calculated as where R 0 (i) and R m (i) are the observed and (reconstructed or predicted) Rs for month k; R 0 and R m are the averages of R 0 and R m , respectively, and max(R 0 ) and min(R 0 ) are the ground-observed maximum and minimum values of R 0 within the time series, respectively.

Results and Discussion
This section examines the reconstruction and prediction performance of all the reconstructed and the predicted R time series from the passive RS surface variables (i.e., NDVI and LST), RS WBVs (i.e., P, ET, S), and their respective standardized forms (i.e., SPI, SETI, and SI). Note that only figures for the RS-observed surface variables averaged over the whole upstream area are shown due to the better performance of RS WBVs than those over Yunnan Province (as shown in Tables 3 and 4). Compared to that of LST and NDVI, the accuracy of RS WBVs for R reconstruction and prediction appeared to be sensitive to the size of our study region used for data averaging, because the accuracy of RS WBVs was highly dependent on the spatial location when compared against the ground-observed data. For instance, research studies yielded a 27-50% relative error for the remotely-sensed P in our study region [61], while a 10-30% relative error for the global inferred ET from MODIS were presented [62,63].  As the baseline results, both the reconstructed and predicted R time series from NDVI and LST captured well the ground-observed R of the TC-CD station pair ( Figure 5). Compared to the ground-observed R, discrepancies in peaks were apparent for the NDVI-reconstructed and predicted R, regardless of the employed models (i.e., linear regression, ANN, and LSTM). This was probably because vegetation index is not a direct water balance variable. Apparent discrepancies in troughs for the LST-reconstructed and predicted R were probably attributable to fewer (more) sunny days during winter (spring), which caused abrupt decreases (increases) in surface temperature. This highlights the deficiency of using LST for R reconstruction and prediction during winter and spring. The LST-reconstructed and predicted R also presented variable time shifts across different years, indicating its unsuitability for runoff reconstruction and prediction. In any case, AI-based models displayed less than 1% improvement when compared to that reconstructed and predicted via linear regression (Table 3). Rs reconstructed and predicted from RS WBVs over Yunnan Province manifested a comparable performance to those of the passive RS surface variables (Table 3), except R reconstructed from S with the lowest PCC value (i.e., 0.809) among RS WBVs. Nonetheless, R reconstructed and predicted from S over the whole upstream area yielded the highest PCC value (i.e., 0.966). R reconstructed and predicted P and ET over the whole upstream area were also consistently better than those over Yunnan Province. This indicated the presence of spatially variable patterns of RS WBVs between the upstream and the downstream due to monsoon variability, while RS WBVs of the whole upstream area appeared to be a better representation of the basin. Moreover, apparent discrepancies were reduced for the P-and ET-reconstructed and predicted Rs when compared to those reconstructed and predicted from the passive RS surface variables (Figures 6 and 7). Although those reconstructedand predicted from AI-based models generated a better performance than those of linear regression, the Rs based on standardized RS WBVs yielded even better outcomes (Figures 6-8) at capturing the peaks and the troughs, because the partly systematic influences were further mitigated by the data standardization process, as also demonstrated by the PCC, NRMSE, and NSE values, when compared to that of linear regression and AI-based models (Tables 3 and 4). In addition, no notable difference was found by comparing the R reconstructed and predicted from MT-CT stations with that of TC-CD stations when RS WBVs of the whole upstream area were being employed, indicating negligible ocean tidal backwater effect. Rs reconstructed and predicted from RS WBVs over Yunnan Province manifested a comparable performance to those of the passive RS surface variables (Table 3), except R reconstructed from S with the lowest PCC value (i.e., 0.809) among RS WBVs. Nonetheless, R reconstructed and predicted from S over the whole upstream area yielded the highest PCC value (i.e., 0.966). R reconstructed and predicted P and ET over the whole upstream area were also consistently better than those over Yunnan Province. This indicated the presence of spatially variable patterns of RS WBVs between the upstream and the downstream due to monsoon variability, while RS WBVs of the whole upstream area appeared to be a better representation of the basin. Moreover, apparent discrepancies were reduced for the Pand ET-reconstructed and predicted Rs when compared to those reconstructed and predicted from the passive RS surface variables (Figures 6 and 7). Although those reconstructedand predicted from AI-based models generated a better performance than those of linear regression, the Rs based on standardized RS WBVs yielded even better outcomes (Figures 6-8) at capturing the peaks and the troughs, because the partly systematic influences were further mitigated by the data standardization process, as also demonstrated by the PCC, NRMSE, and NSE values, when compared to that of linear regression and AI-based models (Tables 3 and 4). In addition, no notable difference was found by comparing the R reconstructed and predicted from MT-CT stations with that of TC-CD stations when RS WBVs of the whole upstream area were being employed, indicating negligible ocean tidal backwater effect.
Overall, most reconstructed and predicted Rs were shown to have PCC values larger than 0.9. Comparing RS WBVs with their standardizations, the reconstructed Rs based on the standardization yielded a 2-6% further increase in terms of PCC, NRMSE, and NSE values. The predicted R based on the standardized S over the whole upstream area resulted in the highest PCC (i.e., 0.978) and NSE (i.e., 0.947), and the lowest NRMSE (i.e., 0.072). This further demonstrated that the standardization was an effective technique in filtering out other environmental influences. After all, careful selection of the time span for the R reconstruction was essential to the passive RS surface variables. Note that the NRMSE values during prediction were, in general, lower than those during reconstruction, because 24 (48) monthly data points were used during prediction (reconstruction) assessment. Moreover, the maximum range between the maximum and minimum runoff values during 2013-2014 was smaller than that during 2009-2012.           Remaining errors might still exist. These include the presence of systematic errors in our model, unknown environmental signals in the RS WBVs, and variable time shifts across different years between the Mekong Basin upstream and the estuary downstream. Moreover, in-situ discharge time series in the river delta is needed for model reconstruction. These are the limitations of this study.

Conclusions
With the passive remotely-sensed (RS) surface variables (e.g., NDVI and LST) along with AI-based models as baseline results, we examined the temporal prediction performance of monthly runoff (R) at the Mekong Delta reconstructed from the upstream remotely-sensed water balance variables (RS WBVs) and their standardizations, which should present more direct causal relationships and mitigate the environmental influences, such as ENSO events and hydrological extremes.
Through linear regression, we found that the reconstructed and predicted Rs based on the RS WBVs were comparable to those based on the passive RS surface variables. The reconstructed and predicted Rs via AI-based models showed less than 1% improvement over those of linear regression. With the presented standardization approach in this study, the reconstructed and predicted Rs based on the standardized RS WBVs yielded a 2-6% further increase in terms of all assessment indicators when compared to the RS WBVs without standardization. The R predicted from standardized water storage averaged over the whole upstream area yielded the highest PCC (i.e., 0.978) and NSE (i.e., 0.947), and the lowest NRMSE (i.e., 0.072). This further indicated the ability of data standardization to mitigate environment influences.
Further study lies in a sensitivity analysis of the standardization process employing different data time periods for assessing the potential influence on the reconstruction and prediction performances under different intensities of monsoonal climate. This in turn provides a more comprehensive assessment of the reconstruction and prediction of R at a higher sampling rate.