Improved Mekong Basin Runoff Estimate and Its Error Characteristics Using Pure Remotely Sensed Data Products

: Basin runoff is a quantity of river discharge per unit basin area monitored close to an estuary mouth, essential for providing information on the ﬂooding and drought conditions of an entire river basin. Owing to a decreasing number of in situ monitoring stations since the late 1970s, basin runoff estimates using remote sensing have been advocated. Previous runoff estimates of the entire Mekong Basin calculated from the water balance equation were achieved through the hybrid use of remotely sensed and model-predicted data products. Nonetheless, these basin runoff estimates revealed a weak consistency with the in situ ones. To address this issue, we provide a newly improved estimate of the monthly Mekong Basin runoff by using the terrestrial water balance equation, purely based on remotely sensed water balance component data products. The remotely sensed water balance component data products used in this study included the satellite precipitation from the Tropical Rainfall Measuring Mission (TRMM), the satellite evapotranspiration from the Moderate Resolution Imaging Spectroradiometer (MODIS), and the inferred terrestrial water storage from the Gravity Recovery and Climate Experiment (GRACE). A comparison of our new estimate and previously published result against the in situ runoff indicated a marked improvement in terms of the Pearson’s correlation coefﬁcient (PCC), reaching 0.836 (the new estimate) instead of 0.621 (the previously published result). When a three-month moving-average process was applied to each data product, our new estimate further reached a PCC of 0.932, along with the consistent improvement revealed from other evaluation metrics. Conducting an error analysis of the estimated mean monthly runoff for the entire data timespan, we found that the usage of different evapotranspiration data products had a substantial inﬂuence on the estimated runoff. This indicates that the choice of evapotranspiration data product is critical in the remotely sensed runoff estimation.


