Improving Stage–Discharge Relation in The Mekong River Estuary by Remotely Sensed Long-Period Ocean Tides

: Ocean tidal backwater reshapes the stage–discharge relation in the ﬂuvial-to-marine transition zone at estuaries, rendering the cautious use of these data for hydrological studies. While a qualitative explanation is traditionally provided by examining a scatter plot of water discharge against water level, a quantitative assessment of long-period ocean tidal e ﬀ ect on the stage–discharge relation has been rarely investigated. This study analyzes the relationship among water level, water discharge, and ocean tidal height via their standardized forms in the Mekong Delta. We found that semiannual and annual components of ocean tides contribute signiﬁcantly to the discrepancy between standardized water level and standardized water discharge time series. This reveals that the long-period ocean tides are the signiﬁcant factors inﬂuencing the stage–discharge relation in the river delta, implying a potential of improving the relation as long as proper long-period ocean tidal components are taken into consideration. By isolating the short-period signals (i.e., less than 15 days) from land surface hydrology and ocean tides, better consistent stage–discharge relations are obtained, in terms of improving the Pearson correlation coe ﬃ cient (PCC) from ~0.4 to ~0.8 and from ~0.6 to ~0.9 for the stations closest to the estuary and at the Mekong Delta entrance, respectively. By incorporating the long-period ocean tidal height time series generated from a remotely sensed global ocean tide model into the stage–discharge relation, further reﬁned stage–discharge relations are obtained with the PCC higher than 0.9 for all employed stations, suggesting the improvement of daily averaged water level and water discharge while ignoring the short-period intratidal variability. The remotely sensed global ocean tide model, OSU12, which contains annual and semiannual ocean tide components, is capable of generating accurate tidal height time series necessary for the partial recovery of the stage–discharge relation. incorporating the ocean tidal components into the stage–discharge relation for a partial relation recovery in the MD. The relation among WD, WL, and ocean tidal height data time series are analyzed via their standardized forms. The ocean tidal components generated from remotely sensed OSU12 global ocean tide model are substituted into the resulting model relation generated from the analysis. The ﬁtted model relation is subsequently applied for estimating WD from ocean tidal height and WL. A quantitative evaluation of the estimated WD against the observed hydrological data is also presented. validation, H.S.F. and H.P.; formal analysis, H.S.F. and H.P.; investigation, H.P.; resources, H.S.F. and L.W.; data curation, H.S.F. and J.G.; writing—original draft preparation, H.P. and H.S.F.; writing—review and editing, H.S.F., H.P., L.W. and J.G.; visualization, H.P.; supervision, H.S.F.; project administration, H.S.F.; funding acquisition, H.S.F. and L.W.


