Upstream GPS Vertical Displacement and Its Standardization for Mekong River Basin Surface Runoff Reconstruction and Estimation

: Surface runoff ( R ), which is another expression for river water discharge of a river basin, is a critical measurement for regional water cycles. Over the past two decades, river water discharge has been widely investigated, which is based on remotely sensed hydraulic and hydrological variables as well as indices. This study aims to demonstrate the potential of upstream global positioning system (GPS) vertical displacement (VD) and its standardization to statistically derive R time series, which has not been reported in recent literature. The correlation between the in situ R at estuaries and averaged GPS-VD and its standardization in the river basin upstream on a monthly temporal scale of the Mekong River Basin (MRB) is examined. It was found that the reconstructed R time series from the latter agrees with and yields a similar performance to that from the terrestrial water storage based on gravimetric satellite (i.e., Gravity Recovery and Climate Experiment (GRACE)) and traditional remote sensing data. The reconstructed R time series from the standardized GPS-VD was found to have a 2%–7% accuracy increase against those without standardization. On the other hand, it is comparable to data that are obtained by the Palmer drought severity index (PDSI). Similar accuracies are exhibited by the estimated R when externally validated through another station location with in situ time series. The comparison of the estimated R at the entrance of river delta against that at the estuaries indicates a 1%–3% relative error induced by the residual ocean tidal effect at the estuary. The reconstructed R from the standardized GPS-VD yields the lowest total relative error of less than 9% when accounting for the main upstream area of the MRB. The remaining errors may be the result of the combined effect of the proposed methodology, remaining environmental signals in the data time series, and potential time lag (less than a month) between the upstream MRB and estuary.