Introduction
River discharge (RD) measured near an estuary mouth is one of the critical water balance components, expressing itself in terms of basin runoff (R) when it is divided by the basin surface area [1]. Continuous monitoring of RD is useful to provide information on flooding and drought conditions, in order to account for potential agricultural losses [2,3]. Unfortunately, insufficient financial budgets have resulted in a number of poorly documented station RD time series over the past four decades [4]. As a result, RD estimated from remotely sensed data products has great potential for supplementing ground measurements worldwide [5,6], highlighting the importance of this study.
At present, the RD (or R) estimated from remotely sensed data can be grouped into the following five categories: (1). Passively sensed gridded data products, such as the normalized difference vegetation index (NDVI), are directly correlated with the in situ water level (or RD) at a nearby gridded point via a linear regression [7][8][9][10][11]. This method is simple and common; however, these data products have no direct causal relationship with RD. (2). Both actively and passively sensed water balance component data products, such as precipitation (P) from the Tropical Rainfall Measuring Mission (TRMM) [12] (hereinafter denoted TRMM-P), evapotranspiration (ET) from the Moderate Resolution Imaging Spectroradiometer (MODIS) [13] (hereinafter denoted MODIS-ET), and terrestrial water storage (S) from the Gravity Recovery and Climate Experiment (GRACE) [14] (hereinafter denoted GRACE-S), are aggregated over the entire basin or sub-basin individually before being correlated with the water level [15] or R [16]. This method is similar to the first method, but the difference lies in the use of the aggregated remotely sensed water balance data over the entire basin or sub-basin. (3). Passively sensed river width is inputted into hydraulic functional models to estimate the RD (e.g., [17][18][19][20][21]). Despite elegant models, high-resolution remote sensing imagery and an accurate in situ measured roughness coefficient are required [22][23][24]. Hence, this method is still dependent on the in situ measurements, which is not desirable. (4). An actively sensed water level (or so-called stage) measured from satellite altimetry is directly correlated with nearby in situ RD either via rating curves (e.g., [25][26][27]) or an ensemble-learning regression technique [28]. This method is direct; however, the nominal size of the radar altimetric footprint (i.e., between 3.5 km and 5 km (e.g., [29])) is normally larger than the river width. This results in an inaccurately sensed water level because the reflected radar signal is contaminated by the riverside (e.g., [30]). (5). Based on the terrestrial or combined land-atmosphere water balance equation, R can be estimated using a combination of a remotely sensed water balance component (RSWBC) and model-predicted data products [31][32][33][34]. This method is independent of in situ measurements; however, R accuracy is highly dependent on reliable modelpredicted data.
Among the five categories, only the last method can be truly independent of the in situ and model-predicted data, largely owing to the remotely sensed S obtained from GRACE since 2002.
GRACE is a satellite mission for observing global time-variable relative gravity, thereby inferring S anomalies with respect to the geoid (i.e., the equipotential surface of the gravity field) with about a three-arc degree gridded spatial resolution at a monthly scale (e.g., [14,35]). This remotely sensed S, which significantly complements the insufficient ground measurements, has gradually become an indispensable component for basin-scale RD or R estimates (e.g., [31][32][33]).
With regard to the combined land-atmosphere water balance equation (e.g., [36]), Syed et al. (2005) [32] first demonstrated its applicability for determining the monthly R estimates of the Amazon and Mississippi river basins by using GRACE-S and a modeled atmospheric moisture flux divergence from the National Centers for Environmental Protection-National Center for Atmospheric Research (NCEP-NCAR) and the European Centre for Medium-Range Forecasts (ECMWF). Similar studies were then conducted for large basins around the world, including the Mekong Basin [37]. The monthly R estimates for the Mekong Basin agree well with the in situ R, with a Pearson correlation coefficient (PCC) of 0.621, a root mean squares error (RMSE) of 46.144, and a Nash-Sutcliffe efficiency (NSE) of −0.521 ( Figure 2 of [33]). Notably, GRACE-S were unavailable before 2003, and their change (∆S) was assumed to be zero [38].
While the lack of GRACE-S data poses an unavoidable uncertainty, ET is another important source of errors for R estimates. Based on the terrestrial water balance equation, Ferreira et al. (2013) [31] estimated R for the Yangtze River Basin using TRMM-P, GRACE-S, and model-predicted ET from the Global Land Data Assimilation System (GLDAS) [39], while remotely sensed MODIS-ET data were used by Chen et al. (2019) [40]. The results of [40] showed better agreement with the in situ data because MODIS-ET data were generated using the Penman-Monteith equation, which is better adapted to cover types of croplands and grasslands than model-predicted ET from GLDAS in the Yangtze River Basin [13,41]. Nonetheless, both the MODIS-ET and GLDAS model-predicted ET should yield different results in different regions.
As a river basin stretching across different latitudes, the Mekong Basin (MB) is subject to a monsoonal climate and is particularly sensitive to seasonality. This basin is affected by the Indian Summer Monsoon, the Western North Pacific Summer Monsoon, and the East Asian Summer Monsoon [42], which intersect and interact with each other [43], as well as El Niño-Southern Oscillation (ENSO) influences [44]. Regardless of the upstream or downstream in the MB, the spatiotemporal pattern of P varies according to the above monsoon interactions, consequently causing the spatiotemporal redistribution of ET and S.
In addition, dams within the MB were actively constructed between 1990 and 2002 [45], owing to reasons of water security concerning different nations [46]. The substantial impact on RD for different seasons and years was quantitatively assessed using an RD time series before and after the dam constructions [47,48]. However, the annual change of RD was insignificant [45,48]. Note that as reliable RSWBC data were only available after 2002, in particular those of GRACE-S, the assessment of the RD change before and after 2002 could be safely avoided because the RSWBC data already recorded the distorted RD. Except for the study of monthly R estimates in the MB using the model-predicted atmospheric moisture flux divergence conducted by Syed et al. (2009) [33], no further study has been conducted.
Given the aforementioned reasons and potential error sources, this study aimed to improve R estimates of the MB based on pure RSWBC data that are truly independent of in situ and model-predicted data. Through the terrestrial water balance equation, the subtraction from among all the RSWBC data averaged over the entire basin should have also mitigated partial systematic effects due to dam operation and ENSO on the R estimates. The R estimates were subsequently compared to those published in [33] in order to demonstrate a marked improvement made in this study. Meanwhile, the mean monthly runoff performance was assessed in order to further explain potential reasons for discrepancies.
This study is organized into three sections: Section 2 presents the data and methodology applied, together with accuracy evaluation indicators; Section 3 presents and compares our estimated R values with a previous published result, along with discussion of the potential discrepancies for a particular RSWBC; and Section 4 summarizes and concludes our findings.