Introduction
Accurate water level (WL) and water discharge (WD) measurements are fundamental to various hydrological applications, including flood forecasting, design and operation of conservancy facilities, as well as water and sediment budget analyses [1,2]. However, due to economy, politics, and topography along a river [3], the spatial distribution of hydrological stations is both sparse and uneven, along with inconsistent and missing datasets [4]. In order to complement the above deficiency of observed datasets, (The topography dataset, called earth_relief_30s, is a derived product of SRTM15+ [38], which is obtainable from http://mirrors.ustc. edu.cn/gmt/data/).

Datasets and Assessment Metrics
In this study, in-situ data from hydrological stations, tidal gauge data, and OSU12 global ocean tide model were analyzed. Table 1 summarizes the essential information about these datasets.

In-Situ Hydrological Data
Station data time series within the MD were obtained from the Mekong River Commission (MRC) (http://www.mrcmekong.org). Acoustic Doppler Current Profiler (ADCP) was applied to gauge flow velocity for deriving precise discharge, according to MRC [39]. To compare between the two main subdivided branches within the MD, Tan Chau and My Thuan stations along the Mekong River, and Chau Doc and Can Tho stations along the Bassac River were used. Situated at the entrance of the MD [27], the Tan Chau and Chau Doc stations are, respectively,~220 and~240 km away from the estuary mouth. Both stations are in the middle between the Tonle Sap Lake and the estuary mouth, where the regulation effect of the lake and the ocean tidal backwater effect are minimized. Being the closest hydrological stations to the estuary mouth, My Thuan and Can Tho stations are subject to the backwater effect caused by landward ocean tidal propagation, which is clearly shown in the data time series [27]. Hence, the comparison between upper and lower station pairs allows us to further quantify the extent of the ocean tidal backwater effect.
Note that WD data of Tan Chau station were missing in 2001, 2002, and 2007. To be consistent with the time span of other WD data, the station time series spanning from January 2003 to December 2006 were selected for investigation, while those from January 2009 to December 2014 were employed for validation. Given the different temporal resolutions among WL, WD, and in-situ ocean tidal data time series and in order to isolate signals unrelated to hydrology, a Butterworth filter was applied to these time series for suppressing periodic fluctuations shorter than 15 days (e.g., diurnal, semidiurnal, etc.). The mean, maximum, and minimum values of those time series are summarized in Table 2.  (Figure 2a-d). Can Tho and My Thuan station time series show a larger ocean tide backwater effect than those of Chau Doc and Tan Chau stations. By comparing WL with WD time series, WD lags behind WL by approximately a month. This fact is more pronounced for stations closer to the estuary mouth (i.e., Can Tho and My Thuan) than their upper counterparts (i.e., Chau Doc and Tan Chau). Obviously, the annual signal is apparent for all station time series, in which the temporal patterns are highly related to not only seasonal variation of watershed runoff, but also the long-period (e.g., semiannual and annual) ocean tidal components, as shown in Figure 2e. As a consequence, external information obtained from the tide gauge or ocean tide model data near estuaries can be potentially used for removing the effect of long-period ocean tidal components, which is the objective of this study.

Sea Level Data from Tide Gauge Station
A tide gauge measures sea level time series at selected locations along the coasts [40]. Vung Tau is the closest tide gauge station to Mekong estuary chosen for relating the long-period ocean tidal variation to WL within the MD. Spanning from 2003 to 2014, the sea level time series at Vung Tau station were recorded on an hourly interval. This dataset is provided by the Hydrological and Environmental station network center in Vietnam and can be obtained from http://www.ioc-sealevelmonitoring.org/ station.php?code=vung. Figure 2e shows filtered and original hourly time series of the tidal gauge data. Fast Fourier transform (FFT) was applied to identify different periodic components of the time series. The highest power spectra are located at both diurnal and semidiurnal ocean tidal components (Figure 3a), which are unrelated to hydrological signals. In order to be consistent with WD and WL time series' daily sampling rate, the hourly tidal height time series are averaged daily after filtering high-frequency components via the Butterworth filter. This process, to a large extent, suppresses or removes the short-period ocean tidal components via the low-pass filtering process, and hence, reducing the effect on long-term ocean tidal components [41][42][43] (Figure 3b). Compared with the unfiltered time series, only semiannual and annual tide components are apparent in the processed time series.

Global Ocean Tide Model Data
A global ocean tide model contains gridded in-phase and quadrature amplitudes (or equivalently amplitude and phase) for major tidal constituents, allowing us to generate ocean tidal height in the absence of tide gauge stations along the coasts [44,45]. Although many remotely sensed ocean tide models (e.g., FES2014, GOT4.8, NAO99.b, TPXO8, EOT11a, DTU10, HAMTIDE, OSU12, etc.) are available for the purpose of our study, the OSU12 model, with a 0.25 • × 0.25 • spatial resolution [46,47], was employed to generate long-period tidal height time series at grid points near Mekong and Bassac river estuaries (Table 3), because it contained long-period tides and was derived purely from remotely sensed satellite altimetry data. Notwithstanding smaller amplitude when compared with semidiurnal and diurnal tides, long-period ocean tidal components are influential to daily and monthly average WL time series. As shown in Figures 2 and 3b, long-period ocean tidal components are likely related to the discrepancies between the pattern of WL and WD time series. It is appropriate to calculate the ocean tidal height time series, TH(t), at time t from the in-phase, H 1 , and quadrature amplitudes, H 2 , of S a , S sa , and M m tides, which can be formulated as: where T i is the period of each long-period ocean tidal component i. Note that both the in-phase and quadrature amplitudes are with respect to Greenwich Meridian with the starting time, 0:00 AM, 1 January 2002 (UTC +0).

Assessment Metrics
To evaluate the estimated WD against in-situ WD time series in Section 4, three assessment metrics, R-Square, the Pearson correlation coefficient (PCC), and the Nash-Sutcliffe efficiency (NSE) coefficient, are employed.
R-Square, ranging between 0 and 1, describes how much the variation of in-situ WD, WD g , is explained by the estimated WD, WD e , generated from the model. The closer the value to 1, the better the model fitting to the WD g . R-Square is equal to the quotient of sum of squares regression (SSR) divided by sum of squares total (SST), and can be defined as: PCC, ranging between −1 and 1, describe how strong the linear relationship between WD e and WD g , which is defined as: where WD e and WD g are the mean of WD e and WD g , respectively. Notably, for the power function relating WL to WD, logarithmic transform is applied to obtain the log-linear relation between the two variables in order to assess their correlation. To highlight the difference, PCC was used to represent the linear relationship between WD e and WD g , while "correlation coefficient" appeared in each figure of this study referred to the log-linear relation between WD and WL, as shown in Equation (6) below. NSE coefficient, ranging from −∞ to 1, describes the gain in the performance of WD e against WD g . The closer the NSE coefficient to 1, the better the performance of the estimation [48]. It is defined as:

Data Analysis and Methodology
This section explores the relations among ocean tidal height, WL, and WD time series over our study region, so as to illustrate the interaction between fluvial and oceanic factors along with their combined effects on WL and WD data. For an ideal hydrological station location where WL and WD are purely influenced by the fluvial process, WL and WD are related by a power function [49] expressed as: where a, b, c are the scaling coefficient, the offset of WL and the exponent of power function, respectively. However, in reality, the stage-discharge phase diagram between WL and WD appears as random data points with trends (i.e., Can Tho and My Thuan stations) and elliptical curves (i.e., Chau Doc and Tan Chau stations) in the MD (Figure 4). The logarithmic transform of Equation (5) allows the conversion into log-linear relation, expressed as: such that Equation (6) measures a linear relationship between ln(WD) and ln(WL − b). All "correlation coefficients" displayed in all stage-discharge phase diagrams were calculated based on ln(WD) and ln(WL − b), as mentioned in Section 2.4. Compared to those of the other two stations, the rating curves between WL and WD of Can Tho and My Thuan stations yield lower correlation coefficients because they are more significantly affected by the ocean tidal backwater.

Data Analysis of Backwater Influence on Water Discharge (WD) and Water Level (WL)
Although the phase diagram between WL and WD in the tide-dominated area appears elliptical, the patterns of the deviation from the rating curves are presumed to be analyzable by different ocean tidal components. Through FFT, the most pronounced periods are 365 days, 182.5 days, and 14.7475 days in both WD and WL time series.
The relative power (to the signal with the largest power) and initial phase of each signal are displayed in Table 4. For an ideal stage-discharge relation (i.e., power function relation), WL and WD are positively correlated. The signals of WD and WL with the same period should have the same initial phase and similar relative power. However, we found that the initial phase of WD and WL time series of the four stations are different from each other. Firstly, annual signals (i.e., 365-day period) of Can Tho and My Thuan present different initial phases between WD and WL, in particular WL, with its initial phases~30 • lower than that of upper counterparts. This indicates that annual tides can cause around a one-month time lag between the lower and upper stations. A similar situation applies to that of the semiannual signal, but to a lesser extent. Secondly, the initial phase of the half-monthly signal (i.e., 14.7475-day period) of WD and that of WL present the phase difference between 160 • and 220 • . This shows that the WD is inversely proportional to WL with an additional time lag. In other words, the WD increase (decrease) when the WL decrease (increase), implying that the half-monthly signal of WL and WD interacts with each other seasonally and alternately. This fact further indicates the half-monthly signal is of two origins: land and ocean, which is supported by physical explanations from Guo et al. (2020) [50] and Jay (1991) [51]. Half-monthly signals of the Can Tho and My Thuan stations yields a much larger relative power than their upper counterparts, indicating the damping effect on the amplitude and changing phase when propagating inland via the estuary mouth. Since these half-monthly signals have different changing ratios for inland propagation direction with annual tide components, a band-pass filter (e.g., Butterworth filter) was applied to remove this signal from tidal-influenced time series for consistency.
To further analyze the interaction between oceanic and fluvial effects, the variation of WD, WL, and TH time series are compared via their standardized forms, x s , expressed as: where x is the average value of xx time series, and NN is the number of data in the time series. The standardized WD, WL and TH (i.e., WD s , WL s , and TH s respectively) are compared for the four stations, respectively, in Figure 5. As shown in Figure 5b,d, it is clear that the standardized WL time series are highly correlated with standardized WD time series, they reach the maximum values in early September and minimum in March and April simultaneously. Influences from ocean tide are minor, and the ocean tidal height series reaches its maximum and minimum values in different months. However, in Figure 5a,c, there exists large deviation between WD and WL time series. In the lower stations, the WL reaches its minimum and maximum value about a month later than WD, consistent with the initial phase difference of around 30 • stated above (Table 4). For most cases, WL (red line) is set between WD (blue line) and TH (yellow line), emphasizing the influence of the ocean on WL. Previous studies attribute this phase difference to floods up and down or a time lag caused by tidal propagation [52]. Since this phenomenon is more apparent in stations closer to the estuaries, we speculate it is mainly caused by the mixing of fluvial-dominated and marine-dominated fluctuations at the annual and semiannual scale. Theoretically, when two signals with the same period (T) are combined, the new signal will have the same period (T) but a different initial phase (φ 3 ), is expressed as: where A 1 , A 2 , and A 3 are three different amplitudes, φ 1 , φ 2 , and φ 3 are three initial phases, and t refers to time epoch. The proof of Equation (8) is listed in Appendix A.
Since annual and semiannual signals are major components in the WL, WD, TH time series, the fluvial-dominated annual (semiannual) signal and marine-dominated annual (semiannual) signal form a mixed annual (semiannual) WL time series in the MD. For both WL and TH fluctuating in the vertical direction, the WL will be potentially corrected by annual and semiannual ocean tidal components if the TH time series are involved in the power function fitting process. This will be further explored in the next subsection. However, Equation (8) does not work for ocean tidal components shorter than half-monthly one, because of possible non-linear interaction among fluvial factors, bottom topography of an estuarine channel and ocean tidal backwater. This leads to the non-linear change of amplitude and phase during the inland propagation process.

Incorporating Long-Period Ocean Tidal Components into Rating Curve
Before incorporating long-period ocean tidal components into the rating curve, short-period fluctuations in both WD and WL time series, including short-period diurnal and semidiurnal ocean tides, have to be removed. As mentioned in Section 2.1, this can be achieved by a Butterworth filter that suppresses all high-frequency signals with a period shorter than 15 days. Consistent with the filtered time series (Figure 2), the rating curve with filtered time series of WD and WL at the four stations are plotted in terms of phase diagrams ( Figure 6). Compared with those in Figure 4, it is clear that the correlation coefficients have been improved significantly ( Figure 6). However, the elliptical loops are still apparent, indicating a time lag between WL and WD time series, as mentioned in Section 3.1.
Given that the relationship between TH s , WD s , WL s has been analyzed in Section 3.1, it is likely that the elliptical looping phenomenon is largely due to semiannual and annual ocean tidal components. For both WL and TH fluctuating in the vertical direction, the TH s time series were applied to separate the tide-induced fluctuation from the WL time series through a fitting process. The WL free from tide influence, WL free , is defined as: where, α is a coefficient that rescales the TH s . Consequently, the relationship among WD, WL, TH s can be represented by: where a, b, c, and α are to be determined from the observed WD, WL, and TH time series. Through a non-linear fitting [53,54], a, b, c, and α can be determined, and WL free is obtainable.

Results and Discussion
The rating curves of WL free and WD are shown (Figure 7), yielding much higher correlation coefficients when compared to the rating curves of original WD and WL time series. Although rating curves of the Can Tho and My Thuan stations still display lower correlation coefficients than their upper counterparts, significant improvement has been observed. Additionally, the elliptical looping phenomenon related to 'time-lag' between WL and WD is also diminished for all four selected stations. As a countercheck, the time series of WD and WL free for the two stations close to the estuary are shown in Figure 8, revealing no apparent time lag.  In the absence of tide gauge data, the TH time series generated from a global ocean tide model would be a viable alternative, because it can provide ocean tidal height components for the global ocean. The method for obtaining model-derived TH time series has been stated in Section 2.3. The model-derived TH series and in-situ gauged tidal height time series are displayed, manifesting high similarity with each other (Figure 9). Employing the above methodology, the rating curves of the four stations have been recovered using model-derived TH time series (Figure 10).  Although the correlation coefficients of the rating curve fitting generated by model-derived TH time series (Figure 10) are slightly lower than by in-situ tide gauge TH time series (Figure 7), the improvement is considerable when compared to the original unmodified rating curves. This is because most global ocean tide models are derived from satellite altimetry, with the model accuracies lower than that of gauge-derived TH, in particular coastal regions [39]. Given the above results, it is appropriate to use an ocean tidal model to partially recover rating curves over tide-dominated regions.
To evaluate the accuracy of the recovered rating curves, the estimated WD time series via Equation (10) generated from both model-derived and in-situ TH time series are compared with the in-situ WD during 2003-2006. Table 5 lists all determined coefficients of Equation (10) along with the assessment metrics (i.e., R-Square, PCC, NSE) that assess the estimated WD against in-situ WD. For both WD estimated based on in-situ and model-derived TH data, all the assessment metrics yield high-correlation values at all the four stations, suggesting our method can partially recover the tide-free WL for estimating WD. Overall, the recovered stage-discharge relation is capable of predicting a relatively reliable WD. These results also validate the data analysis in Section 3. To highlight the importance of tidal separation by the term −α × TH s in Equation (10), the PCC values with different combinations of coefficient b, c, and α were calculated with their best fitted a fixed. Taking Can Tho station as an example (Figure 11), b and c impact the PCC values (i.e., > 0.9) significantly, only if α is~0.16. The same holds for My Thuan station. Therefore, adding the term −α × TH s to the conventional power function (i.e., Equation (5)) is necessary for the improvement of the stage-discharge relation in the MD, which is in the fluvial-to-marine transition zone. In summary, we found that an appropriate α is a prerequisite for the PCC larger than 0.9. To assess the applicability of the determined coefficients of Equation (10) for Can Tho and My Thuan station time series during 2003-2006, these coefficients were directly employed for the analysis of the WL and TH data time series during 2009-2014. The predicted WD were then compared with the monthly in-situ WD (Figure 12), since only monthly WD are available for Can Tho and My Thuan stations. Hence, WL and WL f were monthly averaged before the comparison. For both Can Tho and My Thuan stations, tide-free WL, WL f , leads to higher correlation coefficients and diminishes the looping curve problem to a large degree. This indicates that coefficient α appears to be stable during our study period. Despite a substantial improvement made in this study, small deviations from the ideal power function still exist, particularly for the two stations closest to the estuary mouth. After all, the interaction between fluvial and marine processes are complicated near estuary mouths [55]. Remaining effects cannot be neglected. For instance, WD should pose a non-negligible effect on the tidal propagation along the river channel during the wet season. Overland flows inward or outward from the Tonle Sap Lake would likely be another important factor affecting the stage-discharge relations, because this lake operates as a natural reservoir that regulates Mekong river discharge from the river delta to the coastal ocean [34,56,57]. Erosion and deposition alter hydraulic geometry and increase channel bottom friction and, hence, contribute to the potential instability of the stage-discharge relation. Furthermore, numerous clusters of dams were built along the main stream of Mekong river, which may also alter the stage-discharge relation [21,22]. Sea level rise, which closely connected to salt intrusion and coastal erosion problems may alter the estuarine topography condition, resulting in a secular shift of ocean tidal components [23,26]. Agricultural practices and deforestation also provide additional impact on the evapotranspiration balance of the catchment area. Furthermore, since short-term signals were filtered out or failed to be captured by daily sampling, the short-term variations in WD and WL have not been quantitatively investigated. These considerations represent the current limitations of this study.

Conclusions
Instead of seeking a qualitative explanation of the stage-discharge relation influenced by the ocean tidal backwater effect, this study quantitatively analyzes the relations among water discharge (WD), water level (WL), and ocean tidal components via their standardized forms. We found that annual and semiannual ocean tidal components are significant contributors to the deviation between WL and WD time series. In particular, the annual and semiannual periods of ocean tidal backwater result in the elliptic loop associated with the presence of time lag between WL and WD.
Based on these findings, we adapt the stage-discharge relation to accommodate the effects of annual and semiannual ocean tidal components. It was found that the WD estimated from the de-tided WL yields PCC and NSE values of~0.9. Although the de-tided WL time series generated based on the TH time series from the OSU12 global ocean tide model are slightly less accurate than that of tide gauge data, the ocean tide model is a viable alternative to partially recover the stage-discharge relation for estuaries in the absence of tide gauge stations.
Further improvement lies in identifying remaining effects contributing to the potential instability of the stage-discharge relation, which include the non-negligible effect of seasonal WD on ocean tidal propagation, the Tonle Sap lake regulation effect on the Mekong river discharge, erosion and deposition effects on the hydraulic geometry, and channel bottom friction. The impact of human activities and artificial structure in the Mekong River area, as well as its interaction with climate change, should also be highlighted. Those factors may introduce a long-term change trend into the WL-WD relationship.
The recent remotely-sensed water balance variables with improved temporal resolutions, such as 8-day MODIS evapotranspiration [58], daily TRMM precipitation [59], and daily GRACE terrestrial water storage data products [60], should enable us to compute tide-free WD, which is independent of in-situ measurements based on the water balance equation [36]. This can serve as a countercheck against the in-situ WD for assessing the first two remaining effects.  Acknowledgments: We acknowledge the requested water level and discharge data from the Mekong River Commission (MRC) under an agreement, which is purchased using NSFC Grant No.: 41374010.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The mathematic proof of Equation (8) is shown below: For the convenience of expression, set 2πt T = x.