Introduction
River water discharge (RWD) is among the critical hydrological components of river basins measured near the mouth of estuaries [1,2]. Another representation of the RWD at the estuaries is surface runoff (R), in which it is the RWD divided by the total basin area. In order to prepare for unpredictable losses in agricultural products and economy (e.g., [3][4][5][6]), the continuous monitoring of RWDs is necessary for tracking abrupt hydrological changes (i.e., droughts and floods) in deltaic environments. There is no global coverage of gauging network, however, that monitors the RWDs [7]. Apart from this, the frequency of discharge data acquisition has continuously declined since the late 1970s [8] because of insufficient funds for facility maintenance and upgrade [9]. As a result, indirect methods for the RWD monitoring, such as remote sensing (RS), have recently gained increasing interest.
Traditional RS, such as Landsat Thematic Mapper (TM) and its Enhanced TM Plus (ETM+) images, as well as moderate resolution imaging spectrometer (MODIS), have passively recorded instantaneous surface parameters since the 1990s (e.g., [7]). The surface parameters obtained from RS [10][11][12][13], such as flood area inundation, land surface temperature (LST), normalized difference vegetation index (NDVI), and RS-derived geometric variables (e.g., river width), allow the direct correlation with water level or RWD. Except for the RS-derived geometric variables, the foregoing localized RS data yield indirect relationships to the RWD. Although RS-derived geometric variables can also be employed as inputs to infer the RWD through Manning's equation and its modified form (e.g., [14][15][16][17][18][19][20][21]), the accuracy of the estimated RWD is region-dependent because of the technique's ability to detect changes in rivers with short widths [22] and the regional availability of roughness coefficients [23,24].
More recently, space geodetic techniques, such as satellite radar altimetry and the Gravity Recovery and Climate Experiment (GRACE), have extensively been utilized to correlate with the RWD (e.g., [25,26]). These observations, hereafter, are referred to as space geodetic-observed variables. Satellite radar altimetry is able to directly monitor water level variations over water bodies, such as river and lakes (e.g., [27,28]). The water level relates to the RWD via a power function (e.g., [29,30]); hence, the time series of RS RWD is derived by directly correlating the measured satellite altimetric water level with in situ RWD measurements (e.g., [31,32]). The use of a basin-wide RWD estimation that employs multisatellite altimetric water level data has also been demonstrated in the Mekong River Basin (MRB) [25]. The reflected signal contains radar altimetry footprints are contaminated by land surfaces, however, when the river width is smaller than 5 km, such as that found in [33], thereby significantly lowering their accuracy near the riverbank.
Although satellite radar altimetry can actively record along-track surface oscillations of inland water bodies (e.g., [34]), GRACE can directly measure time-variable gravity, thereby making it possible to calculate terrestrial water storage fluctuations at a global or regional scale (e.g., [35,36]). The terrestrial water storage, being one of the water balance components, can physically relate to the RWD (e.g., [26,37]). Therefore, the RWD can be estimated from the terrestrial water storage. GRACEinferred terrestrial water storage (hereinafter called GRACE-S) also allows the calculation of solid Earth vertical surface deformation (e.g., [38][39][40][41]) that is consistent with the recorded global positioning system (GPS) vertical displacement (VD) (hereinafter called GPS-VD) (e.g., [42][43][44]), as demonstrated and validated in different geographic regions [45][46][47]. Conversely, the GPS-VD can also be utilized to infer the terrestrial water storage (e.g., [48][49][50][51]). Given the foregoing GRACEinferred physical quantity that is comparable to the GPS-VD, it is anticipated that the latter can be a potential alternative for reconstructing the RWD; no paper in this regard has ever been published in recent literature.
Essentially, the GRACE-S and its standardization have been recently demonstrated to have a good correlation with water level [52] and R [53]. The reason for the standardization is that it enhances the regional characteristics of the averaged time series [54] when local means and variances are largely different from the regional one [55], thus mitigating systematic influences due to geographic environment [56]. Given the aforementioned similarity between GRACE-inferred quantities and GPS-VD observations, the GPS-VD and its standardization can presumably achieve a similar quality to that of GRACE in capturing the time series of R and standardized R via correlative analysis, respectively.
The Mekong River Delta (MRD), a geographic region that is vital for food (e.g., [57]) and water security (e.g., [58]) in Southeast Asia, is the MRB downstream region immediately before the freshwater is discharged into the coastal ocean. The R of the MRD is affected by dam operation, which increases (decreases) the discharge flow during drought (flooding) seasons under the principle of no significant annual changes [59,60]. Regardless of whether the dam is in the upstream or downstream MRB, the effects of all dam operations accumulate. The accumulated effects that generate biases should be partly systematic for any specific month, year-on-year [61]. These biases can be partly reduced by the subtraction (or differencing) process in the above standardization procedure. The aforementioned reasons justify the potential use of the GPS-VD standardization that is obtained at the upstream of the MRB to correlate with the R in the MRD statistically.
The potential use of the GPS-VD and its standardization is explored in the upstream of the MRB for statistically correlating with the R at a hydrological station in the MRD on a monthly scale via linear regression. The fitted parameters are then utilized to estimate the time series of R of another location having an in situ time series for external validation. The available RS instantaneous data (NDVI [62] and LST [63]), Palmer drought severity index (PDSI) [64], and GRACE-S and its standardization [65], are employed for the purpose of comparison.
The structure of this paper is arranged as below: Section 2 describes the MRB and MRD geography; Section 3 illustrates the datasets and their processing; Section 4 presents the methodology and evaluation metrics; Section 5 demonstrates the reconstruction and estimation of R based on the GPS-VD and its standardized form, while compared to those based on the NDVI, LST, PDSI, and GRACE-S and its standardization; Section 6 summarizes the conclusions.

