Interactive comment on “ Past terrestrial water storage ( 1980 – 2008 ) in the Amazon Basin reconstructed from GRACE and in situ river gauging data ” by M .

Abstract. Terrestrial water storage (TWS) composed of surface waters, soil moisture, groundwater and snow where appropriate, is a key element of global and continental water cycle. Since 2002, the Gravity Recovery and Climate Experiment (GRACE) space gravimetry mission provides a new tool to measure large-scale TWS variations. However, for the past few decades, direct estimate of TWS variability is accessible from hydrological modeling only. Here we propose a novel approach that combines GRACE-based TWS spatial patterns with multi-decadal-long in situ river level records, to reconstruct past 2-D TWS over a river basin. Results are presented for the Amazon Basin for the period 1980–2008, focusing on the interannual time scale. Results are compared with past TWS estimated by the global hydrological model ISBA-TRIP. Correlations between reconstructed past interannual TWS variability and known climate forcing modes over the region (e.g., El Nino-Southern Oscillation and Pacific Decadal Oscillation) are also estimated. This method offers new perspective for improving our knowledge of past interannual TWS in world river basins where natural climate variability (as opposed to direct anthropogenic forcing) drives TWS variations.


Introduction
Terrestrial water storage (hereafter noted as TWS) is an important component of the global and continental water cycle.TWS refers to the total amount of water integrated over depth, stored in a catchment area.It is comprised of surface waters, soil moisture and underground waters.In some Correspondence to: M. Becker (melanie.becker@legos.obs-mip.fr)regions, TWS also includes snow.Variation of TWS with time t is linked to accumulated precipitation (P ), evapotranspiration (ET), surface and subsurface runoff (R) within a given area or basin, through the water balance, written as: Quantification of TWS variability and change is difficult because limited ground water level and soil moisture observations are available, and often are simply inadequate or inexistent (e.g., Rodell and Famiglietti, 1999;Shiklomanov et al., 2002;Alsdorf et al., 2007;Liu and Yang, 2010;Liu et al., 2009) at basin or smaller scales.Global and regional hydrological models developed for water resources assessment and climate research purposes provide an alternative to missing in situ measurements (e.g., Döll et al., 2003;Milly andShmakin, 2002, Rodell et al., 2004).Some of these models compute the water and energy balance at the Earth' surface, yielding -among other parameters, temporal variations of the total water storage in response to prescribed forcing (solar radiation and precipitation) and variations of near-surface atmospheric conditions.Since 2002, the Gravity Recovery and Climate Experiment (GRACE) space mission provides gridded time series of TWS at monthly or sub-monthly interval, with a resolution of ∼300-400 km (Tapley et al., 2004;Wahr et al., 2004) and a precision of ∼2 cm in equivalent water height (EWH) (in the following EWH refers to spatially distributed water storage while TWS refers to water storage averaged over a given area).GRACE provides a highly valuable new data set that allows studying water storage change over large river basins worldwide, complementary to precipitation, in situ river level and discharge data.However the GRACE lifetime is still short and does not allow studying past decade variability of Published by Copernicus Publications on behalf of the European Geosciences Union.
TWS.Only hydrological model outputs can be used for that purpose.
We have developed a novel method to reconstruct 2dimensional (i.e., gridded) TWS over the past decades (since 1980).The method is similar to that classically used to reconstruct past atmospheric or oceanic fields such as marine sea level pressure (Kaplan et al., 2000), sea surface temperature (Smith and Reynolds, 2003) and sea level (Chambers et al., 2002;Church et al., 2004;Llovel et al., 2009;Calafat and Gomis, 2009).The method developed in the present study combines spatial information on TWS from GRACE (over 2003-2008) with multi-decade-long (over 1980-2008) but sparse river level time series based on in situ gauges.The past reconstructed TWS grids cover the 1980-2008 period.The method is applied to the Amazon Basin and is focussed on interannual variability.
Compared to classical reconstructions applied to atmospheric or oceanic data sets which combine grids and sparse records of the same physical quantity (e.g., surface atmospheric pressure or sea level), here this is not the case as we combine gridded TWS from GRACE with in situ river level records.River level is a component of the total TWS, related in a nonlinear way to inundation extent and thus surface water volumes.However, in the Amazon Basin, previous studies (e.g., Xavier et al., 2010;Vaz de Almeida, 2009) have shown that at seasonal and interannual time scales, river water level fluctuations can locally be correlated to TWS (as observed by GRACE).Such a correlation suggests that, at these time scales, TWS (including underground waters) and surface waters co-vary in a similar way.In the present study, we take advantage of this correlation at the interannual time scale and compute a scaling factor between river level and GRACE-based TWS over the GRACE time span (since 2002).This allows us to construct virtual multi-decade long TWS time series at the gauging sites for further combination of 2-D TWS grids from GRACE (of limited time duration) with sparse but long virtual TWS time series (based on the re-scaled river level time series).The final products are gridded (i.e., 2-D) time series of past TWS.

