GRACE-Derived Time Lag of Mekong Estuarine Freshwater Transport in the Western South China Sea Validated by Isotopic Tracer Age

: Estuarine freshwater transport has a substantial impact on the near-shore ecosystem and coastal ocean environment away from the estuary. This paper introduces two independent methods to track the Mekong freshwater-induced mass transport by calculating the time lag (or equivalently, the phase) between in situ Mekong basin runoff and the equivalent water height (EWH) time series over the western South China Sea from a gravity recovery and climate experiment (GRACE). The ﬁrst method is the harmonic analysis that determines the phase difference between annual components of the two time series (called the P-method), and the other is the cross-correlation analysis that directly obtains the time lag by shifting the lagged time series forward to attain the highest cross-correlation between the two time series (called the C-method). Using a three-year rolling window, the time lag variations in three versions of GRACE between 2005 and 2012 are computed for demonstrating the consistency of the results. We found that the time lag derived from the P-method is, on average, slightly larger and more variable than that from the C-method. A comparison of our gridded time lag against the age determined via radium isotopes in September, 2007 by Chen et al. (2010) revealed that our gridded time lag results were in good agreement with most isotope-derived ages, with the largest difference less than 6 days. Among the three versions of the GRACE time series, CSR Release 05 performed the best. The lowest standard deviation of time lag was ~1.6 days, calculated by the C-method, whereas the mean difference for all the time lags from the isotope-derived ages is ~1 day by P-method. This study demonstrates the potential of monitoring Mekong estuarine freshwater transport over the western South China Sea by GRACE.


Introduction
The Mekong river basin is the largest tropical water resource in Southeast Asia, dominated by the annual flood pulse. Its estuary annually discharges approximately 475 km 3 freshwater and 160 million tons of sediment into the South China Sea (SCS) [1]. Most heavy sediment settles down quickly along the coast, while buoyant freshwater discharged from the estuary forms a river plume that can be transported on the ocean surface hundreds of kilometers away [2]. The near-field river plume located within 300 km off the coast can be modeled under static-state balance via runoff flux, local current, wind, and ocean tides [3][4][5]. Direct field sampling, such as isotope sampling, enables us to derive the freshwater ages from the estuaries of the Amazon river [6], the Mississippi river [7], the Yangtze river [8] and the Mekong river [9] to their adjacent ocean. Partly due to the costs of field trips and reliability issues [10], it remains a challenge to monitor the continuous spatiotemporal variation in Mekong freshwater transport. Dashed boundary (on the right) is our regional study area with selected (red dots) and unselected (grey dots) isotope samplings, along with their sampling codes (above the dots) and age (in month) (below the dots) (adapted with permission from ref. [9], Copyright 2010 American Geophysical Union).

Oceanic EWH time Series from GRACE
In this study, three versions of GRACE monthly data from the Center for Space Research at the University of Texas (CSR) were used: Level-2 Release 05 (CSR05), Level-2 Release 06 (CSR06), and a recent mascon solution based on CSR06 (CSR06-mascon). The first two versions of the relative gravity field, in the form of 60 degree and order of spherical harmonic (SH) coefficients (i.e., GSMs), were accessible at http://icgem.gfz-pots- Dashed boundary (on the right) is our regional study area with selected (red dots) and unselected (grey dots) isotope samplings, along with their sampling codes (above the dots) and age (in month) (below the dots) (adapted with permission from ref. [9], Copyright 2010 American Geophysical Union).
The SCS is a western marginal sea of the Pacific ocean with significant water exchange through the Karimata, Malacca and Luzou straits [34]. Between June and September every year, the predominant northwestward along-shore current affects the tropical river plume extension [4]. Estuarine mixtures are transported to the western SCS when Mekong runoff normally peaks in September. Besides this, the well-known eastward South Vietnam upwelling at about 12 • N and frequent eddies in summer further complicate the sea surface variability [35,36]. The global El Niño-Southern Oscillation (ENSO) events, which are closely associated with the interannual variability of tropical monsoons, also cause anomalous changes in both the Mekong basin runoff and oceanic dynamics in this geographic region [37,38]. As such, Mekong freshwater transport in the western SCS will be affected by the above combined effects.

