Radar Altimetry as a Proxy for Determining Terrestrial Water Storage Variability in Tropical Basins

: The Gravity Recovery and Climate Experiment (GRACE) mission has provided us with unforeseen information on terrestrial water-storage (TWS) variability, contributing to our understanding of global hydrological processes, including hydrological extreme events and anthropogenic impacts on water storage. Attempts to decompose GRACE-based TWS signals into its di ﬀ erent water storage layers, i.e., surface water storage (SWS), soil moisture, groundwater and snow, have shown that SWS is a principal component, particularly in the tropics, where major rivers ﬂow over arid regions at high latitudes. Here, we demonstrate that water levels, measured with radar altimeters at a limited number of locations, can be used to reconstruct gridded GRACE-based TWS signals in the Amazon basin, at spatial resolutions ranging from 0.5 to 3 ◦ , with mean absolute errors (MAE) as low as 2.5 cm and correlations as high as 0.98. We show that, at 3 ◦ spatial resolution, spatially-distributed TWS time series can be precisely reconstructed with as few as 41 water-level time series located within the basin. The proposed approach is competitive when compared to existing TWS estimates derived from physically based and computationally expensive methods. Also, a validation experiment indicates that TWS estimates can be extrapolated to periods beyond that of the model regression with low errors. The approach is robust, based on regression models and interpolation techniques, and o ﬀ ers a new possibility to reproduce spatially and temporally distributed TWS that could be used to ﬁll inter-mission gaps and to extend GRACE-based TWS time series beyond its timespan.