In situ river level data
In this study, we use water level data from the in-situ gauging stations instead of river discharge data because it is one of the TWS components, unlike discharge.Since direct measurements of discharge in river channels can be time-consuming and costly, flow is commonly estimated indirectly by means of a curve relating stage (river level) to discharge (Clarke, 1999).Hence, uncertainties of stage measurements and rating curve method increase the final uncertainty.Using river level data is thus more straightforward (but note that our reconstruction method would also work with discharge data).We considered 58 in situ gauge sites with almost continuous water level time series over 1980-2008.The data are available from the Brazilian water agency ANA -Agência Nacional de Águasnetwork (www.ana.gov.br).The location of the 58 sites is shown in Fig. 1.The river level time series are located on the Amazon River and on some of its tributaries.They cover the period from January 1980 to December 2008 and are given at monthly interval.Only time series with gaps smaller than 2 consecutive years are considered (see Table 1).These gaps are then filled by linear interpolation.This dataset is subjected to outlier analysis in order to identify and remove extreme values that may lead to an incorrect interpretation of the data.In the present study, outliers are detected using the Rosner's test (Rosner, 1975).As we focus here on interannual to multidecadal time scales, we removed the seasonal cycles from in-situ river level data.The seasonal cycles in these data was removed by fitting sinusoids with periods of 12 and 6 months, before filling the gaps.

GRACE TWS data
The GRACE space mission, jointly developed by NASA and DLR (German Space Agency), was launched on 17 March 2002.It utilizes a state-of-the-art technique to measure temporal variations of the Earth's gravity field by tracking through a K-band ranging (KBR) system, the intersatellite range and range rate between two coplanar, low altitude satellites (GRACE A and B) (Tapley et al., 2004).In addition, each satellite is equipped with a SuperSTAR Accelerometer, GPS receiver/antenna, Star Cameras, and Hydrol.Earth Syst.Sci., 15, 533-546, 2011 www.hydrol-earth-syst-sci.net/15/533/2011/  pur, 2007;Flechtner, 2007).Time variable GRACE global gravity solutions are provided by three GRACE data processing centers of the Science Data System (SDS): Center for Space Research (CSR) at the University of Texas at Austin, the Geoforschungszentrum (GFZ) in Potsdam, and the NASA Jet Propulsion Laboratory (JPL).The GRACE solutions are distributed by NASA PODAAC (http://podaac.jpl.nasa.gov/grace/).Other groups external to SDS also provide GRACE solutions, e.g., the Goddard Space Flight Center (NASA; Rowlands et al., 2002), the Delft Institute of Earth Observation and Space Systems (DEOS; Klees et al., 2008), the Groupe de Recherche de Geodesie Spatiale (GRGS; Lemoine et al., 2007), among others.The GRACE monthly (or sub monthly) solutions are generally expressed in the form of spherical harmonic coefficients of the geoid height, up to some maximum degree (typically between 60 and 100, corresponding to wavelengths of ∼400 to 700 km).Gridded time series are also available, in general expressed in terms of EWH.Since the beginning of the mission, different GRACE solutions releases have been made available by the SDS groups, with improved quality from release to release.Here we use GRACE products (release 2) computed by the Groupe de Recherche de Geodesie Spatiale -GRGS (Bruisma et al., 2009).These are monthly EWH solutions provided as 1 • × 1 • global grids from January 2003 through December 2008.As we focus here on interannual to multidecadal time scales, we removed the seasonal cycle from the gridded GRACE EWH by fitting sinusoids with periods of 12 and 6 months at each grid mesh.
We selected GRACE-based EWH data over the Amazon Basin.For small basins, it is necessary to correct for a leakage factor (due to gravitational signal from outside the considered basin; Chambers, 2006).However over the Amazon Basin, the leakage correction is ∼5% of the seasonal signal (Chen et al., 2007;Ramillien et al., 2008;Xavier et al., 2010).To confirm this result, we estimated the leakage error on the Amazon Basin following Klees et al. (2007) and Longuevergne et al. (2009).Outputs from the Global Land Data Assimilation System (GLDAS-NOAH) (Rodell et al., 2004) has been considered as an a-priori information to compute the leakage error over the period [2003][2004][2005][2006][2007][2008].We obtain a leakage error of about 5% of the seasonal cycle, as in previous studies, and around 3% of the interannual signal.As we focus here on the interannual variability, we conclude that the leakage error is negligible.
Figure 2 shows the first three leading modes of the Empirical Orthogonal Function (EOF) decomposition (Preisendorfer, 1988) of gridded TWS based on GRGS GRACE data over 2003-2008.The EOF mode 1 is dominated by a strong positive signal affecting the Negro subbasin and a small area in the southern part of the basin.The EOF mode 2 suggests that the Amazon Basin is divided into two main hydrological zones (west and east).The EOF mode 3 shows the mid-2005 drought that affected the centre and the western part of the basin (Chen et al., 2009).These EOF results are very similar to those obtained by Xavier et al. (2010), who made an EOF decomposition of TWS from the CSR GRACE solutions filtered by Chambers (2006).