Oceanic EWH Time Series from GRACE
In this study, three versions of GRACE monthly data from the Center for Space Research at the University of Texas (CSR) were used: Level-2 Release 05 (CSR05), Level-2 Release 06 (CSR06), and a recent mascon solution based on CSR06 (CSR06-mascon). The first two versions of the relative gravity field, in the form of 60 degree and order of spherical harmonic (SH) coefficients (i.e., GSMs), were accessible at http://icgem.gfz-potsdam.de/ series (accessed on 20 March 2021). The last one used the recently employed regional mass concentration functions (mascons) to parameterize the Earth's gravity field [39]. The CSR06-mascon product with all relevant corrections in a 0.25 • × 0.25 • grid format could be Remote Sens. 2021, 13, 1193 4 of 18 directly used and was accessible at http://www2.csr.utexas.edu/grace/ (accessed on 20 March 2021). [40]. All data products were converted into EWH for our investigation.
Note that GSM represents geopotential coefficients after corrections to errors of commission or omission from the background gravity models [41], in which they were corrected by a priori models of solid and ocean tides, pole tides, non-tidal atmospheric and oceanic loadings, etc. Customarily, to study the oceanic mass variation, part of the non-tidal atmospheric and oceanic loading from the Atmosphere and Ocean De-Aliasing Level-1B (AOD1B)-modeled data product (i.e., the GAD component) was re-added back to the GSM in order to restore the ocean bottom pressure [42] or mass-induced sea level variation [26]. Nonetheless, the AOD1B-modeled data did not account for river flux into the ocean [43]. Therefore, to study the freshwater mass variations in the ocean surface caused by river flux in our present investigation, we directly used GSM.
The following steps ( Figure 2) were employed to derive the EWH time series for GRACE Level-2 CSR05 and CSR06 data products: (1) C 20 coefficients were replaced with that observed via satellite laser ranging [44]; (2) geocenter motion was added to the degree-1 SH coefficient [45]; (3) the SH coefficients of the GGM05C gravity field [46] was subtracted; (4) a de-striping procedure and 350 km Gaussian filtering were applied to reduce the combined noises in the high degrees of the SH coefficients [47]. Notice that the CSR06mascon EWH was free from step (4), because the compact geodesic grid and the design of the regularization matrix had already taken into consideration the conditions of land and ocean separately [39].
Remote Sens. 2021, 13, 1193 4 of 18 dam.de/series (accessed on 20 March 2021). The last one used the recently employed regional mass concentration functions (mascons) to parameterize the Earth's gravity field [39]. The CSR06-mascon product with all relevant corrections in a 0.25° × 0.25°grid format could be directly used and was accessible at http://www2.csr.utexas.edu/grace/ (accessed on 20 March 2021). [40]. All data products were converted into EWH for our investigation. Note that GSM represents geopotential coefficients after corrections to errors of commission or omission from the background gravity models [41], in which they were corrected by a priori models of solid and ocean tides, pole tides, non-tidal atmospheric and oceanic loadings, etc. Customarily, to study the oceanic mass variation, part of the nontidal atmospheric and oceanic loading from the Atmosphere and Ocean De-Aliasing Level-1B (AOD1B)-modeled data product (i.e., the GAD component) was re-added back to the GSM in order to restore the ocean bottom pressure [42] or mass-induced sea level variation [26]. Nonetheless, the AOD1B-modeled data did not account for river flux into the ocean [43]. Therefore, to study the freshwater mass variations in the ocean surface caused by river flux in our present investigation, we directly used GSM.
The following steps ( Figure 2) were employed to derive the EWH time series for GRACE Level-2 CSR05 and CSR06 data products: (1) C20 coefficients were replaced with that observed via satellite laser ranging [44]; (2) geocenter motion was added to the degree-1 SH coefficient [45]; (3) the SH coefficients of the GGM05C gravity field [46] was subtracted; (4) a de-striping procedure and 350 km Gaussian filtering were applied to reduce the combined noises in the high degrees of the SH coefficients [47]. Notice that the CSR06-mascon EWH was free from step (4), because the compact geodesic grid and the design of the regularization matrix had already taken into consideration the conditions of land and ocean separately [39].  Moreover, a linear model with the offset, trend, and periodic terms was used to fill the missing monthly data. To avoid the potential seismic signal contamination before and after the nearby Sumatra-Andaman earthquake in December 2004 [48], the data spanning from 2005 to 2012 were used with missing data in January 2011, June 2011, August 2012 and October 2012, interpolated through the fitting equation (2) below. GRACE data after 2012 were not selected because at least three monthly data sets were missing for each year. All the EWHs were resampled into 1 • × 1 • grids. Following the above procedures, the average gridded EWH time series over 9-15 • N and 110-114 • E in the western SCS (blue regions in Figure 1) are presented in Figure 3. Among the three time series, the EWHs generated from CSR05 and CSR06 matched each other well, but the EWH generated from the CSR06-mascon presented a much lower amplitude (i.e.,~1.5 cm) than those from CSR05 and CSR06 (i.e.,~5 cm). This difference depends on the completely different inversion methodology used for inferring EWHs, in addition to the refined geophysical corrections.
Moreover, a linear model with the offset, trend, and periodic terms was used to fill the missing monthly data. To avoid the potential seismic signal contamination before and after the nearby Sumatra-Andaman earthquake in December 2004 [48], the data spanning from 2005 to 2012 were used with missing data in January 2011, June 2011, August 2012 and October 2012, interpolated through the fitting equation (2) below. GRACE data after 2012 were not selected because at least three monthly data sets were missing for each year. All the EWHs were resampled into 1° × 1° grids. Following the above procedures, the average gridded EWH time series over 9-15°N and 110-114°E in the western SCS (blue regions in Figure 1) are presented in Figure 3. Among the three time series, the EWHs generated from CSR05 and CSR06 matched each other well, but the EWH generated from the CSR06-mascon presented a much lower amplitude (i.e., ~1.5 cm) than those from CSR05 and CSR06 (i.e., ~5 cm). This difference depends on the completely different inversion methodology used for inferring EWHs, in addition to the refined geophysical corrections.