Geographic Environment of the Mekong River Basin and Mekong River Delta
The Mekong River, originating from the three-river headwater region in the eastern Qinghai-Tibetan Plateau, is the transboundary river across the Southeast Asian continent [66]. Water first flows through the Lancang River in Yunnan Province within the boundary of China, followed by Laos, Myanmar, Thailand, Cambodia, and Vietnam. Spanning 25° of latitude, the total surface area of the entire MRB is around 795,000 km 2 [67] (Figure1).
Situated at latitude 21°-29°N and longitude 94°-102°E, Lancang River within Yunnan Province, China (hereinafter abbreviated as LRWY), constitutes the main component of the upstream MRB with an altitude ranging from 1500 to 4000 m and descending from north to south [68,69]. It is climatedriven and affected by the Indian monsoon [70]; its distinct rainy and dry seasons are the principal seasonal characteristics. During the wet season (May-October), the SW monsoons from the Bay of Bengal yield the almost entire annual rainfall. During the dry season (November-April), however, severe drought may occur [71]. Bare karst geology, which is characterized by rocks with low permeability, surrounds rocks with high permeability at the eastern part of Yunnan Province. This further causes rapid water infiltration into the ground, aggravating drought conditions [72].
The hydrological extremes in the upstream area significantly affect agriculture, living environments, and economy in the midstream and downstream areas where half a billion people live within China and the country's transboundary [73]. In addition, the constructed massive dams have raised water conflicts among various countries in Southeast Asia [74], despite the aim to regulate flow during extreme hydrological periods [59,60]. The dam operated at different times of the year modifies the upstream R, which would adversely impact the water availability downstream (e.g., [75,76]). This indicates that understanding the hydrologically related variables upstream is necessary in order to raise an early alert of extreme events that may occur in the downstream MRB, in particular, the MRD.
The MRD is a transition zone that is seasonally affected by both water discharge of the MRB, and ocean tidal processes propagated landward [57] via the Bassac River and Mekong River branches [77,78]. In addition, the regulation effect of Cambodia's Tonle Sap Lake on the total discharge of the MRB is substantial (e.g., [79][80][81]) before discharging into the northern part of Sunda Shelf.

In Situ Discharge and Passive Remote Sensing Data
Given the regulated effect of Tonle Sap Lake and the backwater effect of ocean tides that govern the total discharge of the MRB, the selection of in situ stations near the estuary mouth is a critical task. Both aforementioned effects are to be minimized in the selected stations. The station at Phnom Penh, despite far away from the mouth of estuaries, is not selected because it intersects with several tributaries that would significantly modify the overall temporal discharge pattern. In this study, the Tan Chau and Chau Doc stations, located at the entrance of the MRD and interior limit of the transition zone [78], are chosen where both the abovementioned effects are minimized ( Figure 1). Despite being closer to the mouth of estuaries, the Can Tho and My Thuan stations are also employed to assess the backwater impact caused by ocean tides on the R time series estimation.
The station discharge time series obtained from the Mekong River Commission (MRC) is available at. The selected WD time series data spanning from January 2012 to December 2014 are extracted because of the shorter time span of GPS data. Since half-daily and daily period of ocean tides are the most dominant ocean tidal forcing, the time series of Tan Chau and Chau Doc pair of stations (hereinafter called TC-pair station) and that of My Thuan and Can Tho pair of stations (hereinafter called MC-pair station) are summed up, respectively, serving as a further monthly averaging process to further mitigate the ocean tidal effect.
The observed RWD (in m 3 /s), accounting for converting second into day, meter into millimeter and dividing by the total basin area (i.e., 795,000 km 2 ), can then be converted into a daily R (in mm/day). The R at a monthly scale (in mm/month) is then computed by adding daily R together. Both the TC-pair and MC-pair station time series exhibit similar temporal patterns ( Figure 2). The LST from MOD11C3 and NDVI from MOD13C2 MODIS products are the traditional RS data that can be downloaded at the Land Processes Distribution Active Archive Center. Both datasets are directly employed to compare against the GPS-reconstructed R time series.

Palmer Drought Severity Index
The PDSI [64], downloaded at [82], is a widely employed index to quantify meteorological drought using worldwide precipitation and temperature data time series to model relative dryness. The index ranges from −10 (dry) to +10 (wet) with a 2.5°×2.5° spatial resolution.
While the SC of the JPL RL 05 are expanded up to degree 90, the SC of the CSR RL05, CSR RL06, and GFZ RL05 are expanded up to degree 60. Using equations in [83], the SC of the CSR RL05, CSR RL06, GFZ RL05 and JPL RL05 can be converted into GRACE-S that is interpolated into a 1°×1°grid.
Except for the CSR-mascon, before converting SC into a GRACE-S time series, the degree-one and terms in SC are added and replaced, respectively, to correct the geocenter motion and the geoid [84,85]. In addition, a de-striping procedure is applied. As tested in this study, a Gaussian filter with a 350-km radius is the optimal radius chosen to reduce the uncertainties arising from correlated errors of TWS data in space at finer resolutions [86,87].
The processed monthly GRACE-S are then used to reconstruct R (i) by directly correlating with the in situ R and (ii) by calculating the recently proposed drought index based on GRACE (hereinafter abbreviated as GRACE-SI) [65], in which it is subsequently correlated with the standardized in situ R.
To reduce GRACE-S anomalies, the median of GRACE-S for each month is employed to calculate GRACE-SI [52] as follows: where , and , are the GRACE-SI and GRACE-S for month j of each year i, respectively; and are GRACE-S and its corresponding sampled standard deviation, respectively. The same calculation procedure for standardizing in situ R is performed using Equation (1).