Data and Methods
All the remotely sensed data products and in situ datasets used in this study are summarized in Table 1. Notably, the remotely sensed TRMM-P (i.e., 0.25 • × 0.25 • ), the MODIS-ET (i.e., 0.5 • × 0.5 • ), and the GRACE-S (i.e., 1 • × 1 • or 3 • × 3 • ) data have different spatial resolutions (Table 1). Their differences are attributable to satellite mission purpose, orbit design, and active (i.e., point location sampling with its own light source) and/or passive (i.e., optical image sampling dependent on sunlight) remote sensing (see, e.g., [49]). Direct usage of these RSWBC data in the terrestrial water balance equation would result in rougher R estimates because the TRMM-P and MODIS-ET data have a higher resolution than those of GRACE-S. These RSWBC data ought to be unified to approximately the lowest spatial resolution before usage. Therefore, the TRMM and MODIS data were resampled at a 1 • × 1 • resolution, while the regular GRACE data were interpolated at a 1 • × 1 • resolution, through associated Legendre functions in spherical harmonics. Even when the RSWBC data were unified, we observed that further smoothing could further reduce the discrepancy between R estimates and in situ data. Table 1. Remotely sensed data products and in situ gauge observations used in this study.

River Discharge Data
The RD data time series were obtained from the Mekong River Commission (MRC). The Kratie station (latitude 12.48141 o , longitude 106.01762 o , elevation 26 m) was selected because it is situated~500 km away from the estuary mouth ( Figure 1), so the complexity due to the ocean tidal effect on the Mekong Delta (MD) [50] and the total adjustment effect from Tonle Sap Lake [51] are negligible. Therefore, the river discharge time series of Kratie station was chosen for assessing the performance of the estimated R values in this study. To be consistent with the monthly RSWBC data, the daily RD was transformed into a monthly R by applying a monthly data average divided by the MB surface area down to Kratie station (i.e., 646,000 km 2 ). The runoff time series of Kratie station display a clear seasonal fluctuation pattern during the period from 2003-2012, with troughs relatively flatter and smoother than peaks ( Figure 2). The yearly peaks and troughs appeared in September and April, respectively. This is because the rainfall season is between May and October and is controlled by the Southwest monsoon from the Indian Ocean, while the drought season is between November and April and is influenced by the northeast monsoon from the mainland [52]. The difference in peak values can be potentially attributed to irregular dam operations deployed to reduce the flooding conditions along the MB during the summer, in addition to ENSO variability [44].