Introduction
The Gravity Recovery and Climate Experiment (GRACE) satellites [1] have provided unprecedented information on global water-storage dynamics by estimating variations in the terrestrial water-storage (TWS) variability. GRACE data has been used to evaluate global TWS trends [2], to monitor and characterize droughts (e.g., [3][4][5], estimate groundwater depletion due to anthropogenic activities [6][7][8][9] and improve hydrological models through data assimilation techniques (e.g., [10,11]) towards a better quantification of land water-storage components and their spatial and temporal variability. Since GRACE data provides estimates of the variability of the total amount of water-storage change within a domain, the inability to identify the water layer and state has been a challenge in decomposing TWS signals into its different components (i.e., groundwater, soil moisture, snow and surface water (including rivers, floodplains, wetlands, lakes and reservoirs). Few studies attempted to determine the impact of surface water storage (SWS) on TWS at the basin scale (e.g., [5], [6]). Using state-of-the-art hydrological models, [12] were able to decompose TWS and determine the  [42] Reduced optimal interpolation using GRACE data and in situ river level records Amazon basin 1980-2008 Yes [18] Water budget using GRACE data and model outputs Amazon basin 1948-2010 No [20] Statistical model using reanalysis data Global 1985-2015 Yes [43] Physically based modelling and deep learning India 2002-2017 Yes

Gravity Recovery and Climate Experiment (GRACE) Terrestrial Water Storage (TWS)
GRACE measurements of the Earth's time-variable gravity are processed by different centers that provide monthly gravity solutions using two approaches: spherical harmonics and mass concentration (mascons) blocks [44,45]. During the first years of GRACE observations, spherical harmonics was the standard. However, recent studies have shown that the mascons solution offers several key advantages, including (1) reduced leakage through the land-ocean interface, which increases signal amplitude; (2) removal of north-south stripes is achieved without the need of empirical filters; and (3) finer resolution for smaller areas [45,46]. Here, we used the University of Texas' Center for Space Research (CSR) version RL05 [47]. Although CSR GRACE masons are natively estimated at 120 km-wide mascon blocks and resampled to a 0.5 • spatial resolution grid, the actual resolution of its solution is about 250-300 km along the equator line. Mascons solutions from other GRACE data-processing centers have different native resolution (e.g., three degree equal-area caps). Therefore, our analysis considered four grid sizes: 0.5, 1, 2 and 3 degrees. For this reason, TWS was resampled into coarser resolutions by averaging the 0.5 • data. GRACE data have multiple sources of errors that emanate from instruments measuring through post-processing [48]. A number of attempts to assess GRACE estimates have been undertaken (e.g., [49,50]). GRACE data accuracy has been estimated to be less than 1cm in equivalent water height, when averaged over areas larger than about 4 × 105 km 2 , and errors increase as the area under observation decreases [51]. [49] noted that the errors increase at lower latitudes (~2.6 cm), whereas it is around 0.8 cm near the poles. While evaluating the global CSR RL05 mascons solution, [46] noted that GRACE-based TWS errors were the lowest (<3cm) for large basins (>0.5 million km 2 ), but up to three times higher for small basins.

Radar Altimetry Data
A total of 984 virtual stations are available over the Amazon basin for 2002-2016 on the Theia-Land's Hydroweb platform [29], derived from multiple satellite missions: Jason (11 VS), Jason 2 (88 VS), Envisat (741 VS) and Sentinel 3A (144). The temporal resolution varies as a function of the satellite, i.e., 10 days for Jason, 27 days for Sentinel 3A and 35 days for Envisat. [28] list two main sources of errors in RA data. The first can be due to bounced back radar echoes, captured by the antenna receiver that are caused by solid reflectors (e.g., vegetation, islets, etc). When such echoes are combined with those from water surfaces at a different level, the height inference will be affected. Second, RA retrievals can be biased if the detected surface is off-nadir, because the signal received by the satellite may be dominated by water surface at the edge of the footprint. A tool was developed by [28] to overcome those limitations and identified errors in altimetric time series below 20 cm up to meters.
More specifically, they reported a characteristic error in the Envisat series below 30 cm. More details about the generation of altimetry-based water levels used in this study can be found in [28].

Estimating Terrestrial Water Storage from Radar Altimetry
In this study, water level data was filtered by removing VS containing less than 12 months of data. Gaps in the VS time series were filled using a shape-preserving piecewise cubic interpolation [52]. This method requires at least four neighboring grid points for the interpolation to be performed. To reduce the chances of generating unrealistic interpolations, we only interpolated time series whose gaps do not exceed 30% of its length. For example, a time series spanning 24 months was gap filled only if the number of valid records exceeded 17, i.e., a maximum of 7 gaps ( Figure S1 shows an example of RA time-series interpolation using this technique). Gaps in the tails of the time series were not filled, i.e., we did not perform temporal extrapolation. The filtering process resulted in 839 VS, whose availability is unevenly distributed through time, i.e., a fraction of that figure is available each month ( Figure S2 shows the temporal distribution of VS data availability). Radar altimetry data was converted to the monthly time step to match GRACE-based TWS estimates (TWS GC ) by averaging observations within the same month. GRACE pixels corresponding to each VS were identified and correlations were computed between water levels and TWS. Because some RA-based water level anomalies correlate better with TWS than others, we selected virtual stations with highest correlation (Figure 1). This procedure was repeated for each of the four GRACE resolutions considered in this study. We then used the best correlated points to fit a linear (TWS = a·H + b, where H stands for water level anomalies) and a quadratic (TWS = a·H 2 + b·H + c) regression models that estimate TWS from radar altimetry (Table S1 provides regression equations and corresponding metrics for the 3 • solution). For each GRACE pixel, the regression model with the lowest mean square error was selected. The regression models obtained for the calibration period were used to derive the estimates for the validation period, i.e., during validation, the procedures described above are not repeated. As a result, the number of pixels with regression-based TWS (TWS reg ) estimates varied according to the spatial resolution and month (i.e., finer resolutions result in more pixels covering the domain, hence more VS are used, one for each pixel, wherever data is available) ( Figure S3). The average number of GRACE pixels with water-level data for each resolution was: 388 (0.5 • ), 216 (1 • ), 81 (2 • ) and 41 (3 • ). Note that TWS reg estimates are limited to the pixels containing at least one VS. The remaining pixels do not have RA data with which a correlation could be established. To generate interpolation-based TWS estimates (TWS interp ) in those remaining pixels, three interpolation methods were tested: ordinary kriging (krig), inverse distance weighting (IDW), and nearest neighbor (NN). The method with best performance was applied to produce the final TWS fields. Kriging is a geostatistical method based on the regionalized variable theory, which accounts for the spatial variability of natural properties by means of a stochastic approach [53]. It has been widely used to perform spatial analysis of meteorological and hydrological data (e.g., [54][55][56][57]). Most of these data often present a spatial dependence, i.e., similarities between two points are likely to increase as the distance between them decrease. Thus, one can infer how a given property varies in space between sample points based on an experimental variogram, which describes the spatial variability of a regionalized property or variable.
Here, we used the time-averaged TWS from GRACE to derive a unique experimental variogram for each spatial resolution, from which weighting factors (λ i ), associated with TWS reg at location x i , are derived. The kriging estimate TWS interp at location x 0 is given by: In the IDW method, TWS interp is the weighted average of the known neighbor values; the weights (w i ) are computed as: where d is the distance between the x i and x 0 . Both krig and IDW methods were applied using the R package gstat [58]. The NN method uses the Voronoi tessellation of a limited number of sample points to select, for each grid point (x o ), neighboring data points for interpolation [59]. Weighting is based on the influence area of such data points. We used Matlab's internal 'scatteredInterpolant class' to perform the NN interpolation. Previous studies have applied IDW and NN to (spatially) interpolate rainfall and other environmental data at multiple time scales (e.g., [60][61][62][63][64][65]).
Finally, the different TWS interp estimates were merged with TWS reg to produce the final altimetry-based TWS (TWS alt ). The interpolation method with best performance was applied to produce the final TWS fields.

Evaluation Procedure
For both calibration and validation periods, the mean absolute error (MAE), correlation coefficient (CC) and standard deviation ratio (γ) between our TWS estimates and TWS GC were calculated.
where TWS est can be TWS reg , TWS interp or TWS alt ; σ GC and σ est are, respectively, standard deviation of TWS GC and TWS est time series. We also evaluated the final product during extreme events and by comparison against other products. Regarding the former, we analyzed TWS GC seasonal anomalies between 2003 and 2010; and selected two extreme periods (one wet and one dry) to carry out such comparison. We then compared TWS estimates averaged for the entire Amazon basin and its main tributaries: Solimões, Negro, Tajapós, Madeira, Xingu rivers, and the area where the extreme events occurred.
Finally, we compared our basin-averaged estimates against previous studies that attempted to reconstruct TWS. In particular, we used the estimates produced by [12] as a benchmark for our product. Such TWS estimates are a result of state-of-the-art hydrological models, i.e., the Noah land surface model with multi-physics options (Noah-MP: [66] and the Hydrological Modeling and Analysis Platform (HyMAP) river-routing scheme [67,68], simulating five major water storage components (i.e., groundwater, soil moisture, snow, canopy interception and surface water storage, which includes rivers and floodplains) at the global scale. Anthropogenic activities were neglected in their TWS estimate. The final product is the mean of a 12-member ensemble composed of model runs forced with different meteorological and precipitation datasets. More details on the ensemble generation can be found in [12].

GRACE-Based TWS versus Radar Altimetry
Radar altimetry data correlates well with TWS GC , with CC≥0.8 in most pixels (1st row in Figure 2). CC tends to increase as the spatial resolution degrades. At finer resolutions (0.5 • and 1 • ), high correlations (CC > 0.8) were mostly found in the central Amazon basin. Few points exhibited a weak, although positive, relationship (CC < 0.4) between TWS GC and RA data, and are located near the basin's boundary.

Evaluation of the Regression-Based TWS Estimates
Supported by the good correlation between TWS GC and water-level observations in the study area, we estimated TWS from RA data. TWS alt and TWS GC correlates well and deviations measured by the MAE indicate that errors are low (2 nd and 3 rd rows in Figure 2). Previous studies (e.g., [19]) showed Amazon basin-averaged TWS amplitudes of 15 cm, with values above 170cm over the center of the basin (see Figure S4). Few pixels showed MAE >10cm, and even less reached the maximum error (MAE~20cm). Most estimates had MAE lower than 10 cm, corresponding to <10% of local TWS GC amplitude.
Following the relationship noted between GRACE-based TWS and RA data, our estimates were closer to TWS at coarser resolutions. At 0.5 • and 1 • , the averaged CC were similar (~0.78), although MAE values were~0.4 cm lower at 1 • than those obtained at 0.5 • . At 2 and 3 • , CC values were mostly higher than 0.8 and MAE values lower than 5.6cm. Considering the straightforwardness of the proposed approach, the errors in most pixels were relatively low, indicating that the estimates were satisfactory. Evaluation metrics did not show much change between calibration and validation periods ( Table 2).

Interpolation Evaluation
TWS interp was satisfactory at all spatial resolutions. Comparisons against TWS GC show high CC (>0.8) and low MAE (<10cm) at most interpolated pixels (Figure 3). Correlation was lower for pixels near the basin boundary because of the lower number of VS available within a short radius, reducing the quality of the interpolation. TWS interp also shows increased MAE values towards the basin outlet, especially around the mainstream in the eastern part of region, where TWS alt also exhibited high MAE (>15 cm). Given the high RA data availability in that part of the basin, the low performance of both TWS alt and TWS interp could be related to the sea tide influence on the water level. According to [69], the tide effect can be noticed more than 1000 km upstream from the Amazon basin outlet. This means that the monthly mean water level derived from altimetry could be contaminated by the tide, reducing the correspondence with GRACE-based TWS.
The interpolation results in Figure 3 were obtained by setting the maximum search radius (MSR) varying from 3.5 • to 4.5 • , depending on the grid resolution, and a minimum of three pixels with VS within that radius. To test how the search radius affects the interpolation, we calculated the mean CC resulting from the kriging and NN interpolation as such radius increases. This analysis is necessary because the quality of the estimates decreases as MSR increases, but setting a small MSR means neglecting pixels near the basin boundary. Because the results from IDW and kriging methods were similar, Figure 4 shows the plots for kriging and nearest neighbor (NN) methods only. We noted that, using kriging (and IDW), correlation is above 0.8 for all MSR and about 100% of pixels were filled with MSR ≥4.5 • for 3 • spatial resolution and MSR ≥3 • for finer spatial resolutions. For the NN method, the gain of spatial coverage comes with a high cost to the quality: CC decreases faster than that using kriging as MSR increases, with CC falling below 0.8 for MSR ≥2.5 • and 10% of empty pixels (Figure 3). Based on that evaluation, we limited the search radius at 3.5 • for the 0.5 • grid, 4 • for the 1 • grid, and 4.5 • for 2 • and 3 • grids.

Basin-Average TWS
We aggregated all altimetry-based TWS estimates (TWS reg + TWS interp ) and averaged the final product (TWS alt ) over the Amazon River basin and basins of its main tributaries. Among the droughts and floods that have occurred in the region, we selected two extreme events, one during the dry season (MAM 2004) and one wet (DJF 2009), that covered practically the same area extent in the northern part of the Amazon Basin ( Figure 5). The time series for such 'extreme' area are shown in Figure 5, as well as selected metrics (MAE, CC and γ) between TWS GC and TWS alt for each basin and resolution.
TWS alt compared well with TWS GC over the study period (October 2002-July 2010), especially at 1 • spatial resolution, which was used to plot the time series in Figure 5. The amplitude, phase and shape of the estimated time series are very close to GRACE data. TWS alt for the 'Extreme' area is also similar to TWS GC , indicating that our estimates are not significantly affected by such extremes and can provide reliable TWS estimates under different scenarios (regular, droughts and floods).
Although the difference between the seasonal anomaly detected by TWS GC and TWS alt were relatively large for some pixels during MAM 2004, such deviations did not affect the basin average time series. For the wet period (DJF 2009), the seasonal anomalies are well represented both at the pixel scale as well as in the basin-average time series. As the spatial resolution degrades, we can see three different trends in MAE: an increase in the Negro and Solimões River basins and in the selected extreme event region; a decrease in the Madeira, Tapajós and Xingu River basins; and little change in the Amazon River basin ( Figure 6).
Except for Xingu at 0.5 and 1 • , and Negro at 3 • , MAE between TWS alt and TWS GC time series was below 6 cm. With a few exceptions, CC was above 0.95 for all resolutions and basins. The lowest CC was found for the Negro River basin at 3 • . At 3 • and compared to other basins, TWS alt for the Negro River basin was not as well represented: MAE >6cm and CC <0.9. This limitation is linked to the fact that, at such resolution, two among the few pixels in the Negro River basin have low CC between RA and TWS GC (see Figure 1). For finer resolutions, the number of well correlated virtual stations is higher, enhancing the quality of our estimates in that basin. The standard deviation ratio (γ) was the only metric with a similar (negative) trend among the analyzed cases, i.e., it decreased for all sub-basins as the spatial resolution degrades. Regardless of the basins and spatial resolution, γ > 1 show that the amplitude of TWS alt signal is greater than that from TWS GC . In most basins, γ decreases as spatial resolution degrades from 1 • to 2 • . At 3 • , γ converged to a range between 1.5 and 1.2 for the Amazon, Madeira, Negro and Tapajós basins, and below 1.1 for Extreme, Solimões and Xingu River basins.

Discussion
TWS alt is competitive compared to previous studies that have reconstructed or predicted TWS time series (e.g., [17,18,21]). As an attempt to reconstruct gridded TWS using the Global Land Data Assimilation System (GLDAS; [70]), [18] found that TWS estimates are poorly correlated with GRACE data, with CC <0.5 in about 50% of the basin. They also proposed a linear regression model to predict the basin-averaged TWS signals that satisfactorily captured the general variation patterns found in TWS GC (no metric was reported for that product). Based on a visual comparison between their plots and ours confirm the competitivity of our product. Not only our time series depicted the variation patterns of TWS GC , we show a good agreement for several sub-basins of the Amazon basin in terms of both time series and spatial distribution.
Water budget-based TWS time series by [20] was also derived for the entire Amazon basin. Performance metrics relative to less than a four-year period during which GRACE data was available are not reported in his paper. The largest difference reported between TWS GC and its estimates occurred in 2003, ranging from five to ten cm. Using seven time series, as shown in Figure 6, each~eight-year long, we quantified the deviations to the basin-averaged GRACE-based TWS. Only seven months (3%) had deviations larger than 15 cm, whereas 300 months (~47%) had deviations smaller than 5 cm ( Figure S5).
Reference [41] found correlations ranging from -0.31 to 0.77 between their reanalysis-based TWS and GRACE estimates for a number of basins across North America, Europe, Asia and Australia. Using water level data at 58 in situ stations, [19] reconstructed TWS signal over the Amazon using in situ water-level data and obtained correlations of about 0.9 between their time series and that from the basin-averaged GRACE data. [21] reconstructed global 0.5 • TWS time series over~30 years and obtained correlations (0.94-0.97) and MAE (2.3 cm-5.1cm) similar to ours, over the Amazon basin and sub-basins. In terms of standard deviation ratio, their product overperformed ours: while we found 1.02 < γ <1.17, that metric ranged from 1.02-1.06 with their dataset.

Conclusions
Here, we presented a simple and efficient method to estimate monthly gridded TWS over the Amazon Basin, at different spatial resolutions, fully based on remote sensing data. This study offers a straightforward option for gap-filling GRACE-based TWS data and reconstruct past TWS over regions where surface water accounts for an important part of TWS signal, such as the Amazon basin [6]. Unlike previous approaches that relied on ancillary information (e.g., hydrological model outputs, meteorological forcing and reanalyzes), the proposed method only requires radar altimetry time series. It consists of using regression coefficients between GRACE-based TWS and RA-based water levels, to predict TWS solely from radar altimetry data. Radar altimetry data is available from 1992 to present and can be acquired near real time, allowing a low-latency TWS estimate, as opposed to GRACE's several-month latency.
We showed that altimetry-based TWS compares well with GRACE-based TWS, as indicated by the performance metrics, and that selected extreme events can be satisfactorily represented. Although altimetric data is not available with enough density to provide TWS estimates all over the domain at finer resolutions (i.e., 0.5 • and 1 • ), we were able to reproduce spatially continuous gridded products, at such grid sizes, by means of interpolation of regression-based TWS estimates. Results also demonstrated that TWS time series can be reconstructed with high accuracy at coarser spatial resolutions over the Amazon basin with as few as 41 virtual stations.
Despite the promising results, future analyses are needed in order to determine the applicability of the method over other regions and, ultimately, at the global scale. The main limitations that need further investigation could be found in regions with low RA data availability or low influence of surface water on TWS, in particular snow-dominated regions. Although the radar altimetry dataset adopted in this study provides a large number of virtual stations over the Amazon basin, data gaps limited its use to the 2002-2010 period. As a consequence, the recent development of tools to automatically extract multi-mission radar altimetry data (e.g., [43]) could contribute to increase data availability. Nonetheless, our approach has a great potential to provide high-temporal and -spatial estimates of TWS in a near future, with the launch of the Surface Water and Ocean Topography (SWOT) mission, scheduled for 2021.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/21/2487/s1, Figure S1. Representative water level (WL) anomaly time series containing gaps and result of the adopted gap filling method; Figure S2. Number of altimetry data per month; Figure S3. Pixel count through time for each spatial resolution; Figure S4. Spatial distribution of GRACE-based terrestrial water stoarge amplitude over the Amazon basin. Amplitudes result from the difference between average climatological maximum and minimum values estimated for 2002-2017; Figure S5. Histogram of differences between the time series of GRACE-based (TWSGC) and altimetry-based (TWSalt) terrestrial water storage (TWSGC-TWSalt,) at 1 • spatial resolution; Table  S1