GPS Vertical Displacement and its Standardization
To obtain a daily GPS-VD in the upstream of the MRB, GAMIT FORTRAN software [88] is employed to preprocess the raw observations of 33 GPS stations that monitor Yunnan Province. The observations during 2012-2014 are available from the Crustal Movement Observation Network of China (CMONOC). This is performed by a stochastically constrained network solution (i.e., assigned with a 5-cm uncertainty) [89] aligned to 24 IGS stations in the ITRF2008 coordinate reference frame that surrounds China [90]. The Earth orientation parameters are also constrained to a priori values listed in the International Earth Rotation Service (IERS) Bulletin B.
Standard correction procedures are applied (i) to constrain the orbits to the IGS final ephemeris; (ii) to make corrections for the first three-order terms of the ionospheric delay in GAMIT [91]; (iii) to make corrections for the tropospheric delay using Vienna mapping function 1 [92] and the global modelof pressure and temperature applied in Geodesy [93]; (iv) to apply the antenna offsets given by the IGS filed antenna correction; (v) to remove the non-tidal atmospheric loading using the MIT correction data files; (vi) to remove the ocean tidal loading by choosing the FES2004 model [94] option in GAMIT [88]; (vii) to remove tides due to the solid Earth and pole according to the IERS standard [95]. Moreover, to yield a purely hydrological signal similar to that of GRACE, vertical crustal displacement caused by ocean bottom pressure changes in the GPS-VD should be corrected using half the daily-modeled data downloaded at the Global Geophysical Fluid Center. Outliers that are larger than twice the standard deviation are discarded. To reduce the aliasing effect and draconitic errors (e.g., ~351 d [96,97]) in the seasonal signal, spectral filtering is applied after the fast Fourier transformation (FFT) of the GPS-VD time series. After the spectral filtering, the first peak is apparently closer to 1 cpy (Figure 3). The 33 filtered GPS spectra are then transformed back into the time domain via the inverse FFT. In general, the height extremes of the GPS time series are reduced after filtering.
Explored at a monthly scale, the daily GPS-VD are averaged every month to yield the monthly VD. To obtain the monthly standardized GPS-VD data, the same calculation procedure is performed using Equation (1).

Correlative Analysis and Runoff Standardization
In this study, all RS and space geodetic-observed variable data are averaged separately over the bounded square covering Yunnan Province and LRWY (Figure 1) before reconstruction and estimation of R in the MRD. Note that smoothing is applied to all the data before the correlative analysis.
The R reconstruction is achieved by a direct linear model fitting with a constant offset, b, and a slope, a, which is given by: where t y and t x are the in situ R (or standardized R) at the selected stations (MC-pair and TCpair) and remotely sensed variables (or its standardized form or indices) at monthly epoch t, except a one-month time shift was applied to NDVI and LST to improve the reconstruction performance. In addition, the GPS-VD, being negatively correlated with GRACE-S and in situ R, was multiplied by a negative one. Note that a and b are the parameters to be estimated during the model fitting process. The reconstructed R when quantitatively evaluated against the same in situ R used for the reconstruction are referred to as internal evaluation, whereas the above estimated parameters determined from the R reconstruction that are subsequently used for the R estimation and evaluated against at other stations in the MRD are referred to as external evaluation; both evaluations assess the methodology employed in this study.
The overlapping time spans of the NDVI, LST, GRACE-S, and GPS-VD during 2012-2014 are employed to correlate the observed R, because the time series of all the remotely sensed variables share similar seasonal patterns to the observed R. However, notable differences are observed. For instance, both the GRACE-S and GPS-VD averaged at the upstream of the MRB against the R time series present a variable time lag over the entire time span (Figure 4). This can be attributable to a slightly different hydrological condition between the upstream and the downstream of the MRB each year under climatic variability (e.g., [98]). Regardless of averaging GPS-VD over the bounded square of the entire Yunnan Province or LRWY, GPS-VD yields an abnormal pattern between January and April of 2014. This can be due to GPS-VD being more sensitive to local hydrological variations or events when compared to GRACE-S [51].
To investigate the improvement resulting from the standardization, the observed R is standardized and compared to the standardized GRACE-S (hereinafter called GRACE-SI) and GPS-VD (hereinafter called GPS-DSI) using Equation (1). The above model fitting procedure is applied to correlate the standardized R with other standardized variables (i.e., PDSI, GRACE-SI, and GPS-DSI).
The estimated parameters and their corresponding standard deviation of Equation (2) is displayed in Table 1. Large offsets (i.e., b) and standard deviation are shown for the in situ R fitted with NDVI and LST. Hence, it is expected that the reconstructed R from NDVI and LST are not fitted well, as shown in Figure 5. The in situ R fitted with GRACE-S and GPS-VD appears to be better than that of NDVI and LST. This linear fitting results preliminarily indicate that our proposed space geodetic-observed variables (i.e., GRACE-S and GPS-VD) are better in R reconstruction.