Water Balance Gridded Data Products
The monthly quarter-degree gridded P data from the TRMM 3B43 version 7 [53] (i.e., TRMM-P; available at https://disc.gsfc.nasa.gov/datasets/TRMM_3B43_V7/summary, accessed on 16 January 2021) with a near-global coverage were used. Though parts of these datasets were bias-corrected with the Global Precipitation Climatology Centre (GPCP) gauges [54], their accuracy has not yet been confirmed.
The monthly time-variable gravity data from GRACE enabled a conversion into globally gridded S data (i.e., GRACE-S). The regular release (RL) 05 and 06 (hereinafter abbreviated RL05 and RL06) and mascon solutions generated by the Center for Space Research (CSR) were used. Though the monthly gridded RL06 mascon data (http://www2 .csr.utexas.edu/grace/RL06_mascons.html, accessed on 16 January 2021) at a 1 • × 1 • resolution could have been used directly, the RL05 and RL06 datasets (http://icgem.gfzpotsdam.de/series, accessed on 16 January 2021) were provided in terms of gravity Stokes coefficients (GSCs) up to the degree of 60, which is equivalent to a spatial resolution of about three arc degrees. Post-processing steps were necessary for the RL05 and RL06 data, including adding the degree-1 term and replacing the second zonal coefficient, C 20 , of the GSCs to account for the improved geocenter motion and the geoid, respectively [56,57]. Destriping and Gaussian filtering (with a 350-km radius) were applied to mitigate correlated noise in space [58]. Subsequently, the monthly gridded S time series were obtained using Equation (14) in [59] divided by the water density.
The monthly one-degree gridded S data were generated from the GLDAS soil moisture data [39], which are the sum of 0-10, 10-40, 40-100, and 100-200 cm layers. Four landsurface models, including the Common Land Model, Mosaic model, variable infiltration capacity (VIC) model, and the Noah model [39], were used.
While both remotely sensed precipitation and evapotranspiration data have been compared with in situ data previously [40,41,60], remotely sensed water storage data could only be compared to the model-predicted water storage data (i.e., GLDAS models) due to the lack of in situ water storage measurements. All the above water balance gridded data products were visualized for a specific month in Summer (i.e., August 2008) in the MB (Figure 3), revealing the spatial heterogeneities of all RSWBCs. For instance, both TRMM-P and MODIS-ET presented larger values in the east of the main river channel than in the west. However, the GLDAS CLM-ET yielded a different spatial pattern when compared to MODIS-ET. All GRACE-S data presented similar patterns, except for the downstream MB for the GRACE mascon solution. These spatial heterogeneities potentially yield high uncertainties in R estimates when subtraction among RSWBCs through the terrestrial water balance equation, as is further discussed in Section 3.

Methodology
The remotely sensed basin runoff, R, for the MB bounded area up to Kratie station was determined from the terrestrial water balance equation [61] (Figure 4) as where P, ET and ∆S are the precipitation, evapotranspiration, and the change of terrestrial water storage (in terms of mm/month), respectively. The change of terrestrial water storage in Equation (1) was defined in the following form: where S i is the terrestrial water storage anomaly with respect to the geoid of month i. Note again that the RSWBC data products were presented with different spatial resolutions, as stated above. These gridded data products were unified to a 1 • × 1 • resolution. Gridded values of each data product within the MB were spatially summed up individually within the MB bounded area up to Kratie station (see Figure 1) each month to form a monthly time series for the period from 2003-2012. Due to the spatial heterogeneities for all RSWBCs within the entire MB, rougher R estimates generated from Equation (1) were obtained, and hence, further smoothing of the R estimates was required. Three-month, five-month, and seven-month moving averages were applied and compared against the in situ data. It was found that the three-month moving average was the best, because the five-month and seven-month moving averages caused a time-shift in the R estimates, resulting in larger discrepancies with the in situ data.

Accuracy Evaluation Indicators
The estimated time series of remotely sensed R in the MB were subjected to accuracy evaluation against the in situ R obtained from the Kratie station. Consistent with the accuracy evaluation indicators employed in classical and operational hydrology, the Pearson correlation coefficient, the Nash-Sutcliffe efficiency model coefficient, and the root mean square error and its normalized form (NRMSE) are conventional assessment indicators (as in, e.g., [62,63]). The mean absolute percentage error (MAPE) was also employed. PCC, with its value within ±1, measures the strength of collinearity between two variables. It is formulated as ( Commonly employed in the field of hydrology, the NSE model coefficient, with its value ranging between −∞ and 1, yields the efficient gain in the performance of the estimated R against the in situ R time series [64]. The best performance is achieved when the NSE equals one. It is formulated as The RMSE is another commonly employed evaluation indicator for comparison of the estimated quantity against the in situ one. It is formulated as The NRMSE is a normalized version of the RMSE, yielding a relative RMSE. It is defined as The MAPE represents the average of the absolute error in terms of a fraction [65], calculated as: where R 0 (i) and R m (i) Ymare the in situ and remotely-sensed R values for each month i; R 0 and R m are the mean of R 0 (i) and R m (i); and max(R 0 (i)) and min(R 0 (i)) are the maximum and the minimum of the in situ R, respectively.

Results and Discussion
The R estimates ( Figure 5) based purely on RSWBCs (namely TRMM-P, MODIS-ET, and several versions of GRACE-S) were compared with a combination of remotely sensed and modeled data products (namely GLDAS CLM-ET and several versions of GRACE-S, and MODIS-ET and several versions of GLDAS model-predicted S) (Figures 6 and 7), using the same TRMM-P data ( Table 2). Our analysis indicated that the R estimates based purely on RSWBCs yielded better results than the combination of remotely sensed and modeled data products. No matter the data products, all the estimated R values from these combinations yielded significantly higher PCCs, RMSEs, and NSEs than those of [33] (PCC of 0.621, RMSE of 46.144 mm/month, and NSE of −0.521). The best-estimated R at Kratie station, based purely on RSWBCs (i.e., MODIS-ET and GRACE-S CSR Mascon), resulted in a PCC of 0.836, RMSE (NRMSE) of 27.465 mm/month (0.169), and NSE of 0.665. A similar situation applied to other evaluation metrics in terms of ranking among different estimations of R. By evaluating the best-estimated R in terms of PCC, we were able to show a remarkable improvement when compared to the PCC of 0.621 obtained in [33].    Despite the remarkable improvement, all the R estimates were not as smooth as the in situ runoff, in particular the peaks and troughs. We observed a consistent discrepancy in troughs during dry seasons when compared to the in situ runoff. Using MODIS-ET and different GRACE-S data products, the R estimates showed an underestimation of the runoff during the dry seasons ( Figure 5), while the one based on the CLM-ET yielded an overestimation ( Figure 6). This highlights the unreliability of ET during the dry seasons, irrespective of whether observed or model-predicted ET data was used. As mentioned, GLDAS CLM-ET has a different spatial pattern from MODIS-ET, in addition to their differing magnitudes ( Figure 3). As a result, ET can be one of the main error sources.
Employing MODIS-ET and several types of GLDAS model-predicted S data, the R estimates showed an underestimation (or overestimation) of the runoff during the dry (respectively, wet) seasons when compared to the in situ runoff (Figure 7). This indicates that the GLDAS model-predicted changes of terrestrial water storage (i.e., ∆S) are always overestimated (underestimated) during the dry (wet) seasons. Moreover, variable time lag between the R estimates and the in situ runoff was present during the time period. These predicted changes in terrestrial water storage imply unrealistic terrestrial water storage conditions generated from GLDAS. This is the reason why the R estimates from MODIS-ET and different GRACE-S data products were better when compared against the in situ runoff (Table 2).
Overall, the in situ runoff presented a smaller range of fluctuations than our estimates. This is attributable to the result of the operation of cascaded dams covering the entire MB, which increases (or decreases) the runoff during drought (respectively, flooding) seasons with nearly no annual change (see, e.g., [45,66]). As a result, the in situ runoff presented smoother troughs than peaks. We speculate that the dam operations are not managed well during summer.
Note also that the runoff time series presents alternately high and low peak levels during 2010-2012, when El Niño and La Niña events occurred alternately with moderate intensities (https://ggweather.com/enso/oni.htm, accessed on 16 January 2021). This intensified the abnormal ranges of the precipitation and evapotranspiration, in addition to their spatiotemporal variations, and hence caused a larger runoff during the 2011 summer. This increased the discrepancy in peaks between the R estimates and the in situ R.
To mitigate the external influences due to climate variability, such as ENSO, on the R estimates, a three-month moving average was further applied to each RSWBC before R estimation using Equation (1), as mentioned. It can be seen that all estimated R values in this study were further improved, with the PCC of the best estimated R based purely on RSWBCs increasing from 0.836 to 0.932 (Table 2). This was mainly owing to the reduction of the discrepancy against the in situ R during dry seasons and for the peak during 2011 (Figure 8). To further understand the error characteristics among the estimated R values from different data combinations, Taylor's diagram was employed to visualize the correlation coefficient, the RMS difference (RMSD), and the standard deviation between the estimated and in situ values of R in a quarter pie chart [67]. Notice that the RMSD considers the mean offset between the estimated and in situ runoff, which is different from the RMSE.  (Table 2). In addition, large discrepancies between the estimated R values generated by the GLDAS model-predicted S and the in situ R indicate that the GLDAS model-predicted S data products are not suitable for runoff estimation when compared to GRACE. After applying the three-month average for each data product before estimation, further improvement was achieved for all the estimated R values (Figure 10).  To further examine whether ET was influential on the estimated R values or not, the discrepancy between the estimated R values and the in situ runoff was evaluated with the estimated mean monthly R for MODIS-ET with several versions of GRACE-S and for GLDAS CLM-ET with several versions of GRACE-S, after applying a three-month average for each data product (Figures 11 and 12). Comparing the mean monthly R values between the two dataset combinations, substantial differences can be observed (Figures 11 and 12). Though the resulting R generated purely from the RSWBCs best matched the in situ R, the estimated mean monthly R presented various degrees of discrepancies for different seasons. Figure 11. The mean monthly remotely sensed runoff generated from the MODIS evapotranspiration data and several versions of the GRACE terrestrial water storages data after applying a three-month moving average and the in situ runoff at Kratie station. Figure 12. The mean monthly remotely sensed runoff generated from the GLDAS CLM evapotranspiration data and several versions of the GRACE terrestrial water storage data after applying a three-month moving average and the in situ runoff at Kratie station.
Apparently, the estimated mean monthly R for MODIS-ET with the GRACE CSR mascon was inconsistent with that of CSR RL05 and RL06 for the two dataset combinations (Figures 11 and 12) because of the different data pre-processing strategies of RL05 (or RL06) and the mascon solutions. In the case of MODIS-ET with CSR mascon, the estimated mean monthly R was consistent with the in situ R for most months, except for being underestimated in August, September, and October ( Figure 11). However, in the case of MODIS-ET with GRACE RL05 and RL06, the estimated mean monthly R values were underestimated (or overestimated) between September and March (respectively, May and July). As for the case of GLDAS CLM-ET with the CSR mascon, the estimated R was overestimated for all months except August and September ( Figure 12). In the case of GLDAS CLM-ET with CSR RL05 and RL06, the estimated mean monthly R was overestimated for most months of the year, except for being matched well in August and September. Overall, the estimated mean monthly R generated from MODIS-ET with the CSR mascon exhibited smaller discrepancies than those from GLDAS CLM-ET with the CSR mascon. No matter whether MODIS-ET or GLDAS CLM-ET was used, unavoidable uncertainties arose during the remotely sensed R estimation process. After all, a number of research studies have complained about the uncertainty of ET data products at global or regional scales (e.g., [68,69]).

Summary and Conclusions
Previous study estimated the monthly runoff in the Mekong Basin according to the water balance equation, based on a hybrid use of the remotely sensed water balance component and modeled data products. In this study, we estimated the runoff purely using the RSWBC. After applying a three-month moving average to the RSWBCs, the best runoff estimate yielded a PCC of 0.932 when compared to the in situ runoff. This was a remarkable improvement when compared to the PCC of 0.621 obtained by Syed et al. (2009) [33]. We also found that dam operations and ENSO had a substantial impact on the R estimates when compared against the in situ runoff.
To further examine the error characteristics, we computed and compared the estimated mean monthly runoff. We found that, given the same TRMM precipitation and GRACE terrestrial water storage data, our best runoff estimate obtained purely from the RSWBCs was consistent with the in situ runoff for most months, whereas the runoff estimate based on GLDAS CLM evapotranspiration overestimated the runoff for most months. This finding indicates that the use of different ET data products has a substantial influence on the remotely sensed runoff estimates, which has not been previously investigated.
Further improvement lies in the separate estimation of upstream and downstream runoffs and in examining the natural boundaries (e.g., topography) before the estimated basin runoff, as there are potential differences between upstream and downstream water balance component properties, in addition to a potential time lag between the upstream and downstream regions of a basin.

Data Availability Statement:
The download link for the remotely sensed precipitation, evapotranspiration, terrestrial water storage data, along with the GLDAS-modeled predicted data are provided in Section 2.2 within the text. The in-situ Kratie station time series data is supplied by the Mekong River Commission (http://www.mrcmekong.org, accessed on 16 January 2021).