Relationship between GRACE-based TWS and in situ river levels
The × 1 • GRACE pixel (in order to not weight over pixels that contain more than one in situ station).The actual GRACE resolution is closer to ∼300 km (e.g., Schmidt et al., 2008) but tests made by averaging the in situ data in grid mesh of 3 • × 3 • did not show significant difference in the reconstruction results.Computing the regression slope between the two data sets is equivalent to computing the ratio (called scaling factor) between the standard deviations of both river level and GRACE EWH time series (after removing their mean).The correlations and the scaling factors for each of 55 in-situ river level stations are gathered in Table 1 (see Fig. 1 for their location).Figure 3 shows in situ water level and GRACE EWH for a subset of 16 sites, chosen for their wide distribution across the basin (note that when 3 sites are located inside a single 1 • 1 • GRACE mesh, only mean location and mean river level time series are shown in Fig. 3; in these plots, river level time series are expressed in EWH using the computed scaling factor over span [2003][2004][2005][2006][2007][2008].From Table 1 and Fig. 3, we observe significant correlation on interannual time scale between the two data sets (river level and total TWS) at a given location.We also checked that the time series of the EOF decomposition of both GRACE and gauges are consistent over the GRACE period (not shown here).
Using the scaling factor computed over the validation period 2003-2008, TWS virtual records were then reconstructed backward in time (back to 1980) from the river level time series.In the reconstruction we only use in situ river levels that do verify a correlation higher than 0.7 with the closest GRACE data point.Only 23 sites verify the correlation among a total of 55.This may result from a different hydrodynamic behaviour between total water storage and surface waters (possibly because of the presence of seasonal floodplains; Alsdorf et al., 2000; or confluence with a tributary).In situ virtual TWS records used to reconstruct the past 2-D TWS fields will have in common that they are dominated by their surface water component variability.For the Negro subbasin Frappart et al. (2008) confirmed that the surface water component is not negligible in the interannual TWS, which is almost equally partitioned between surface water and the combination of soil moisture and groundwater.We cannot exclude that in other regions of the Amazon Basin, the relationship does not hold.At first sight, this could bias the reconstruction process.Actually this is not the case, as we will see below (see method in Sect.4) since the functions (EOFs) used to interpolate the virtual records in the reconstruction process are statistical modes of the total TWS variability that intrinsically take into account the co-variability of the different layers of the soil at different places (the method optimally interpolates virtual TWS records that are all dominated by their surface waters).In other (not correlated) areas, reconstructed spatial patterns will be based on the statistical information contained in GRACE TWS grids (provided that these patterns are stationary; see discussion below).

Method
The method used to reconstruct past (over 1980-2008) TWS in 2-dimension over the Amazon Basin is based on the reduced optimal interpolation described by Kaplan et al. (2000) and used by Church et al. (2004) to reconstruct past sea level.This method has 2 steps.In the first step an EOF decomposition (Preisendorfer, 1988;Toumazou and Cretaux, 2001) of the GRACE-based TWS grids is done.This decomposition allows to separate the GRACE signal (here represented by a matrix H, with m lines for each spatial point and n columns for each date) into spatial modes (EOFs) and their related temporal amplitude as follow:, H(x,y,t) = U (x,y) α(t). (2) In this equation U (x,y) stands for the spatial modes and α(t) for their amplitude over the GRACE period.Assuming that the spatial modes U (x,y) are stationary in time, we deduce that the reconstructed TWS field of the Amazon Basin over the long period 1980-2008 (called here H R (x,y,t)) has an EOF as follow: where α R (t) represents the new amplitudes of the EOFs over 1980-2008.
The second step consists of computing the new amplitudes α R (t) over the whole period 1980-2008 thanks to the in situ virtual TWS records.It is done through a least squares optimal procedure that minimizes the difference between the reconstructed field and the in situ virtual TWS records at the in situ gauge locations.
In the first step, the EOF modes and amplitudes of the GRACE data set matrix H are computed through a singular value decomposition approach, such that: where U (x,y) still stands for the EOF spatial modes, S is a diagonal matrix containing the singular values of H and V represents the temporal eigen modes.At this stage the amplitude α(t) of the EOF modes can be simply written as α(t) = SV T .Conceptually, each EOF k (k-th column of U (x,y) multiplied by the k-th line of α(t): U k (x,y)α k (t)) is a spatio-temporal pattern of TWS variability that accounts for a percentage of the total variance of the TWS signal H.
As stated in Eq. ( 2), the computation of the EOFs is purely statistical and has no information on the layers of the soil involved in the local variance of the EOFs.This information is actually randomly distributed in each EOFs so that each one of them carries variability that comes from all the layers.Obviously if at a location (x s , y s ) the total signal H(x s , y s ) is dominated by its surface water component, then for each EOF k, the signal U k (x s , y s ) will also be dominated by the surface water.But at the same time, for other points (x j , y j ) where this is not the case, then U k (x j , y j ) will carry variability from other layers.The U k (x,y) modes intrinsically take into account the co-variability of the different layers of the soil at different locations.
The low-order EOFs (eigenvectors of the largest singular values) explain most of the variance and contain the largest spatial scales of the signal.The higher-order EOFs contain smaller spatial scale patterns and are increasingly affected by noise.Besides, their amplitude is decreasingly well resolved by the least squares procedure because the sparseness of the set of in situ gauges does not allow to resolve too small scale patterns.Consequently, to be efficient, the TWS reconstruction over the Amazon Basin uses only a subset of the M lowest-order EOFs (the best fit between maximum variance explained and minimum noise perturbation led us to choose M = 3, which accounts for 79% of the total variance of the GRACE data).
Consequently; the data matrix H can be written as where α(t) = S M V T M is the matrix of the amplitude of the M lowest EOFs.Following Kaplan et al. (2000), in the second step, we compute, at each time step over the time span of the in situ records, the amplitudes α R (t) by minimizing the cost function: In Eq. ( 6) H 0 is the in situ observed TWS, P is a projection matrix equal to 1 when and where in situ records are available and 0 otherwise and is a diagonal matrix of the M largest eigenvalues of the covariance matrix.R is the error covariance matrix accounting for the data error covariance matrix (instrumental error) and the error due to the truncation of the set of EOFs to only the first M EOFs.The second term on the right hand side of the function is a constraint on the EOF spectrum of the solution.It prevents the least squares procedure to be contaminated by high-frequency noise (it filters out non significant solutions that display too much variance at grid points without nearby observations).The least squares procedure is then applied to the virtual in situ TWS records (H 0 ) deduced from the in situ water level records.It provides the reconstructed amplitude α R of the EOFs.Since virtual TWS records are all relative to their own local datum that are not cross-referenced over the basin, this solution may be polluted by spatial variability of the in situ TWS reference surface not necessarily consistent with the GRACE TWS reference surface.To cope with this problem we solved Eq. ( 3) for changes in TWS between adjacent steps following Church et al. (2004).Once changes in amplitude have been obtained at each time step, the amplitudes themselves have been recovered integrating backward in time.The integration constants are chosen to equal the reconstructed EOF amplitudes mean to the GRACE EOF amplitude mean over the GRACE measurement period, ensuring consistency between both sets of EOFs.Finally the reconstructed field of TWS over the Amazonian Basin is obtained by multiplying the first three EOFs with their reconstructed amplitude:

Stationarity of the spatial patterns
Since the spatial structure of EOFs is sensitive to noise in the observational dataset, we use the longest GRACE dataset available.The GRACE record used in this study spans only 6 years between January 2003 and December 2008 (this is somewhat less but not significantly different from the 9 years of TOPEX/Poseidon altimetry data used by Church et al., 2004, to reconstruct past 50 years sea level).However, as full coverage of the Amazon Basin's TWS is only available for the GRACE period, there is no other way to determine TWS spatial patterns.As an alternative, EOFs could have been computed from long hydrological model outputs.But while hydrological models agree rather well in terms of TWS spatial average, they much differ when looking at the spatial patterns.
One possible drawback of using GRACE to compute spatial EOFs is linked to the question of (non) stationarity in time of the spatial patterns.Here, we assume that the EOF spatial patterns of the GRACE TWS capture most of the interannual/decadal variability.To test this stationary assumption (inherent to the method used here), we consider a dataset composed of 23 in situ gauges (black dots in Fig. 4, see below), spanning over the period 1980-2008.We computed the river level EOFs over three interval subsets (2003-2008, 1990-2008 and 1980-2008).Corresponding first 3 leading modes (>80% of the total variance) are shown in Fig. 4. Very little difference is noticed between the three cases.It is striking to see how the EOFs are similar despite the different time span over which they have been computed.
In addition, we computed the EOFs of the monthly 1 • × 1 • gridded precipitation produced at the Global Precipitation Climatology Center (GPCC) (Rudolf, 1995) over the same three interval subsets (2003-2008, 1990-2008 and 1980-2008).As we focus here on interannual time scale, we removed the seasonal cycle from the gridded GPCC precipitation data.The first 3 leading modes for each case are shown in the Fig. 5. Again, very little difference is noticed between the three cases, suggesting that the shortest time span captures well most of the interannual variability.These results suggest that the main spatial structures are present in the in situ data and in precipitation data over the three periods considered (2003-2008, 1990-2008 and 1980-2008).We conclude that surface water patterns (hence TWS fields, because of the reported correlation) and precipitation patterns are quasi stationary with time (at least as far as the early 1980s).Thus the basic assumption of our reconstruction method (stationarity of the spatial EOFs) holds.

Reconstructed TWSR averaged over the Amazon Basin
Figure 6 shows past TWSR  averaged over the Amazon Basin (a 5-month smoothing was applied to the time series; TWSR holds for reconstructed TWS).For comparison, we also show the reconstructed TWSR (with same smoothing) over 1990-2008.For the latter case, a larger number (36) of in situ stations well correlated (r > 0.7) with GRACE TWS were available.Hence a larger number of virtual TWS time series could be used.Comparing the two reconstructions shows very little difference, suggesting that the number of virtual stations is not critical, provided that they are well distributed across the basin (see Fig. 1).In Fig. 6 is also shown GRACE-based mean TWS, which superimposes well with the reconstructed curves.We compared mean TWSR over 1980-2008 with mean TWS computed by the ISBA-TRIP hydrological model.ISBA is a hydrological model that uses the force-restore method to calculate the time variation of the surface energy and water budgets (Noilhan and Planton, 1989).The soil hydrology is represented by three layers: a thin surface layer (1 cm) included in the rooting layer and a third layer to distinguish between the rooting depth and the total We filtered out the high frequencies using a simple 5-month running mean.r is the correlation coefficient (with a significance level higher than 99%).Fig. 7. Basin-averaged of the TWSR comparison with ISBA-TRIP.The basin-averaged of the TWSR is the black line and its error bars are in grey.The error in TWSR computed here is the sum of the error due to the least squares method and the error of the in situ records.This signal accounts for 60% of the total reconstructed signal variance.ISBA-TRIP is in red line.The El Niño events are represented by arrows.We filtered out the high frequencies (>1 yr −1 ) using a simple 12-month running mean.r is the correlation coefficient (with a significance level higher than 99%).soil depth.Several surface water storage compartments as dams, lakes, or groundwater storage are not simulated in the ISBA-TRIP model.The soil water content varies with surface infiltration, soil evaporation, plant transpiration and deep drainage.The infiltration rate is computed as the difference between the through-fall rate and the surface runoff.The through-fall rate is the sum of rainfall not intercepted by the canopy, dripping from the interception reservoir and snowmelt from the snow pack.ISBA also uses a comprehensive parameterization of sub-grid hydrology to account for the heterogeneity of precipitation, topography and vegetation within each grid cell (Decharme and Douville, 2007).TRIP -Total Runoff Integrating Pathways-was developed by Oki and Sud (1998).It is a simple runoff routine model used to convert the daily runoff simulated by 1 • by 1 • resolution.The ISBA-TRIP version used in this study is driven by prescribed atmospheric forcing using monthly precipitation data from the Global Precipitation Climatology Center (GPCC) Full Data Product V4 (Alkama et al., 2010;Decharme et al., 2010).The temporal resolution of this data set in 1 month.
A fair comparison with GRACE observations requires that ISBA fields be spatially filtered in a similar way.To do this, ISBA-TRIP TWS gridded fields were expanded in spherical harmonic (SH) functions.SH coefficients were truncated at degree and order 50 (corresponding to a spatial resolution of ∼400 km).Then, new gridded ISBA-TRIP TWS were computed with the ISBA-TRIP SH truncated at degree 50.TWS based on ISBA-TRIP was averaged over the Amazon Basin, and as for GRACE data.The seasonal signal was removed data were smoothed with a 12-month running mean filter.ISBA-TRIP-based mean TWS for the 1980-2008 period is shown in Fig. 7.In Fig. 7, mean TWSR (as shown in Fig. 6) is superimposed.Looking at Fig. 7, we note that the two time series are well correlated, both in amplitude and timing.The correlation between the two time series amounts to 0.9 (significant at the 99% level), giving confidence in the mean TWSR, since the ISBA-TRIP simulation is based on a totally independent approach.Conversely, this mutual agreement provides another validation of the ISBA-TRIP model (see Alkama et al., 2010).Non account of ground water storage in ISBA-TRIP could explain the small differences that remain between the two curves.We observe TWS minima during ENSO events (indicated by arrows in Fig. 7), as expected from previously reported correlation between ENSO and rainfall and in situ discharge data over the Amazon Basin (Marengo, 2004;Molinier et al., 2009;Ronchail et al., 2002).The lack of trend in both curves (reconstruction and ISBA-TRIP model) suggests no net gain or loss in total water storage over the Amazonian Basin between 1980 and2008.In Fig. 7, the grey zone represents uncertainty in TWSR.This uncertainty is based on the sum of errors due to the least-squares inversion (as presented in Sect.4) and errors of the in situ records.This latter error is estimated from a bootstrap method for standard errors of the in situ water levels for each month (significant at the 95% level).Actually, the reconstruction method includes additional uncertainties.For example, the scaling factor between in situ river level and GRACE-based TWS may introduce some uncertainty.Indeed, river level may be more closely tied to the surface wetness condition than to groundwater (or TWS).Another source of error possibly be due to precipitation events in upstream areas, may affect downstream water levels with some time delay due to water transport.GRACE measurements have also their own uncertainty.As a matter of fact, GRACEbased TWS is not a point measurement (as it is the case for river level) but represents an average over a much larger region (of about 300 km size).All these uncertainties will affect the precision of the reconstruction but they are very difficult to estimate.They have been neglected here.Thus the grey zone likely underestimates of the actual uncertainty.

Spatial patterns of TWSR
We have analysed the spatial patterns of the 2-D TWSR fields over 1980-2008 and 1990-2008, through an EOF decomposition approach.Results are shown in Fig. 8.The spatially constant mode (EOF0, around 60% of the total variance for each TWSR) which represents the interannual variability of the mean TWSR over the Amazon Basin from 1980-2008 is that shown in Figs. 6 and 7.
Figure 8 shows the first three EOF modes of the TWSR fields over the two time spans (1990-2008 and 1980-2008).Recall that the 1990-2008 reconstruction is based on a larger set of in situ station than the 1980-2008 one (36 versus 23).Looking at the spatial pattern maps and at the temporal curves, we note quite good agreement between the two reconstructions (1990-2008 and 1980-2008), as previously noticed for the mean TWSR.On EOF1 (26% and 25% of the total variance) temporal curve, we have superimposed the Southern Oscillation Index (SOI; mean sea-level pressure difference between Tahiti and Darwin), a proxy of ENSO.Correlation between SOI and TWS mode 1 is high (r = 0.7 with a SL >99%).In the north-eastern region, including the Negro sub-basin, and in the Madeira basin, where EOF1 exhibits a strong positive anomaly in the spatial map, minima in local TWS correspond to negative SOI (warm phase of ENSO).In contrast, a negative correlation can be observed locally between SOI and TWS mode 1 in the Tapajos and the Xingu sub-basins.Studies have shown that precipitation over the northeastern part of the Amazon Basin is largely controlled by ENSO (Zeng, 1999;Marengo, 1992;Espinoza et al., 2009).Thus, a positive correlation between mean TWS and SOI is expected, as our result indeed shows.
Studies provided diagnostic evidence of ENSO low frequency modulation (McCabe and Dettinger, 1999;Gutzler et al., 2002), in particular by the Pacific Decadal Oscillation (PDO) (Mantua et al., 1997).The PDO index is defined as the leading principal component of the North Pacific monthly sea surface temperature variability.The PDO is a long-lived El Niño-like pattern of Pacific climate variability.In the cool (warm) phases of PDO, the central and northwest Pacific sea surface temperature (SST) is warm (cold) and the SST at the coast of the North America is of cold (warm) SST.In Fig. 8, is displayed the EOF2 temporal curve of TWS with the PDO index superimposed.Correlation between PDO and TWS EOF2 is significant (r = −0.5 with a SL >99%).In the whole eastern Amazon Basin, where the spatial EOF2 mode shows a negative pattern, local TWSR temporal variability is negatively correlated to the PDO index, while in the Negro and the Solimões sub-basins a positive correlation is observed.
The spatial pattern and temporal curve of the EOF3 in Fig. 8 reflect the recurrent droughts that affect the centre of the Amazon Basin, in particular the main river.These are indicated by arrows on the temporal curve.The well publicized 2005 drought is clearly visible.It has been reported to be one of the most intense droughts of the past 100 years in the Amazon Basin (Marengo et al., 2008;Zeng et al., 2008;Chen et al., 2009), but the TWSR mode 3 seems to indicate that the 1980 and 1998 droughts had even stronger impact in terms of basin total water storage depletion.Drought events in 1983and 1998occurred during El Niño years, unlike in 1980, 1992, 1995and 2005.

Conclusions
The present study has established for the first time direct, observation-based, estimate of TWS spatiotemporal variability over the Amazon Basin for the past ∼3 decades.It is only based on TWS observations from GRACE and (re-scaled) in situ water levels.TWSR appears able to depict interannual to multidecadal variability in water storage and associated spatial patterns.Such an approach is complementary to traditional methods, e.g., moisture convergence method (Zeng, 1999;Masuda et al., 2001) or PER method, where the key input variables are observed precipitation P and runoff R and estimated evaporation ET.In these methods, the authors apply the basin water budget equation to diagnose the long-term variability of total TWS (Zeng et al., 2008).The observed timing and regional distribution of the TWSR over the Amazonian Basin during ∼30 years shows no long-term trend and confirms the dominant influence of ENSO and PDO.In addition, recurrent drought events affecting the centre of the basin are also well reproduced.The approach developed in this study offers interesting perspective for improving our knowledge of past TWS in many river basins over the world where climate variability is the main driver of TWS change.However, it will be less easily applicable in river basins which have been strongly affected by anthropogenic forcing over the GRACE period, such as in northern India.Indeed, in these basins, GRACE measurement may be impacted by the recent ground water mining used for crop irrigation (Rodell et al., 2009;Tiwari et al., 2009).Thus GRACE EOFs will exhibit the spatial signature of this recent anthropogenic forcing and will be inappropriate to reconstruct past TWS variability in such regions.Some caution should then be highlighted in order that imprint of such events in GRACE data be not erroneously extrapolated backward in time.

Fig. 1 .
Fig. 1.Amazon River watershed with its main sub-basins.Location of the in situ river stages is indicated by dots.The red dots are the stations used in the reconstruction over 1980-2008 and correspond to a correlation coefficient ≥0.7 with GRACE TWS and in-situ level data over 2003-2008.The stations with a correlation coefficient in the range (0.5-0.7) are in green dots and in grey dots for a correlation coefficient <0.5.

Fig. 2 .
Fig. 2. EOF analysis of TWS from GRACE data over 2003-2008.Spatial patterns of the EOF decomposition of TWS.The three modes are arranged from top to bottom.The principal components are given in the right.

Fig. 3 .
Fig. 3. Scaled EWH for in situ station.The figure shows the scaling result for some in situ stations (see red dots on Fig. 1 for locations).The TWS from GRACE over 2003-2008 period is plotted in black lines and the in situ river levels scaled in EWH over the same period are in blue lines.r is the correlation coefficient and SL its significance level.

Fig. 4 .
Fig. 4. Spatial patterns of the third EOFs computed from in situ gauges.The black dots correspond to the 23 in situ gauges (the same 23 used in the reconstruction) used for the EOF analysis over 2003-2008 (left panel), 1990-2008 (middle panel) and over 1980-2008 (right panel).

Fig. 6 .
Fig. 6.Basin-averaged of the TWSR.The basin-averaged of the TWSR over 1980-2008 is the black line.The basin-averaged of the TWSR over 1990-2008 is the red line.The TWS GRACE data is superimposed in dot line over 2003-2008.We filtered out the high frequencies using a simple 5-month running mean.r is the correlation coefficient (with a significance level higher than 99%).

Fig. 8 .
Fig. 8. EOF analysis of TWSR.The left panel shows the EOF analysis for the TWSR over 1990-2008.The locations of the 36 in situ stations used for the reconstruction are indicated by black dots.The right panel shows the EOF analysis for the TWSR over 1980-2008.The locations of the 23 in situ stations used for the reconstruction are indicated by black dots.The middle panel shows the EOFs' time series computed on the TWSR over 1990-2008 in red line and in black line over 1980-2008, and superimposed (a)the scaled SOI index (dash line)and the El Niño events (arrows), (b) the inverse scaled PDO index (dash line) and (c) the drought events (arrows).We filter out the high frequencies (>1 yr −1 ) using a simple 12-month running mean.

Table 1 .
Characteristics of the in-situ river level stations: location (latitude, longitude), name of the sub-basin, percentage of gaps in the record over 1980-2008 and maximum of consecutive months missing in the record (in the brackets, "-" means that no month is missing).Laser Retro Reflector.The GRACE Science Data System uses measured inter-satellite range and range rate data, along with ancillary data, to estimate monthly (or sometimes submonthly) time series of global Earth's gravity fields (Bettad- River, and in the Madeira and Negro sub basins, in particular at interannual time scale.We computed the regression function between in situ river level and GRACE-based EWH time series (annual cycle removed), averaging all in situ river level records available in each 1 • Xavier et al. (2010)alysis consists of expressing interannual river level data in terms of local EWH time series, taking advantage of the local correlation existing between in situ river level and GRACE-based TWS measurement at the river gauge location.Xavier et al. (2010)investigated this correlation and showed it is quite significant along the Amazon Hydrol.Earth Syst.Sci.,15,[533][534][535][536][537][538][539][540][541][542][543][544][545][546]2011www.hydrol-earth-syst-sci.net/15/533/2011/