4.2.Assessment Metrics
Three performance evaluation metrics, including the Pearson correlation coefficient (PCC), normalized root mean square error (NRMSE), and Nash-Sutcliffe efficiency (NSE) model coefficient, are employed to evaluate the reconstructed R against the observed R at in situ stations.
The PCC, which describes the degree of collinearity between two variables, is defined as follows: The NRMSE, being the RMSE divided by the maximum fluctuating range of observations, is defined as follows: The NSE model coefficient, ranging from −∞ to 1, is a performance indicator for evaluating the predictive power of the estimated R compared to the in situ R [99]. When the NSE model coefficient is closer to 1, the performance of the estimated R is better. It is defined as follows: where

Evaluation and Discussion
This section illustrates the accuracy performance of the R time series' reconstruction and estimation for internal and external evaluations, respectively. These time series are from the traditional RSD (NDVI and LST), space geodetic-observed variables (GRACE-S and GPS-VD), and their corresponding standardizations (GRACE-SI and GPS-DSI), including the drought index (PDSI). The internal and external evaluations of results are both applied to the MC-pair and TC-pair stations, so that the estimated R from these two pairs can be compared against each other to assess the residual ocean tidal effect at the estuary. Note that the MC-pair station is closer to the estuary mouth than the TC-pair station with the former being located ~100 km and the latter ~220 km from the estuary mouth. The combined internal and external evaluations could quantify a portion of the systematic ocean tidal backwater effect on both reconstructed and estimated Rs.
Serving as baseline results, both NDVI-and LST-reconstructed Rs from the LRWY exhibit temporal patterns similar to the observed Rs from both the MC-pair (Figure 5a,b) and TC-pair ( Figure  5c,d) stations. However, relatively large differences are presented in peaks and troughs forthe NDVIand LST-reconstructed Rs when compared to against the in situ R (Figure 5c,d). While no apparent time lag for the NDVI-reconstructed R against the in situ R, the in situ R is lagged behind the LSTreconstructed R from March 2013 to September 2014. We speculate that substantial differences between meteorological conditions (e.g., LST) at the upstream and hydrological conditions (e.g., R) at the downstream of the MRB (i.e. MRD) might exist during the above period. Note also that LST is the localized RS quantity that has no direct relationship with the R, as mentioned in the introduction of this study. This is the reason that it is served as a baseline result.
The R reconstructions from GRACE-S and GPS-VD exhibit better results than those of the NDVI and LST-reconstructed Rs ( Figure 6) (Tables 2), whereas their respective standardizations demonstrate even better performances (Figure 7) in capturing the peaks and troughs because a portion of the systematic effects is reduced by the standardization process. A similar situation is also observed in the TC-pair station time series with further reductions in the differences in peaks and troughs (Figure 7d,e,f), as consistently shown by the increase in PCC and NSE values of the reconstructed R in the downstream MRB (Tables 3). The PDSI-reconstructed R exhibits results that are comparable to those of GRACE-SI and GPS-DSI because the PDSI is a hydrometeorological index generated from temperature and precipitation that captures the relative dryness of river basins in relation to river discharge variations [100]. Among all GRACE-reconstructed Rs, the one reconstructed from JPL RL05 yields the best result. The one reconstructed from RL05 and RL06 of CSR show similar performances, whereas that of CSR-mascon solution yields a comparable performance.
By evaluating the differences of the MC-pair and TC-pair estimations in the assessment metrics between Tables 2 and 3 or between Tables 4 and 5, it is found that the usage of TC-pair account for a 1%-3% decrease in the relative error, when compared to that using MC-pair close to the estuary mouth. This is the remaining backwater effect due to ocean tides of the MC-pair in the estuary. The standardization process yields a 2%-7% increase in accuracy, no matter which pair station is used. The R reconstructed from the GPS-DSI yields the lowest NRMSE value when accounting for the LRWY of the upstream MRB only, indicating that it remains subject to the total relative error of less than 9%. It is speculated that the remaining errors may be caused by our methodology, remaining environmental signals in the data time series, and potential time lag of less than a month between the upstream MRB and the MRD.
Overall, the estimated Rs are slightly less accurate than the reconstructed Rs, whereas their relative rankings remain practically the same. This could be partially caused by internal errors that are introduced by the reverse process. The proposed methodology that employs the upstream standardized GPS-VD (i.e., GPS-DSI) is proven to be a viable alternative to the estimation of R in the MRD located at the estuary mouth of the MRB. However, the limitations of this study are that one in situ discharge time series in the river delta or estuary is required, and the GPS stations should be situated on the bedrock surface for observing the elastic deformation due to seasonal water storage changes.