Mekong Basin Runoff from In Situ Hydrological Stations
Along the two main tributaries of the Mekong estuary mouth, the Chau Doc and Tan Chau hydrological stations are about 200-300 km away from both the inland Tonle Sap Lake and the estuary mouth ( Figure 1). Therefore, the effects of lake regulation and the backwater due to ocean tides on the runoff were minimized [49]. Note that we assumed that sub-surface discharge was negligible when compared to the main flow. Then, the observed discharge time series of the two stations were summed up in order to represent the whole basin runoff and average out the ocean tidal backwater effect due to semidiurnal and diurnal ocean tides. Both in situ discharge time series can be requested from the Mekong River Commission (MRC) (http://www.mrcmekong.org accessed on 20 March 2021).

Mekong Basin Runoff from In Situ Hydrological Stations
Along the two main tributaries of the Mekong estuary mouth, the Chau Doc and Tan Chau hydrological stations are about 200-300 km away from both the inland Tonle Sap Lake and the estuary mouth ( Figure 1). Therefore, the effects of lake regulation and the backwater due to ocean tides on the runoff were minimized [49]. Note that we assumed that sub-surface discharge was negligible when compared to the main flow. Then, the observed discharge time series of the two stations were summed up in order to represent the whole basin runoff and average out the ocean tidal backwater effect due to semidiurnal and diurnal ocean tides. Both in situ discharge time series can be requested from the Mekong River Commission (MRC) (http://www.mrcmekong.org accessed on 20 March 2021).
To be consistent with the GRACE data time span, the daily discharge time series spanning between January 2005 and December 2012 were selected. To convert the unit of the daily discharge (in m 3 /s) into runoff (mm), it was divided by the total Mekong river basin drainage area (i.e., 795,000 km 2 ), while adding up daily values of each month. In general, the Mekong basin runoff generally reaches its maximum and minimum in September and April, respectively. During the strong La Niña event during 2010-2011, the runoff reached the lowest peak value at~76 mm in 2010 and the highest peak value of over 100 mm in 2011 ( Figure 4). Besides this, to be consistent with the GRACE EWH anomalies, runoff time series were subtracted from the mean of the entire study period. the daily discharge (in m /s) into runoff (mm), it was divided by the total Mekong river basin drainage area (i.e., 795,000 km 2 ), while adding up daily values of each month. In general, the Mekong basin runoff generally reaches its maximum and minimum in September and April, respectively. During the strong La Niña event during 2010-2011, the runoff reached the lowest peak value at ~76 mm in 2010 and the highest peak value of over 100 mm in 2011 ( Figure 4). Besides this, to be consistent with the GRACE EWH anomalies, runoff time series were subtracted from the mean of the entire study period.

Isotope-Derived Age Based on Radium Isotope Data Measured in 2007
Ref. [9] investigated the distributions of the 223 Ra, 226 Ra, and 228 Ra isotopes, and salinity variations, on the ocean surface of the western SCS during a one-month field trip between August and September 2007. Based on an empirical two-end-member mixing model [6,50], the apparent age of the Mekong river plume freshwater was determined (at dots shown in Figure 1), serving as a ground truth to validate our GRACE-derived time lag. Note that the ground truths are located more than 800 km away from estuary mouth. Hence, the land-to-ocean leakage effects of the GRACE data should be negligible. The isotope-derived apparent age (t) of each sample can be calculated by: where and represent 223 Ra activity (in disintegrations per minute (dpm)/100 L) and the dimensionless estuarine fraction of each sampling (modeled with an input of the 228 Ra and 226 Ra activities and salinity), respectively; represents the 223 Ra activity of estuarine end-member sampling (with a default constant of 0.74 dpm/100 L), and

Isotope-Derived Age Based on Radium Isotope Data Measured in 2007
Ref. [9] investigated the distributions of the 223 Ra, 226 Ra, and 228 Ra isotopes, and salinity variations, on the ocean surface of the western SCS during a one-month field trip between August and September 2007. Based on an empirical two-end-member mixing model [6,50], the apparent age of the Mekong river plume freshwater was determined (at dots shown in Figure 1), serving as a ground truth to validate our GRACE-derived time lag. Note that the ground truths are located more than 800 km away from estuary mouth. Hence, the land-to-ocean leakage effects of the GRACE data should be negligible. The isotope-derived apparent age (t) of each sample can be calculated by: where 223 Ra obs and f es represent 223 Ra activity (in disintegrations per minute (dpm)/100 L) and the dimensionless estuarine fraction of each sampling (modeled with an input of the 228 Ra and 226 Ra activities and salinity), respectively; 223 Ra es represents the 223 Ra activity of estuarine end-member sampling (with a default constant of 0.74 dpm/100 L), and λ 223 is the decay constant of 223 Ra, which was 0.061 d −1 . The age, t, of the Mekong river plume samples is depicted in Figure 1.
Considering the high uncertainty and sensitivity resulting from the end-member coefficients in the model, extremely anomalous samplings were not chosen where (1) the contribution of the Mekong river plume to the surface water, f es , was less than 0.2 or greater than 1.0, and (2) the calculated age was less than 0.2 months (i.e., 6 days). Under these criteria, 17 out of 25 isotope samples remained (at red dots shown in Figure 1), whose derived ages would be further compared to GRACE-derived time lags. All Mekong river plume sample positions, f es , and ages are listed in Table A1, with unselected sample codes and corresponding extreme values underlined.

Methodology
Time lag analysis was conducted to determine the time lag between in situ runoff near the estuary mouth and the freshwater component of the GRACE EWH time series on the ocean surface in our study region. Two methods were employed. One is a harmonic analysis to determine the annual phase of in situ runoff and GRACE EWH time series, followed by taking the phase difference between them (hereinafter called the P-method); the other is performing the cross-correlation analysis via an empirical shift in GRACE EWH time series against the in situ runoff in our study region until reaching maximum cross-correlation (hereinafter called the C-method). Each method is independent of the other, except that the latter one used the detrended time series generated from the former one for calculating cross-correlation.

Phase Analysis (P-Method)
To obtain the phase of semi-annual and annual signals of EWH, a conventional harmonic analysis is employed, which is expressed as: where a and b are the offset and the linear trend; A i and ϕ i represent the corresponding amplitude and phase with predefined period T i for annual (i = 1) and semiannual (i = 2) signals, respectively. Six parameters (i.e., a, b, A 1 , A 2 , ϕ 1 , and ϕ 2 ) are to be determined via least-squares estimation. The same model has also been used for monthly in situ runoff and filling the missing data of monthly EWH. According to the regional averaged result from the western SCS from 2005 to 2012 (Table A2), the annual amplitude is several times larger than the semiannual one, in terms of both runoff and EWHs. Ignoring the semi-annual one, the freshwater time lag from the Mekong estuary (specifically at two hydrological stations) to our study region could be calculated as the annual phase (ϕ 1 ) difference between EWH and runoff.
Considering the variance calculated for the lower frequency, at least three complete cycles of data should be used in order to be more representative in terms of statistics [51]. We also evaluated the subsequent difference in our results by sequentially adding one more year of data span for data fitting in our study area starting from 2005 [52] (See Appendix B). We found that the difference (in terms of root-mean-square error (RMSE)) was stabilized when a three-year data time span was employed for all three GRACE time series. The RMSE converges to less than 10 mm/year ( Figure A1). Therefore, a three-year data time span was adopted as a data window rolling month-by-month during 2005-2012.

Cross-Correlation Analysis (C-Method)
The cross-correlation between EWH and runoff time series allows us to empirically determine the spatial distribution of time lag, representing freshwater transport on the ocean surface. The cross-correlation is, in fact, a Pearson correlation coefficient (PCC) introduced with a time shift, δτ, expressed as: where R i and E i+δτ are the runoff and EWHs at month i and month i + δτ, respectively, and R and E are their mean values. Finally, when the PCC(δτ) reaches the maximum, the corresponding δτ is considered as the optimal freshwater time lag. In this method, the monthly detrended runoff and EWHs time series were interpolated to a daily interval based on a bi-cubic interpolation algorithm in MATLAB. Note that interpolation is necessary for daily δτ despite the fact that this process would introduce unavoidable uncertainty into our results. In addition, the cross-correlation is less dependent on the length of time series. Note that whenever time shift δτ is applied, the number of valid fitting months in the formula, N, is slightly less than the given time period. Therefore, the three-year rolling window method was still applied in order to be consistent with the P-method.

Regional Study Results by Phase Analysis (P-Method)
Employing the three-year rolling window, the three-year trend, annual amplitude and annual phase variation in runoff and EWHs are calculated for each three-year period during 2005-2012 ( Figure 5). The time lag is derived by the annual phase difference between EWHs and in situ runoff. It should be noted that the values of CSR06-mascon annual amplitude are multiplied by two for a better visual expression of the annual amplitude in Figure 5b. Table 1 displays the PCC between runoff and EWH time series in terms of three-year trend, annual amplitude and annual phase variation, along with the calculated mean lag and corresponding standard deviation (STD). Within the time span shown in Figure 5b, the annual amplitude of runoff ranges from 36.15 to 42.28 mm, with the mean annual amplitude of ~49.9 mm and ~61.0 mm for CSR05 and CSR06, respectively, which values are approximately four times larger than that of  During 2005-2008, the runoff and the CSR05 and CSR06 EWH trends increase slowly. However, no apparent change manifests for the CSR06-mascon EWH (Figure 5a). The threeyear trend in EWHs and runoff rapidly decline when the rolling window starts to cover the time series of 2008 and 2009, reaching a minimum of~-0.5 to -0.8 mm/month in 2009. This abrupt downward trend is speculated to be attributable to the moderate El Niño event during 2009-2010, which reduced seasonal rainfall, storage and runoff in the Mekong basin, as concluded from [37,53], thus weakening the Mekong river plume and its transport on the ocean surface. After this anomaly, the three-year trend in runoff and EWHs gradually increases and reaches the highest values, 0.5 mm/month and 1 mm/month, respectively. Among the three decomposed components of the EWH compared to those of the runoff, the highest PCC of 0.92 is attained for the CSR06 EWH in terms of three-year trend.
Within the time span shown in Figure 5b, the annual amplitude of runoff ranges from 36.15 to 42.28 mm, with the mean annual amplitude of~49.9 mm and~61.0 mm for CSR05 and CSR06, respectively, which values are approximately four times larger than that of CSR06-mascon (i.e.,~12.7 mm). These numerical differences in GRACE are mostly due to the different versions of the non-tidal oceanic and atmospheric model (i.e., AOD1B), corrected for CSR05 and CSR06, and the design of the regularization matrices for continental and oceanic areas for CSR06-mascon. The peaks and troughs of three EWHs are about 1~2 months later than the runoff time series. The PCCs between EWH and runoff in terms of annual amplitude are all higher than 0.60 (Table 1), and CSR05 EWH attains the highest PCC of 0.78.
On the contrary, the annual phase variations for the three EWHs are shown to be inconsistent with that of runoff, particularly in the ENSO event in 2010. During the El Niño event, the annual phase variation of the runoff was relatively small (Figure 5c), indicating the time lag determination by the annual phase difference would depend more on the EWH time series themselves, rather than on the in situ runoff. In general, the PCCs between all the three EWHs and the in situ runoff in terms of the annual phase were relatively lower than that of the three-year trend variation and annual amplitude (Table 1), showing the instability of the annual phase variation in freshwater mass as derived from GRACE in the open ocean. This inconsistency potentially reflects the time-varying influence of the ocean on Mekong freshwater transport. The mean and STD of the regional time lag are calculated in the last column of Table 1. Overall, the time lag estimated from CSR05 had the longest mean value (20 days) and the lowest STD (3.1 days).

Regional Study Results by Cross-Correlation Analysis (C-Method)
Using the C-method, the freshwater transport from estuary to our study region in the ocean has been determined by a proper δτ while maximizing the PCC(δτ) value. Again, employing the three-year rolling window, the original PCC with no lag (hereafter called oPCC), and the improved PCC (hereinafter called iPCC) with the proper time lag, δτ, are shown in the three EWH time series for each three-year period ( Figure 6). Their corresponding mean values and STD are listed in Table 2. among which CSR06 yields the best correlation with runoff, with its iPCC reaching 0.96. The mean values of time lag are slightly smaller than that from the P-method, along with a smaller STD (i.e., ~2 days). Besides this, the time lag differences among the three EWHs are not apparent, with the highest time lag of 16 days for CSR05 and the lowest one of 14 days for CSR06-mascon. During the entire period, the mean and STD of the lag (Table 2) are all lower than those retrieved by the P-method (Table 1), which implies that the estimation by C-method was temporally smoother.   The oPCC values for three EWHs are all above 0.75 during the time span in our study, indicating their high relationship despite the time lag existing between EWHs and runoff. The mean improvement from oPCC to iPCC is~0.05 in terms of three EWHs (Table 2), among which CSR06 yields the best correlation with runoff, with its iPCC reaching 0.96. The mean values of time lag are slightly smaller than that from the P-method, along with a smaller STD (i.e.,~2 days). Besides this, the time lag differences among the three EWHs are not apparent, with the highest time lag of 16 days for CSR05 and the lowest one of 14 days for CSR06-mascon. During the entire period, the mean and STD of the lag (Table 2) are all lower than those retrieved by the P-method (Table 1), which implies that the estimation by C-method was temporally smoother.

Discussion
Former isotope-derived ages in the western SCS measured between August and September 2007 [9] were employed to evaluate our time lag results for the two presented methods in the same region.

Gridded GRACE-Derived Time Lag Estimation by Two Methods
For evaluation, the time lags were determined using the P-method and the C-method in a three-year rolling window, with the median time epoch at 1 September 2007. Figure 7 presents the time lag results (in month) of each individual grid in the western SCS for the three EWHs derived by the P-method (Figure 7a-c) and the C-method (Figure 7d-f).
the lag between runoff and the three EWHs.

Discussion
Former isotope-derived ages in the western SCS measured between August and September 2007 [9] were employed to evaluate our time lag results for the two presented methods in the same region.

Gridded GRACE-Derived Time Lag Estimation by Two Methods
For evaluation, the time lags were determined using the P-method and the C-method in a three-year rolling window, with the median time epoch at 1 September 2007. Figure  7 presents the time lag results (in month) of each individual grid in the western SCS for the three EWHs derived by the P-method (Figure 7a-c) and the C-method (Figure 7d-f).  Among the six results, the gridded time lag from CSR06 using the P-method (Figure 7b) shows a clear eastward freshwater flowing path from the Mekong estuary along the Vietnam coastline to the ocean. The propagation was the integrated results of the seasonal offshore upper current [54,55], the monsoon [56], and ENSO events [57]. For the GRACE products, the gridded time lag results from CSR05 and CSR06 (Figure 7a,b,d,e) display several longitudinal panels parallel with the coastline, which were not evident from CSR06mascon (Figure 7c,f). The existence of filtering for the Level-2 GRACE data might contribute to these differences. In terms of methodology, the spatial distribution of time lag using the C-method was smoother than that using the P-method, in addition to there being a narrower range of determined time lag. This indicates that the spatial distribution of the time lag retrieved from the EWHs and runoff time series depends not only on the GRACE data, but also on the methodologies themselves. Figure 8 illustrates the difference between the isotope-derived age in September 2007 and the GRACE-derived time lag in the corresponding time span by the two presented methods. The time lag at the boundary of the GRACE grids is interpolated. Based on the scatter plots of time lag versus age, the time lag determined using the C-method fits the isotope-derived age better than that derived using the P-method, partially because the former method yields a narrower range of time lag. This result is particularly evident for the results from CSR06-mascon (Figure 8f), in which the lag clusters were approximately 0.5 months. However, in this case, the higher the age, the larger the discrepancy between the GRACE-derived time lag and the CSR06-mascon.

Comparison between GRACE-Derived Time Lag and Isotope-Derived Age
Remote Sens. 2021, 13, 1193 13 of 18 Figure 8. Difference between GRACE-derived time lag and isotope-derived age for 17 selected isotope samplings. Note that the GRACE-derived time lag was calculated from CSR05, CSR06 and CSR06-mascon using the P-method (a-c) and the C-method (d-f), respectively. The subgraphs (bottom left) plot the scatterings of the corresponding time lag and age from the different versions of GRACE data, derived using the P-method and the C-method.  . Difference between GRACE-derived time lag and isotope-derived age for 17 selected isotope samplings. Note that the GRACE-derived time lag was calculated from CSR05, CSR06 and CSR06-mascon using the P-method (a-c) and the C-method (d-f), respectively. The subgraphs (bottom left) plot the scatterings of the corresponding time lag and age from the different versions of GRACE data, derived using the P-method and the C-method.
The differences between the GRACE-derived time lag and the isotope-derived age are displayed in colored dots at 17 selected stations (Figure 8). When the residual is larger than 0.2 months (i.e., 6 days) or smaller than -0.2 months, we consider the error was positively or negatively large, respectively. Based on this criterion, samplings with negatively large errors were observed mostly at locations away from the estuary (e.g., Y05, Y06, Y35), and so too were those with rather positively large errors (e.g., 2Y94, Y26, Y98). The high uncertainty of the end-member model used for calculating age by radium isotope partially contributes to these extremes. Besides this, a uniform time lag estimation by long-term time series analysis (i.e., three-year rolling window) is unable to detect short-term anomalies in the ocean. One example is the sampling station Y12, which displays an extremely old age because an eddy occurred nearby when the radium isotope was measured [9].
The mean and STD of age and time lag corresponding to all isotope samplings are listed in Table 3, in terms of GRACE data products and methods. Among the three GRACE data products, the mean of time lag derived from CSR05 approximates the mean of the isotope-derived age the most, with the differences of 0.03 months (i.e., 1 day) and -0.10 months (i.e.,~3 days) for the P-method and C-method, respectively. In terms of the methodology, the mean of time lag given by the P-method was slightly larger than that given by the C-method, which is consistent with our regional study results in Section 4.1. Additionally, the STD of the time lag minus age is~0.3 months (i.e., 9 days), indicating that the quantitative comparison between time lag and age is greatly influenced by extreme values, partly because of the high uncertainty in the end-member model. Despite all these differences, the GRACE-derived time lag still fits properly with most isotope-derived ages, where the absolute errors are less than 0.2 months (i.e., 6 days). The paired t-test (Table A3) does not reject the null hypothesis of their difference equaling zero at the 1% significance level (except results given by the C-method for CSR06-mascon). According to measurements from Argo profiling floats during 2007, the mean surface velocity of the along-shore current was about 49 cm/s [58]. This means estuarine freshwater following this path would take about half a month to reach the center of our study area (i.e., 700 km away from the estuary). This estimation was largely shown to be in agreement with our results and isotope-derived ages.
Overall, the GRACE-derived time lags are slightly smaller than the isotopic tracer ages. This could be partially caused by the deficiencies of our methods. For instance, the P-method depends on the length of data time span used for phase determination, whereas the C-method requires interpolation to the daily interval for the time shift in the cross-correlation analysis. The ages inferred from the isotope tracer analysis also suffer from the in situ measurement uncertainties, which are related to nearby ocean eddies and the highly sensitive end-member value determination process in the mixing model derived from [6]. Nonetheless, the results of both the methods presented in this study are more stable than those from the isotopic tracer analysis.

Conclusions
In this paper, the in situ runoff and GRACE product-derived freshwater mass in the ocean (EWHs) were used to determine the freshwater transport, in terms of time lag, from the Mekong river estuary to the western SCS via two independent methods: phase analysis (i.e., P-method) and cross-correlation analysis (i.e., C-method). A rolling data window was further applied to determine the time lag variation every three years during the period 2005-2012.
Among the three components (i.e., three-year trend, annual amplitude and phase) modeled in the P-method, EWHs correlated with runoff well in terms of trend and annual amplitude, with the minimum Pearson correlation coefficient (PCC) averaged for all threeyear windows reaching 0.71 and 0.64, respectively. Inconsistency in the annual phase between runoff and EWHs potentially reflects the ocean's influence on freshwater transport. Besides this, the differences resulting from three GRACE data products are apparent when climate extremes are encountered, for instance, a moderately strong El Niño event during the period 2009-2010. For the P-method, the average determined time lag of our study region was about 20, 17 and 15 days for CSR05, CSR06 and CSR06-mascon, respectively. In terms of the C-method, the mean PCC improvement was~0.05 for three EWHs, and the determined time lag and its STD were relatively smaller compared with that derived from the P-method. The EWH value that correlated the best with runoff was given by CSR06, in which the mean PCC was improved from 0.91 to 0.96. Compared to the P-method, the fluctuation in time lag retrieved from the C-method was slightly lower in terms of both mean values and STD, thus weakening the difference caused by GRACE products.
To evaluate these two methods, the age results derived from radium isotope samplings in the western SCS by Chen et al. [9] in September 2007 were compared with our gridded time lag estimates whose median epoch in the rolling window was in the same month. Among the two methods, the time lag derived from the P-method showed a more distinguishable spatiotemporal pattern of freshwater propagation from the Mekong estuary toward the east, in particular CSR06. For the quantitative comparison against the GRACE-derived time lag, the isotope-derived age ranges were wider, showing a relatively high STD. This is partly due to high-frequency events (e.g., eddies) influencing isotope measurements, which cannot be recorded by GRACE. However, the isotope-derived age and GRACE-derived time lag still present a considerable agreement with each other, with over half the samplings showing small differences (i.e.,~6 days), regardless of the method. Overall, time lag determination by GRACE demonstrates a good potential for monitoring estuarine freshwater transport in the ocean. However, more experiments in different estuaries are required to further validate the universal feasibility of our presented analysis.  Appendix A Table A1. The information of radium isotope sampling for the Mekong river plume group adapted with permission from ref. [9], Copyright 2010 American Geophysical Union.

Appendix B
The root-mean-square error (RMSE) for monthly EWH time series from 2005 to Y is calculated by [52]: where N is the number of grids in the study area, and v Y n is the fitted long-term trend of EWHs from 2005 to Y for each grid. Note that Y is up to 2011 because the valid GRACE time series in this study span from 2005 to 2012. Figure A1 shows the RMSE against the increasing data fitting period.

Appendix B
The root-mean-square error (RMSE) for monthly EWH time series from 2005 to Y is calculated by [52]: 1 where N is the number of grids in the study area, and is the fitted long-term trend of EWHs from 2005 to Y for each grid. Note that Y is up to 2011 because the valid GRACE time series in this study span from 2005 to 2012. Figure A1 shows the RMSE against the increasing data fitting period. Figure A1. The variance in RMSE when the fitting period ranges from 1 year to 7 years.