Conclusions
In lieu of employing traditional remote sensing (RS) data for surface runoff (R) reconstruction, the potential use of upstream GPS vertical displacement (GPS-VD) and its standardization (GPS-DSI) for R reconstruction at estuaries on a monthly temporal scale is explored. It is found that the reconstructed R time series at the Mekong River Delta (MRD) from the Mekong River Basin (MRB) upstream GPS-VD are comparable to those from the GRACE terrestrial water storage (GRACE-S) and traditional RS data. All reconstructed R time series from the standardized variables, including GRACE-SI, GPS-DSI, and PDSI, are found to have a 2%-7% increase in accuracy in terms of the NRMSE when compared to those without standardization. The reconstructed R time series from the spatially averaged of the upstream GPS-DSI is shown to be comparable to that reconstructed by the PDSI, but better than those obtained by traditional RS data and GRACE-S. The internal evaluation also demonstrates that the reconstructed R based on GPS-DSI attains a PCC of 0.97 and NSE of 0.94 for both MC-pair and TCpair stations. Despite being slightly less accurate than the reconstructed R, the estimated R exhibits an accuracy that is similar to the above as externally validated via another station location.
The comparison of R reconstruction and estimation from the MC-pair and TC-pair stations indicates that the remaining backwater effect induced by ocean tides yields a1%-3% in the relative error on the estimated Rs in this study. The R reconstructed from the GPS-DSI yields the lowest NRMSE value of less than 9% when accounting for the main upstream area of the MRB (i.e., Lancang River within Yunnan Province). This reveals that the best reconstructed R from the GPS-DSI remains subject to the total relative error of ~9%. This may be caused by our methodology, the remaining environmental signals in the data time series, and the potential time lag (less than a month) between the upstream MRB and the MRD.
Overall, the proposed methodology, which employs the upstream GPS-VD and its standardization, is proven to be a potential alternative for reconstructing and estimating R in the MRB. It is anticipated that the proposed GPS-VD and its standardization can also be applied to basinwide discharge estimations and potentially replace the function of GRACE-S in terms of water balance. Higher temporal resolutions of the reconstructed and estimated Rs can also be achieved via the GPS because the use of daily GPS-VD solutions have become a standard practice, and also offers acceptable accuracy.
Author Contributions: H.S.F. initiated the experimental design, performed data collection, conducted data preprocessing and result interpretation, and wrote the manuscript. L.Z. conducted an entire experiment, including post-processing of remotely sensed data. Y.L. performed GPS data pre-processing. Z.M. visualized data. Y.C. proofread the revised manuscript and provided critical comments. All authors have read and agreed to the published version of the manuscript.