Impact of Pacific and Atlantic sea surface temperatures on interannual and decadal variations of GRACE land water storage in tropical South America

We analyze 10 years of Gravity Recovery and Climate Experiment (GRACE) terrestrial water storage anomalies (TWSAs) over tropical South America along with seven climate indices linked to equatorial Pacific and tropical Atlantic oceans sea surface temperatures (SSTs) using a multichannel singular spectrum analysis and lagged cross correlations. We focus on the interannual, nonlinear modes of covariability between TWSAs and SSTs. By comparing the relative distributions of the leading modes, we identify teleconnections between TWSAs, Pacific and Atlantic SSTs at different time periods. Thus, the northern and northeastern regions of tropical South America are mainly influenced by Pacific SSTs, while the central and western Amazon regions are more influenced by Atlantic SSTs. The former regions are more sensitive to central Pacific SSTs than to eastern Pacific SSTs. A quasi‐biennial mode explains the largest part (27%) of the residual, interannual cross covariance and is found both in the El Niño–Southern Oscillation and in the Atlantic meridional mode. A trend‐like mode explains the second largest part (24%) of the residual cross covariance and may be caused by the following: (1) the decadal variability in the North Pacific climate, as expressed by the negative trend in the Pacific decadal oscillation and by increased water storage in northern and northeastern South America, (2) the melting of Andean glaciers in Peru and Bolivia due to man‐induced increase in land surface temperatures, and (3) the land use/cover changes after deforestation leading to increased runoff and groundwater recharge, expressed by increased water storage in southern Amazon regions.


Introduction
[2] The Amazon River basin and neighboring basins are tightly connected to Earth's climate through strong interactions with the atmosphere and the ocean. For example, the Amazon basin is a major driver of atmospheric CO 2 variations [Bousquet et al., 2000] and provides almost one fifth of global continental fresh water discharge [Dai and Trenberth, 2002]. The hydroclimatology of tropical South America is characterized by a strong annual cycle forced by the South American summer monsoon. Annual water storage changes in equatorial South America east of 70°W are mostly correlated to moisture transport from the Atlantic Ocean, while the western equatorial South America is correlated to moisture advection from the Pacific Ocean [Liu et al., 2006]. Moisture transport from the Pacific and Atlantic oceans to the South American continent can be seen as the atmospheric link between sea surface temperatures (SSTs) and land total water storage. Interannual to decadal perturbations of the Walker circulation (e.g., by the El Niño-Southern Oscillation in the Pacific Ocean) or of the Hadley circulation (e.g., the meridional motion of the Intertropical Convergence Zone) impact moisture advection, rainfall, and eventually, terrestrial water storage, as precipitation is the first-order driver of water storage changes in equatorial regions [e.g., Syed et al., 2008, Crowley et al., 2008.
[3] At interannual time scales, the El Niño-Southern Oscillation (ENSO) affects land monsoon precipitation worldwide [Zhou et al., 2008]. Its correlations with Amazon rainfall [e.g., Ropelewski and Halpert, 1987;Zeng, 1999;Liebmann and Marengo, 2001;Ronchail et al., 2002;Molinier et al., 2009], streamflow [e.g., Richey et al., 1989;Dettinger et al., 2000;Dai et al., 2009;Molinier et al., 2009], and mass balance of tropical Andes glaciers [Rabatel et al., 2013] have been highlighted. In addition to the equatorial Pacific, the tropical Atlantic has been identified as another oceanic pool influencing Amazon rainfall at interannual to decadal time scales. The ITCZ position has been found sensitive to the strength and sign of the meridional SST gradient in the tropical Atlantic [e.g., Liebmann and Marengo, 2001;Ronchail et al., 2002;Molinier et al., 2009;Yoon and Zeng, 2010] and consequently has an impact on the intensity of the South American summer monsoon [e.g., Bird et al., 2011]. Anomalously high SSTs in the tropical North Atlantic have been responsible for the recent 2005 and 2010 droughts [Marengo et al., , 2011Zeng et al., 2008a;Yoon and Zeng, 2010;Espinoza et al., 2011]. Since both Atlantic and Pacific SSTs have an influence on the climate of South America, they have been used to simultaneously forecast fire severity in the Amazon basin [Chen et al., 2011].
[4] From decadal to multidecadal time scales, the Pacific Decadal Oscillation (PDO) [Mantua et al., 1997;Mantua and Hare, 2002] influences the climate not only in North America [McCabe et al., 2004] but also in the southern hemisphere, as evidenced by tree ring and reef coral paleorecords [Evans et al., 2001], and especially in South America, according to precipitation and surface temperature records [Dettinger et al., 2000;Mantua and Hare, 2002;Marengo, 2004] or as shown by reconstructions of land water storage from river stage data .
[5] In addition to natural climate variability, human-driven climate change has an impact on the climate and hydroclimatology of the Amazon basin. Models predict more frequent droughts, due to warming North Atlantic SSTs and a reduction of northern hemisphere aerosol pollution [Cox et al., 2008], as well as an increase in precipitation extremes [Kitoh et al., 2013]. Also, the frequency and strength of Central-Pacific type of El Niños [Yu and Kao, 2007;Kao and Yu, 2009] have increased in the last 30 years [Lee and McPhaden, 2010] and are likely to do so in the future [Yeh et al., 2009;Kim and Yu, 2012].
[6] Assessing the level of ecological stress of the tropical forest requires the knowledge of the hydrological stress, i.e., the ground water storage [Toomey et al., 2011], the net balance between cumulative rainfall, evapotranspiration and runoff. Terrestrial water storage is an important component of the global water cycle, but, especially for the Amazon basin, is poorly represented by land surface models that do not represent groundwater extraction by tree deep roots and have a limited storage capacity. The former weakness may affect model predictions at interannual and decadal scales, while the latter may be an issue at decadal and longer time scales, as ground water storage multidecadal variability may significantly exceed both seasonal and interannual variability [Zeng et al., 2008b].
[7] The knowledge of the state of hydrological stress in basins worldwide relies on the durability of the Gravity Recovery and Climate Experiment (GRACE) mission. Since its launch in April 2002, GRACE has provided unprecedented gridded observations of time-variable terrestrial water storage anomalies (TWSAs) over land [Tapley et al., 2004;Wahr et al., 2004]. The decade-long GRACE record enables the study of interannual to decadal variability in land water storage. Interannual variability of land water storage from GRACE has been highlighted at different spatial scales [Andersen et al., 2005;de Viron et al., 2006;Morishita and Heki, 2008;Seitz et al., 2008;Chen et al., 2010b;Sasgen et al., 2010;García-García et al., 2011;Llovel et al., 2011;Boening et al., 2012;. In most of the previous studies, ENSO was found to be significantly correlated to water storage [de Viron et al., 2006;Sasgen et al., 2010;Morishita and Heki, 2008;Chen et al., 2010b;García-García et al., 2011;Boening et al., 2012]. In the Amazon basin, the 2005 drought in the central and western parts of the Amazon basin has been identified using GRACE [Chen et al., 2009;Xavier et al., 2010], as well as the 2002-2003 drought and the 2009 floods, which have been linked to ENSO [Chen et al., 2010a]. Impacts of the 2005La Niña and the 2006 El Niño have been also highlighted [Morishita and Heki, 2008;Xavier et al., 2010], as well as the transition from the 2009-2010 El Niño to the 2010-2011 La Niña on terrestrial water storage and global ocean mass from GRACE [Boening et al., 2012]. However, while the link between Amazon basin's water storage and Pacific SSTs has been identified in previous studies, the influence of Atlantic SSTs on land water storage has been widely ignored so far. One of the goals of this study is to highlight this influence.
[8] As climate interannual variations are intrinsically nonlinear, so are those of induced hydrological observables such as water storage. Methods that are able to capture such nonlinear behaviors in short and noisy time series should therefore be preferred to a traditional harmonic analysis. Principal component analysis (PCA) extracts variability modes that are orthogonal in both space (empirical orthogonal functions, EOFs) and time (principal components, PCs). The PCA [Preisendorfer, 1988] has been applied to up to 6 years of GRACE data [de Viron et al., 2006;Rangelova et al., 2007;Schmidt et al., 2008;Xavier et al., 2010;. de Viron et al. [2006] found a good correlation to the Southern Oscillation Index worldwide. Schmidt et al. [2008] highlighted the presence of a quasi-biennial oscillation worldwide as well as in the Amazon basin. However, the short records used in the aforementioned studies did not include many El Niño-La Niña cycles. Moreover, traditional, noncomplex PCA does not have any frequency resolution, neither does it account for phase lags between grid cells. Also, it suffers from mathematical artifacts due to the orthogonality constraints, which need to be overcome by rotating either the EOFs or the PCs [Richman, 1986;Dommenget and Latif, 2002]. Recently, the independent component analysis (ICA) method has been used to extract statistically independent components from GRACE and has performed a better separation than that obtained by rotating either EOFs or PCs . A representation of propagating oscillations can be achieved by the complex PCA, also called complex EOF analysis [García-García et al., 2011], as well as the multichannel singular spectrum analysis (MSSA) [Plaut and Vautard, 1994;Ghil et al., 2002].
[9] In this study, we use the MSSA for its ability to isolate trends and nonlinear oscillations of distinct periodicities in noisy and short data sets. We want to test whether the interannual to decadal pseudo-oscillations and trends in GRACE-TWSAs can be linked to known climate variability.
The MSSA is indeed capable of accounting for lags due to teleconnections between SSTs and TWSAs, as well as between TWSAs of different regions worldwide. While Rangelova et al. [2012] applied this technique globally, we choose to work in the spatial domain, and focus on tropical South America. Our MSSA also includes different SST-and SLP-based climate indices as additional channels, in order to highlight modes of covariability and teleconnections between GRACE-TWSAs on one hand, and SSTs/SLPs in different pools of the Pacific and Atlantic oceans on the other. After presenting the GRACE-TWSAs and climate indices in section 2, we present in section 3 our methods based on the MSSA technique and on univariate lag correlations between TWSAs and SSTs/SLPs. Results are presented in section 4 and discussion follows in section 5.

Study Area and GRACE Total Water Storage Anomalies
[10] Our study area spans 82.5°W-46.5°W, 22.5°S-13.5°N, encompassing the Amazon River basin, the neighboring basins, and the coastal and mountainous areas of Ecuador and Peru (see Figure 1 for more detailed geographical locations).
[11] We use the GRACE-Release 2 solutions computed by the Groupe de Recherche en Géodésie (GRGS) [Bruinsma et al., 2010]. The period of the study spans 10 years, from August 2002 to July 2012. These solutions have an approximate spatial resolution of 400 km at the equator and represent the sum of surface water, soil moisture, and groundwater, expressed as equivalent water height deviations toward the mean value for the study period and named terrestrial water storage anomalies (TWSAs). Uncertainty of retrieved TWSAs is about 3 cm of equivalent water height (given by the RMS of the GRACE variations over the open-ocean domain within the same latitude range). No postprocessing (e.g., decorrelation, smoothing) is performed because the regularization applied in the inversion of the GRACE observations aims at reducing the noise without the need of additional filtering techniques [Bruinsma et al., 2010]. The 10 day fields are linearly interpolated in the gaps and decimated to a monthly time series of 120 months. In order to reduce redundant information as well as to avoid undersampling, the TWSAs solutions are computed on a 3°× 3°grid (~330 km × 330 km at the equator) from the Stokes coefficients of harmonic degrees 1-50. The resulting 112 time series are plotted in Figure 1a.
[12] As we are interested in the correlation between land water storage and climate, using scaled GRACE products (as those available on the GRACE Tellus website: http://grace.jpl.nasa. gov/data/) would not change our results. The scaling process does not add any temporal frequency to the unscaled time series, because the scaling factor is frequency-independent by construction [Landerer and Swenson, 2012]. However, localized precipitation/water storage gradients, e.g., due to high topography gradients, will be missed at such a low resolution, e.g., in the coastal regions of Peru where a strong gradient of precipitation occurs between the dry plains and the wet plateaus.

Climate Indices
[13] We use eight monthly climate indices available from the NOAA's Earth System Research Laboratory Web site The main watersheds' contours are plotted in black as a background: the Amazon River basin in the center, surrounded by (counter clockwise) the Araguaia-Tocantins River basin to the east, the Maroni, Courantyne, and Essequibo River basins to the northeast, the Orinoco River basin to the north, the Magdalena River basin to the northwest, the coastal and mountainous areas of Ecuador and Peru to the west, and the Lakes Titicaca and Poopo System to the southwest. Main rivers and tributaries are plotted in light blue. Gray shading represents glaciers from the World Glacier Inventory database.
(www.esrl.noaa.gov/psd/data/climateindices/list/) and based on instrumental sea surface temperatures or pressure anomalies. They are related either to the equatorial Pacific, the North Pacific or the tropical Atlantic. They are plotted in Figure 2a and described in Table 1

. El Niño-Southern Oscillation
[14] The El Niño-Southern-Oscillation (ENSO) is a global-scale climate oscillation involving a coupling in the ocean-atmosphere system and resulting in a zonal displacement of the Walker circulation. El Niño (La Niña) conditions correspond to anomalously warm (cold) sea surface temperatures in the eastern equatorial Pacific accompanied by weaker (stronger) easterly winds. The characteristic time scales of ENSO lie in the 2-7 year band, at two dominant periods: a quasi-biennial period of~2.5 years and a quasi-quadrennial period of~4.8 years [Rasmusson and Carpenter, 1982;Barnett, 1991;Ghil et al., 2002;Yu, 2005;Kao and Yu, 2009]. Two different types of ENSO have been distinguished: a central Pacific (CP)-type and an eastern Pacific (EP)-type [Yu and Kao, 2007;Kao and Yu, 2009]. Since the late 1970s, the EP-ENSO has been characterized by a dominant quadrennial periodicity while CP-ENSO has exhibited a dominant biennial period. In addition to having different dominant periodicities, CP-ENSOs and EP-ENSOs have different structures, dynamics, and teleconnections [Kao and Yu, 2009;Li et al., 2011]. Since the 1990s, CP-ENSOs have been more frequent than EP-ENSOs [Lee and McPhaden, 2010], possibly because of a strengthening of the Hadley circulation in the central Pacific [Yu et al., 2012a]. One of the objectives of this study is to determine whether or not GRACE TWSAs are sensitive to the different ENSO types. We consider different Niño indices computed by the NOAA's Climate Prediction Center (CPC), corresponding to SSTs averaged on different overlapping oceanic pools or Niño regions, ranging from the eastern tropical Pacific (Niño 3 or EP-ENSO region) to the central tropical Pacific (Niño 4 or CP-ENSO region). In addition to the oceanic aspect of ENSO, we also consider its atmospheric component by considering the Southern Oscillation Index (SOI, also from NOAA's CPC), that is, the sea level pressure gradient across the Southern Pacific Ocean, between Darwin and Tahiti.

Pacific Decadal Oscillation
[15] The Pacific Decadal Oscillation (PDO) index was computed at the Joint Institute for the Study of the Atmosphere and Ocean at the University of Washington (http://jisao.washington.edu/pdo/PDO.latest) and was downloaded from the NOAA's Earth System Research Laboratory Web site. In the North Pacific, it is characterized by a dipole pattern in SST, opposing the eastern regions to the central and western regions [Mantua et al., 1997]. In the eastern equatorial Pacific, PDO shares the same spatial pattern as the canonical El Niño SSTs, with warmer temperatures during a warm phase [Mantua et al., 1997]. Thus, PDO is correlated to ENSO at interannual scales, but contains more variability at decadal scales ( Figure 2a). It has been found indeed that warm (cold) PDO phases favor El Niños (La Niñas) . PDO exhibits periodicities ranging from 50 to 70 years and from 15 to 20 years [Minobe, 1997;Chao et al., 2000;Mantua and Hare, 2002]. Even though the GRACE period is too short to detect a full PDO cycle, a phase shift can still be detected. In late 2007, simultaneously to the initiation of a La Niña in the equatorial Pacific, the PDO that had been mainly positive since the late 1970s switched sign. It has remained negative most of the time since then, thus exhibiting a clear negative trend during the GRACE period ( Figure 2a [Enfield et al., 1999]. TNAI exhibits interannual as well as multidecadal variations. It is highly correlated (r 2 = 0.60) with the Atlantic Multidecadal Oscillation (AMO) [Enfield et al., 2001] that represents SST changes in the whole North Atlantic Ocean. AMO/ TNAI phases having been mainly positive during the GRACE record, we only consider a possible influence of AMO/TNAI on TWSAs at interannual scales.

Atlantic Meridional Mode
[17] The Atlantic Meridional Mode (AMM) is a mode of interannual to decadal variability in the coupled atmosphereocean system in the tropical Atlantic Ocean [Chiang and Vimont, 2004]. It reflects variability of the meridional sea surface temperature gradient whose TNAI is the oceanic northern component. The AMM spatial pattern in SSTs is a dipole, with larger amplitudes in the eastern tropical North Atlantic and South Atlantic. Stronger convergence (divergence) of surface winds and more (less) rainfall happen in the warmer (colder) hemisphere due to a meridional shift of the ITCZ toward the warmer (colder) hemisphere [Chiang and Vimont, 2004]. Details of how this index is computed are provided in Table 1 from Chiang and Vimont [2004]. The index was originally computed at the University of Wisconsin (http://sunrise.aos.wisc. edu/~dvimont/MModes/Data/AMM.txt). As for the TNAI, no decadal variability in AMM may explain the long-term variation in GRACE-TWSAs.

Multichannel Singular Spectrum Analysis
[18] We use the multichannel singular spectrum analysis (MSSA) [Plaut and Vautard, 1994;Ghil et al., 2002]. This purely statistical method extracts modes of covariability among multivariate time series called channels. A technical description of the method is provided in the Appendix A.
Here are listed some properties of the MSSA that are especially interesting for this study: [19] 1. MSSA is a data-adaptive, nonparametric method that does not introduce a priori information on the input signal.
[20] 2. Most of the covariance can be explained by a few leading modes, while the remaining modes are considered as noise.
[21] 3. The extracted modes are orthogonal (i.e., uncorrelated) both in time and space, although the input channels may be correlated to each other.
[22] 4. The extracted modes can be nonlinear, such as amplitude-or phase-modulated oscillations, which is consistent with the intrinsically nonlinear climate variability.
[23] 5. MSSA is particularly well suited to the study of short data records for which a spectral analysis inevitably creates artifacts.
[24] We consider 119 channels, among which 112 channels are TWSAs time series covering the study area and seven channels are climate indices. Each channel is N = 120 month (10 years) long. For each channel, the temporal mean over the GRACE period is subtracted to each monthly value. By using climate indices as MSSA channels in addition to the GRACE-TWSAs, we aim to extract the part of each climate index variability that is correlated to the GRACE-TWSAs (with possible lags), thus to highlight possible teleconnections between SSTs and TWSAs at different time scales. Each climate index is rescaled so that its new variance equals the mean variance of the TWSAs channels over the studied area.
[25] The only user-defined parameter is the embedding dimension or lag-window length M that enables the method to account for phase lags between channels. A rule of thumb assumes that MSSA is able to separate periods in the range [M/5 M]. M also limits the spectral resolution that is defined as max(1/M, 1/(N À M + 1)), N being the number of time samples of each channels. M is set to M = N/2 = 60 months in order to reach the best spectral resolution, i.e., 1/M = 1/60 months À1 . Thus, periods in the range (1-5) years can be isolated from each other. Outside of this range, small periods (e.g., the semiannual period) are likely to be mixed with long periods. Should periods larger than 5 years be found, their statistical significance would be questionable, because less than two full cycles would be observed.
[26] MSSA provides spatiotemporal modes ranked by decreasing order of their explained cross covariance. PCs do not carry phase information, whereas EOFs do. Reconstruction is needed to reintroduce the phase information into the PCs, following equation (A5). Summing up every reconstructed modes leads to the initial time series. MSSA can be used to filter the input channels, by selecting some modes and summing them up while discarding the others. We discard the annual and subannual modes and their overtones, as we are only interested in the interannual variability. We then focus on the remaining leading modes in the eigenspectrum, as the tail essentially contains noise.

Lagged Cross-Correlation Analysis Between TWSAs and Climate Indices
[27] For each grid cell a lagged cross-correlation function is computed between the MSSA-filtered TWSAs and each of the MSSA-filtered climate indices for different leads/lags comprised between ±8 months, at a monthly step. Lagged cross correlations between climate indices are similarly computed in order to highlight correlations among indices. We compute the 5% significance level and the 95% confidence interval for every lagged correlation coefficient r. We look for local extrema of the correlation function for a given index leading or lagging TWSAs. If a local extremum is found and r is significant, the bound of the 95% confidence interval closer to zero is plotted. Further in the analysis, we introduce a threshold of ±0.7 (i.e., more than 50% of the variance of TWSAs is explained by a given index) for r. Such a high threshold is used because the data have been low-pass filtered as described in section 3.1.

Decomposition Into Pseudo-Oscillations and Trend Modes
[28] The singular spectrum of the spatiotemporal modes ( Figure 3) can be split into three parts. The first part is composed of a pair of the two leading modes having equal variances and explaining the main part (79%) of the total cross covariance. The second part consists of a suite of modes of smaller variance (between 4.1% and 1.9%, up to mode 13) reaching a plateau. From mode 14 onward, the variance constantly decreases, forming the spectrum tail.
[29] The space-time PCs are shown in Figure S1 (supporting information). Reconstructed components/modes (RCs), obtained by introducing the phase information carried by the EOFs (following equation (A5), Appendix A), are shown in Figure S2 (supporting information) and Figures 5 and 6. MSSA modes that are related to a steady or transient oscillation come by pairs. If each autocorrelation function of the PCs displays at least two side lobes, a pair is considered as an oscillation. The period of the oscillation equals the lag between two consecutive peaks of the lagged autocorrelation function of each PC (Figure 4). An oscillation is considered as steady if at least the two side lobes of the autocorrelation function have an amplitude larger than 0.5. Otherwise, the oscillation is transient, i.e., its amplitude may change with time. According to this criterion, modes 1-2 (annual), 11-12 (semiannual), 7-8, and 10 and13 (~1.5 years) are clear steady oscillations unlike modes 3 (trend/decadal) and 6 (long-period transient) and modes 4-5, 9, and 14-20 (quasi-oscillations).
[30] Although with amplitudes smaller than 0.5, the side lobes of the autocorrelation functions of PCs 4-5 are well defined at lag ±30 months, corresponding to a quasi-biennial period of 2.5 years. The rather small side lobes reflect a modulation of the oscillation amplitude, as confirmed by the increasing amplitude of the reconstructed modes from 2006 onward (Figures 5, 6, and A2). This is also found in some extent with modes 9 and 14 (quasi-period of 3 years, with increasing amplitude) and 17-18 (quasi-biennial mode at 2.4 years, with decreasing amplitude).
[31] Modes 4-5 are well isolated-along with mode 6from the other modes in the singular spectrum ( Figure 3). Even though their individual ranks are lower than mode 3, their combined explained covariance is actually the second largest (4.6%), after that of the annual oscillation (79%) (Figure 4 and Table 2). Their quasi-biennial periodicity has been found in ENSO [Rasmusson and Carpenter, 1982;Barnett, 1991;Ghil et al., 2002;Yu 2005;Kao and Yu, 2009]. It has been detected in global precipitation patterns [e.g., Ropelewski and Halpert, 1987]  streamflow observations at Manacapurú [Richey et al., 1989] and in the Magdalena River [Restrepo and Kjerfve, 2000]. Dettinger et al. [2000] showed that tropical South America's streamflows are correlated to ENSO over all time scales, but with enhanced coherency in the quasi-biennial, quasi-quadriennial, and decadal bands. Similarly, 30 year long water storage simulations by the WaterGap hydrological model in the Amazon basin showed a 24-30 month spectral peak that was associated with ENSO [Guentner et al., 2007]. The quasi-biennial oscillation explaining most of the cross covariance at interannual time scales means that TWSAs and SSTs/SLPs are highly coherent at this period. As CP El Niños (in 2002-2003, 2004-2005 have been more frequent than EP El Niños (in 2006 during the GRACE record [Lee and McPhaden, 2010;Yu et al., 2012b], only the CP-ENSO's characteristic period shows up in our analysis while the EP-ENSO's quadrennial period [Kao and Yu, 2009] is not detected. [32] Mode 3 is well isolated in the singular spectrum with a 4.1% explained variance. It is not an oscillation but rather a linear trend as suggested by the two negative lobes. A 14 year period decadal oscillation might explain the reconstructed mode ( Figures 5, 6 and A2). However, edge effects may be responsible for such a distortion of the linear trend, as the ratio of the record's length to the analyzing window length is small (10/5 = 2). Mode 6 comes along with the quasi-biennial oscillation (modes 4-5) in the eigenspectrum. When reconstructed, mode 6 looks like a long-period transient oscillation vanishing around 2006, when the amplitude of modes 4-5 starts to increase (Figures 5, 6 and A2). For each channel, similar trends are found in modes 3 and 6, but the former only occurs in the first half of the GRACE period.
[33] Since we only consider interannual to decadal variations in this study, we filter out the annual mode (1-2), the semiannual mode (11-12) and its overtones (modes 10 and 13). Modes with a rank higher than 19, forming the singular spectrum's tail, are also filtered out, as they are likely to be contaminated by noise. Eventually, only 12 modes (forming seven mode combinations) are kept to reconstruct a low-pass filtered version of the initial data (full circles in Figure 3). Together, they explain 17% of the total cross covariance, against 82.7% explained by the annual and semiannual signals and overtones ( Table 2). The reconstructed, filtered time series are plotted in Figures 1b and 2b and A2) and the GRACE-TWSAs (Figures 6a and A2) based on the seven interannual mode pairs only. In each channel   Figure S1 of the supporting information) grouped by pairs whenever they represent a (quasi-)oscillation. In each plot, the modes are colored in numerical order (blue, green, and red). For each quasi-oscillation, the characteristic period per mode combination is estimated from the lag between two consecutive peaks (above the threshold of 0.5 shown by the horizontal dotted line) and is given in Table 2. The fraction (in percentage of the total) of cross covariance explained by each (quasi-)oscillation or mode is provided in parenthesis.
(climate index or GRACE), we look at the relative importance of each mode combination, as given by the ratio of its variance to the total variance of the combined 12 modes (Figures 5b and  6b). By definition, for a given channel, the different EOFs are orthogonal, whereas the reconstructed modes are not, because the phase information carried by the PCs is added back during the reconstruction. So the sum of the all the variance fractions is not 100% in Figures 5b and 6b. Trend/decadal (3) Quasi−biennial (4−5) Long−term transient (6) 1.5 yr (7−8) 3 yr (9&14)  (3) Quasi−biennial (4−5) Long−term transient (6) 1.5 yr ( Table 2 for the characteristic period of each pseudo-oscillation. separation between interannual and longer-period changes. In the Atlantic-related indices, the quasi-biennial mode reproduces the warm events of 2010 and the cold events of 2003, 2009, and 2012, but modes 9 and 14 (3 year quasi-oscillation) are needed to modulate the amplitude of the events (especially in 2010) and reproduce the 2005 warm event. The quasi-biennial mode explains the most part of each index filtered variance ( Figure 5). Its relative contribution is the highest in TNAI. In Pacific Ocean SSTs, its relative variance fraction increases as the equatorial Pacific pool moves eastward from the CP region (Niño 4) to the EP region (Niño 3), while the relative fraction of the trend/decadal mode (mode 3) decreases and that of the long-period transient (mode 6) is constant.
[36] The trend/decadal mode (mode 3) is significant only in the Pacific SSTs and consists in a decreasing trend ( Figure 5a). As expected, its relative importance is the largest in the PDO (Figure 5b). Similarly, the long-period transient (mode 6) does not contribute much to the variance of Atlantic SSTs and its relative importance is the largest in the PDO. On the other hand, the relative contribution of the 3 year quasi-oscillation (modes 9 and 14) is more important in the Atlantic SSTs than in the Pacific SSTs.
[37] A given mode of variability is not specific to only one climate index but is usually split among different indices, due to correlations among them. For example, the quasi-biennial mode is found in all indices, but with a different phase. What makes an index different from another is its distribution or spectrum of the fractions of explained variance. The spectrum of Pacific SSTs is indeed very different from that of Atlantic SSTs, because in the latter spectrum, the contribution of the trend mode (mode 3) is lower and that of the quasi-biennial mode (modes 4-5) is higher. By looking for similarities in the mode, relative distributions among channels may enable us to identify teleconnections between indices and TWSAs within a region.
[38] Following this reasoning, we consider Figure 6, the counterpart of Figure 5 for GRACE-TWSAs. The trend/decadal mode (mode 3) is particularly strong in the southern half and coastal regions of Peru, in the southern Amazon regions and in the downstream part of the Tocantins basin and moderately strong in the northeastern drainages, northeastern Amazon regions (Rio Branco subbasin), and in the Magdalena River basin. It always pairs up with the quasi-biennial oscillation (modes 4-5) and the long-period transient (mode 6), although the opposite is not true. It is noticeably small in the western, central, and southeastern Amazon regions as well as in the Orinoco River basin and along the Pacific coast from Colombia to Ecuador. In all these regions, the variance is dominated by the quasi-biennial oscillation. Remaining quasi-oscillations (modes 7-8, 9 and 14, 15-16, and 17-18) explain relatively tinier parts of the total variance and do not exhibit clear spatial patterns (Figures 5 and 6).
[39] To summarize, the northeastern Amazon and the lowlands of the Amazon delta exhibit a variance fraction distribution similar to that of ENSO-related indices, which indicates a favored teleconnection with the equatorial Pacific Ocean. Coastal regions of Peru, the mountains of southern Peru, and the south of the Amazon basin show more similarities with the PDO's variance spectrum, i.e., exhibit longer periods, suggesting a preferred teleconnection with the North Pacific Ocean, or simply ongoing trends, as it will be discussed later in section 5.2. The western, central, and southeastern Amazon Figure 6. Same as Figure 5, but for GRACE-TWSAs. For every grid cell, the Y axis range matches the grid cell height. It varies from cell to cell in Figure 6a, while it is fixed to [0 70]% as indicated for the lower right cell. In Figure 6a), X axis spans the study period August 2002 to July 2012. Refer to Figure 5 for the name of each mode and to Table 2 for the characteristic period of each pseudo-oscillation. The main watersheds' contours are plotted in black as a background: the Amazon River basin in the center, surrounded by (counter clockwise) the Araguaia-Tocantins River basin to the east, the Maroni, Courantyne, and Essequibo River basins to the northeast, the Orinoco River basin to the north, the Magdalena River basin to the northwest, the coastal and mountainous areas of Ecuador and Peru to the west, and the Lakes Titicaca and Poopo System to the southwest. Gray shading represents glaciers from the World Glacier Inventory database. regions as well as the Ecuador, the Columbian Pacific coast, and the Orinoco River show a variance fraction distribution that is closer to those of the TNAI and AMM, suggesting a favored teleconnection with the tropical Atlantic Ocean. An analysis of the correlation between TWSAs and each climate index will aim at quantifying these qualitative findings (see section 4.3).

Time Evolution of Reconstructed TWSAs Seasonal Averages
[40] Figure 7 shows the quarterly averages of the MSSAfiltered TWSAs for each year of the 10 year period, which are obtained by summing up the different reconstructed modes ( Figure S3). When only the quasi-biennial mode is reconstructed ( Figure S3a [41] During El Niños, negative TWSAs occur in boreal winter and spring, lasting two seasons; during La Niñas, positive TWSAs appear in boreal winter too, but are more persistent, lasting up to 1 full year (e.g., in 2008-2009 and in 2011-2012, even though its effect was combined with that of a weak La Niña). This reflects the well-known asymmetry between the two phases of ENSO. During ENSO events, extreme TWSAs are found in the northern and northeastern basins and in the lower Amazon. However, during the only EP-El Niño of the GRACE period, in 2006-2007, the negative TWSAs are centered further southwest, affecting only the Amazon basin, in comparison to those found during the CP-El Niño events of the study period ( Figure S3c). This suggests an intrinsic difference in the effects of the two El Niño types on TWSAs, as highlighted in rainfall and temperature patterns by Li et al. [2011]. Their results are consistent with our findings of negative TWSAs located further south in the Amazon basin during EP-El Niños.  Figure S3 for the mode-wise reconstruction.
[42] During AMM events, extreme TWSAs happen in the central western parts of the Amazon basin in boreal summer and fall, but can also persist after the end of the events: for example, TWSAs in 2011 are still negative because of the memory effect of the previous dry year.
[43] Finally, after the trend/decadal mode (mode 3) is added (Figure 7), anomalies both at the beginning and at the end of the GRACE period become stronger: the northeastern drainages, the lower Amazon basin (Rio Branco, and lower Rio Solimões), as well as the southern part of the Amazon basin, are wetter in 2012 than they are in 2002, while the western and central parts of the Amazon River basin (headwaters of the Ucayali, Juruá, Madre de Dios, and Purus Rivers) and the southern Peru (Andes' Cordillera Oriental, centered at 72°W-15°S) have become drier, at least until the end of 2011. Thus, adding the trend mode strengthens a trend already existing in the interannual modes. A clear separation between the long-term trend and interannual variability is not achieved because both time scales share the same spatial patterns as shown in Figure 8: the trend/decadal mode (Figure 8b) is spatially correlated with the quasi-biennial mode (Figure 8c) and the long-period transient mode (Figure 8d). A longer record is therefore desirable in order to better disentangle interannual and decadal variabilities.

Teleconnections Between Reconstructed Climate Indices and TWSAs
[44] We now study the lag correlation between the MSSAfiltered TWSAs channels and each of the MSSA-filtered climate indices. For each channel, MSSA filtering consists in summing up the seven reconstructed mode combinations (made up of 12 modes) studied above, which leads to the low-pass filtered time series of Figures 1b and 2b. Lead/lag correlation is performed as explained in section 3.2. Largest correlation coefficients and corresponding lags (i.e., TWSAs lagging a climate index) are displayed on the map when a statistically significant local extremum is found in the correlation function (Figure 9). We focus on the |r| ≥ 0.7 values to compute the statistics provided in Table 3.   Table 3. See Figure S4 (supporting information) for lagged correlations with GRACE-TWSAs leading the climate indices.
Largest correlation coefficients with negative lags (i.e., TWSAS leading a climate index, Figure S4 of the supporting information) are usually smaller with longer lags than when TWSAs lag a climate index, except with the PDO as explained in section 4.3.2.

TWSAs Teleconnection With ENSO
[45] Teleconnections with the different ENSO-related indices are shown in the first four subplots of Figure 9. Overall positive (respectively negative) correlation coefficients are found between SOI (respectively the Niño indices) and GRACE, meaning that there is larger (smaller) total water storage during La Niña (El Niño) events. Cells with the largest correlation are found in the northeastern part of the Amazon River basin (Rio Branco and lower Rio Solimões), in French Guyana, Suriname, in the Essequibo River Basin (Guyana), lower Orinoco River basin (Venezuela), and Magdalena River basin (Colombia). It is noteworthy that no extremum of the lag correlation function is found in the western and central parts of the Amazon basin. Smaller, negative (respectively positive) correlations to SOI (respectively Niño indices) are found in the southern Amazon, and small positive (negative) correlations are found in the Peruvian southern coastal regions and Andes' Cordillera Oriental. Our spatial patterns of correlation in the northeastern drainages and northeastern Amazon basin are consistent with those found in, e.g., Ropelewski and Halpert [1987], Mantua and Hare [2002], Landerer et al. [2008], and Becker et al. [2011].
[46] Among the ENSO-related indices, fewer cells pass the |r| ≥ 0.7 threshold with the Niño 3 index (20 cells versus 24-26 cells for the other indices). Mean correlation between GRACE-TWSAs and Niño 3 is also the smallest (À0.73 versus À0.77, À0.79, or 0.80). The significance of the difference between two r mean values was tested using a two-sample, right-sided t test with unknown mean and standard deviation (null hypothesis: |r 1 | À |r 2 | = 0, against the alternative: |r 1 | À |r 2 | > 0, where r 1 and r 2 are the mean r values obtained with the tested pair of indices, with |r 1 | > |r 2 |). Except for the mean r found with Niño 4 that is not significantly different with 95% confidence from those obtained with Niño 3.4 and SOI (p = 0.0508 and 0.32 respectively), all other mean r differences are significant (p < 0.01), especially those obtained between Niño 3 and any other index. This latter result suggests that during 2002-2012, SSTs in the equatorial CP region are more correlated than those in the EP region to tropical South America TWSAs.
[47] The spatial average of the time lag between TWSAs and the ENSO-related indices over the cells with |r| ≥ 0.7 equals 2-3 months (Table 3). However, as shown by Figure 9b, the time lag varies spatially, with the largest lags (5-7 months) found in the lower Amazon while it decreases to 0-1 month northwestward. Since equatorial Pacific SST peak in early boreal winter (December) during ENSO, the peak in TWS occurs in December (no lag) in the Magdalena basin (Colombia), in February (2 month lag) in the northeastern drainages and Orinoco basin, and finally in May-July (5-7 month lag) in the lower Amazon, as confirmed by the time evolution of seasonal averages in Figure 7. According to the timing of the rainy season in those regions [e.g., Restrepo and Kjerfve, 2000;Ronchail et al., 2002], the peak in TWS due to ENSO happens between the end of the wet season and the beginning of the dry season, lagging the peak of precipitation by~2 months. Therefore, ENSO's maximal impact on TWS occurs when TWS is high (wet season for TWS). This means that, in the northeastern regions, El Niño's impact on vegetation may become obvious later, during the dry season. However, the asymmetry between the ENSO phases (i.e., La Niña events lasting longer than El Niño events, as seen in Figure 7) may introduce some deviation toward these time lags, which are averaged lags for both phases of ENSO.

TWSAs Teleconnection With PDO
[48] Teleconnections with PDO are shown in the fifth subplot of Figure 9. Areas with negative correlation between GRACE-TWSAs and PDO are the same regions that show positive correlation to SOI and vice versa. SOI and PDO are indeed strongly anticorrelated (r = À0.88 at 3 month lag, Table 4), especially at interannual scales. PDO containing more decadal variability than SOI correlation coefficients are expected to be larger with PDO than with SOI, when the ratio of TWSAs' decadal variability to interannual variability is high (see section 5.2).
[49] Twenty-four grid cells have r ≤ À0.7, the mean correlation coefficient being significantly (p < 0.05, two-sample, For each index, number of cells, mean correlation coefficient, and mean lead or lag time are averages computed for the grid cells with |r| ≥ 0.7 within the 95% confidence interval, at the largest correlation. Values in italics correspond to the case of the index lagging GRACE. one-sided t test) larger (À0.83) than that found with SOI, Niño 4, and Niño 3.4 (0.80, À0.79 and À0.77, respectively). On average over these cells, PDO leads TWSAs bỹ 2 months. Largest time-lags (5 months) are found in the lower Amazon while they decrease to zero northward. PDO being maximal in January, the TWS response occurs from January to May. This is consistent with the timing of the TWS response to ENSO (section 4.3.1) and with the fact that PDO lags ENSO by 2-3 months (Table 4). Mantua and Hare [2002] found similar anticorrelation in northern South America between PDO and precipitation/streamflow, indicating a likely influence of the North Pacific decadal variability on the hydroclimatology of these regions. In the Amazon basin, decreasing (increasing) rainfall in the north (south) were found to be associated with positive PDO phases and more frequent El Niños [Marengo, 2004].
[50] On the other hand, two cells, located in the headwaters of the Ucayali River in the Cordilleras of east-central Peru and in the coastal part of southern Peru (at 72°W-15/18°S), show a large positive mean correlation coefficient (0.77) to PDO, with no lag. Since no very strong correlation to ENSO is found in these cells (Figure 9), we conclude that decadal variability dominates over interannual variability there, as confirmed by the dominant decreasing trend observed in the TWSAs time series in those cells (Figure 1a) (see section 5.2).
[51] Fourteen grid cells over three distinct areas (coastal Guyana and Suriname, Magdalena River basin, and southern Amazon basin) display strong anticorrelation with PDO (À0.86 on average), but with GRACE leading PDO by 1-3 months (1 month on average) as shown by Figure S4 (supporting information) and Table 3. Similarly, small lags with same anticorrelation are found when GRACE lags PDO in regions that are in the vicinity of those where GRACE leads PDO. Also, the peak of the lagged crosscorrelation function is very broad, making the choice of the lag time at the minimum correlation not very characteristic. We therefore conclude that there is no significant difference among those regions in their teleconnection to PDO. After including those regions in the statistics of Table 3, average correlation between TWSAs and PDO is higher (À0.84) and more spread out (38 cells) than with any ENSO-related index.
[52] Large correlations (0.74 on average) are found over four grid cells in the Jurua River, upstream part of Purus River, and coast and Cordilleras of central Peru. However, the 5-month average lag time is significantly larger than the 0-month lead time found further south in the Cordilleras of east-central Peru and in the southern coast of Peru, indicating different mechanisms for the teleconnections of these regions to PDO.

TWSAs Teleconnection With TNAI and AMM
[53] Teleconnections with TNAI and AMM are shown in the sixth and seventh subplots of Figure 9, respectively. Correlation patterns are similar among both indices, with strong anticorrelations in the western and central parts of the Amazon basin. During AMM and TNAI positive phases, the ITCZ lies anomalously far northward, creating an anomalous moisture divergence and droughts in the south and central parts of the Amazon basin. Correlations are similarly strong for TNAI (À0.77 over 13 cells) and AMM (À0.75 over 11 cells), although not significantly different (p = 0.11) according to a two-sample one-sided t test, which is due to the fact that both indices are highly correlated as shown in Table 4. Note that by construction, the impact of ENSO on AMM has been removed in order to have an index that is uncorrelated to ENSO at zero lag (see Table 4). Therefore, the correlation patterns for ENSO and AMM almost do not overlap.
[54] On average, TNAI leads GRACE by 3 months, while AMM leads GRACE by 1.4 months. Deviation toward these mean lag times ranges from 1 (0) to 5 (3) months for TNAI (AMM), and increases southward. The different mean lags are consistent with the fact that AMM lags TNAI by 2 months ( Table 4) and suggest that TNAI should be favored for TWSAs forecasting. TNAI and AMM peaking in February and April respectively, the associated response in TWS in the central and western Amazon peaks in May (3 month lag), at the very end of the wet season [Ronchail et al., 2002]. So, in these regions, a positive TNAI (or AMM) leads to lower-than-usual water storage in the dry season, which may have dramatic consequences on vegetation with increased fire frequency [Chen et al., 2011[Chen et al., , 2013.

Comparison With Other Proxies
[55] The influence of ENSO on Amazon rainfall was shown by Zeng [1999], Molinier et al. [2009], Yoon and Zeng [2010] and Espinoza et al. [2011] who found a 3-4 month lag between ENSO and rainfall. A 6-7 month lag was found between ENSO and discharge at Manacapuru and Obidos [Richey et al., 1989;Zeng 1999], while Zeng [1999] found a Largest correlation coefficients are highlighted in bold font. b Lag times are in brackets and are positive (negative) when the index in a row lags (leads) the index in the column.
a 1-2 month lag between rainfall and basin-averaged soil moisture in the Amazon. These results are consistent with the 4-6 (i.e., 6-7 minus 1-2) month lag we find between TWSAs and SOI in the lower Amazon. Within the basin, Molinier et al. [2009] found the highest correlations between rainfall and SOI in the northeastern Amazon, with rainfall and discharge lagging SOI by 3 and 6 months respectively. Lags of a few months are usually expected between fluxes and stores, so that land water storage (store) lags rainfall (incoming flux) and leads discharge (outgoing flux). Our results (Figure 9 and Table 3) indicate that in the northeastern Amazon regions, lag times between ENSO and TWSAs (5-7 months) are quite close to those between ENSO and discharge (6 months, as found by Molinier et al. [2009]). This may be due to the fact that surface water is the largest storage component in these regions.
In the Magdalena basin, however, no lag is found between TWS, SOI, and discharge [Restrepo and Kjerfve, 2000], which is consistent with the very small lags found between GRACE-TWSAs and SOI there (Figure 9).
[56] Li et al. [2011] found different anomalies and timings in rainfall and land surface temperatures between the two El Niño types, with negative rainfall anomalies being more spread out southward and affecting a larger part of the Amazon basin during EP El Niños. Even though our observation of EP El Niños with GRACE is limited to only one event (in 2006-2007), we find a significant difference in the spatial patterns of the negative TWSAs between the two El Niño types, which is consistent with the findings of Li et al. [2011].
[57] Our results indicate on average a 3 month lag between TNAI and GRACE-TWSAs. Yoon and Zeng [2010] found that basin-averaged rainfall was correlated to TNAI mostly in fall with no delay, while discharge at Obidos lagged TNAI by 5 months, which once again shows that TWS lags rainfall and leads discharge.
[58] Chen et al. [2011] computed lagged correlations between interannual fire intensity in South America and Pacific and Atlantic SSTs. They found spatial patterns of correlation similar to those found with TWSAs: the eastern Amazon regions being more correlated to the Pacific SSTs while the southern and southwestern regions were more correlated to Atlantic SSTs. Their time lags are a few months larger than ours, which is consistent with the fact that the hydrological drought precedes the ecological drought, thus suggesting a long chain of consequences starting with warm SSTs and ending with the degradation and destruction of the tropical forest [Chen et al., 2013;Saatchi et al., 2013]. Ecosystems have a long memory, so that the droughts' impacts on water storage and vegetation are still visible after the droughts end [Saatchi et al., 2013], which leads eventually to lower correlations to SSTs in comparison to rainfall.
[59] Correlation patterns and associated time lags between SSTs and the different observables mentioned above (rainfall, TWSAs, streamflow, and fire activity) should not necessarily be identical, as those observables are more or less spatially integrated, or that they are measured at a coarse spatial resolution. Because the Solimões River drains such a large land surface area, streamflow in the downstream Amazon regions integrates runoff from all the upland watersheds as a response to rain falling farther upstream. So streamflow is expected to have different patterns of correlation with SSTs and longer time lags [Dettinger et al., 2000] than those found with rainfall. Surface water being a large constituent of GRACE-TWSAs, more complex patterns of correlation between GRACE and SSTs are also expected.

Origin of Observed GRACE Trend: Anthropogenic or Climate Decadal Variability?
[60] Differences between correlations of SOI and PDO with GRACE are due to the stronger decadal variability embedded in the PDO. Interannual variability is overall dominant in SOI, even though it shows a negative trend because of more frequent La Niña events toward the end of the GRACE period (Figures 2 and 5). On the opposite, decadal variability dominates the PDO as evidenced by the clear negative trend during the GRACE period (Figures 2  and 5). Therefore, cells that are more correlated to PDO than SOI are likely to display a dominant trend-like/decadal variation. In this case, if the correlation to PDO is positive (negative), land water storage is expected to have decreased (increased) during the study period.
[61] However, PDO and ENSO interannual variabilities are correlated, so we cannot unequivocally attribute the origin of TWSAs interannual variations to PDO or ENSO. Because PDO moreover has decadal variability, correlation to PDO may be due to coherence at either time scale. When high correlations are found between GRACE and ENSO, similar correlations are found with PDO, while the opposite is not always true (Figure 9).
[62] A way to separate interannual and decadal variability in GRACE-TWSAs is to consider their correlation with each MSSA mode reconstructed in the PDO (Figure 10). Since mode 3 represents linear/decadal variability and modes 4-6 represent interannual and transient variability, respectively, better correlations with mode 3 are expected if TWSAs exhibit a dominant decadal variability.
[63] By comparing the subplots of Figure 10, we find that the cells located in the headwaters of the Ucayali and Madre de Dios Rivers in the Cordilleras of east-central Peru as well as in the coastal part of southern Peru (12/15°S 72°W, region A in Figure 10) are more correlated to the trend/decadal mode than to the quasi-biennial or the long-period transient reconstructed in the PDO. This suggests that water storage in those cells is experiencing a decadal mass loss. This can be confirmed by studies showing the retreat of Peruvian and Bolivian Andes glaciers in the Cordillera Oriental, such as the Quelccaya Ice Cap in southern Peru, or the Chacaltaya glacier in Bolivia that disappeared in 2009 [Jacob et al., 2012;Rabatel et al., 2013]. On the other hand, decadal water storage increase is observed for the cells located in the southern Amazon regions (Bolivian Amazon and neighboring Brazilian state of Rondônia, region B in Figure 10) where larger anticorrelations are found only with the trend/decadal mode reconstructed in PDO. Land cover changes due to agriculture and cattle farming development have been found responsible for decreased evapotranspiration and increased runoff, streamflow, and floodings in the savanna regions of Brazil (Tocantins-Araguaia basins) [Coe et al., 2011]. We may also expect an increase in groundwater recharge and storage, like that observed in Niger [Favreau et al., 2009]. Similar impacts are also expected in the southern Amazon regions corresponding to the Brazilian states of Mato Grosso and Rondônia and the El Beni department in northeastern Bolivia, the Amazon country that exhibited the highest deforestation portion in 2011 (In the Amazon, threats to an arc of wilderness, Washington Post, August 30, 2012, http:// www.washingtonpost.com/world/in-the-amazon-threats-to-anarc-of-wilderness/2012/08/30/44354750-f309-11e1-adc6-87dfa8eff430_graphic.html). We therefore suspect that the positive trend in TWSAs (Figure 8) as well as their anticorrelation to mode 3 reconstructed in PDO (Figure 10), both especially large in Rondônia and El Beni, are a consequence of deforestation and land cover change, although a contribution from climate variability or an indirect atmospheric feedback to precipitation (due to reduced evapotranspiration) cannot be excluded [Davidson et al., 2012]. However, the latter effect is not expected to reverse the sign of annual mean discharge changes in the Madeira basin, as shown by coupled land-atmosphere model simulations with different future deforestation scenarios [Coe et al., 2009].
[64] Finally, in the northeastern Amazon regions, correlations of similar magnitude as those found with PDO are obtained with the quasi-biennial and long-period transient (modes 4-6), showing that interannual variability is dominant in these regions (region C in Figure 10). Interdecadal variations of the teleconnections between ENSO and hydrology have been highlighted and associated to PDO decadal variability. In the Amazon, decreasing (increasing) rainfall in the north (south) were found to be associated with positive PDO phases and more frequent El Niños [Marengo, 2004]. Dettinger et al. [2000] showed a strengthening of streamflow/ENSO teleconnections during PDO cold phases. An interdecadal PDO phase shift from hot to cold may have happened in 2007 and might be responsible for wetter N-NE South America (triggering more La Niñas) at the end of the GRACE record compared to its beginning. However, transient phase shifts may happen without any change in the North Pacific climate according to SSTs' canonical pattern [Bond et al., 2003;Overland et al., 2008].

Impact of Climate Change
[65] According to Lee and McPhaden [2010], frequency and amplitude of CP(EP)-type El Niños have increased (decreased) in the past 30 years while the intensity and frequency of occurrence of La Niñas have been kept unchanged, leading to warmer global temperatures. Kim and Yu [2012] found the same results with coupled climate models. Such a predicted increase in frequency of CP-type El Niños is consistent with our findings, i.e., TWSAs are better correlated to SSTs in the Niño 4 region than in the Niño 3 region. However, in a long 1903-1985 record of Rio Solimões discharge, Richey et al. [1989] found more energy at 2.4 years than at 4 years, indicating that CP-ENSO already has a stronger influence on Amazon streamflow than EP-ENSO maybe because of the integrative character of streamflow in such a large basin. This would imply that the predicted increase in the frequency of CP events would not change the natural hydrological response of the Amazon. An increasing El Niño intensity would lead to more droughts in areas affected by ENSO, whereas more floods are observed due to more frequent La Niñas in the last 4 years of the GRACE record. The increasing amplitude of the quasi-biennial oscillation that we found (Figures 5, 6, and S2) suggests that both El Niño and La Niña events have become stronger, whereas Lee and McPhaden [2010] suggest that only El Niños have. The ongoing decline in the PDO may be the natural cause of the recent observed variability in GRACE-TWSAs.
[66] Li et al. [2011] found different rainfall anomalies and timings between the two El Niño types and showed that the EP-type led to a more severe water stress and a negative net ecosystem production (and subsequently a positive net carbon exchange to the atmosphere). Using a biogeochemical model, they could not predict a decrease in the net ecosystem production (and subsequently in the carbon sequestration) due to an increased frequency of CP El Niños with climate change, but rather suggested that the Amazon will remain a carbon sink in the future. More frequent and intense droughts such as the 2010 drought may therefore come from the association of warmer SSTs in the tropical North Atlantic and El Niño events in the equatorial Pacific [Ronchail et al., 2002;Espinoza et al., 2011].

Conclusion
[67] TWSAs are influenced by both Atlantic and Pacific SSTs. By performing univariate correlation analyses, we found distinct domains of sensitivity for each oceanic pool: the central-western Amazon regions for the tropical North Atlantic pool versus the northeastern Amazon regions and northeastern drainages of South America for the equatorial Pacific pool. Mean lag times of 2.5 months (ENSO) and 3.2 months (TNAI) were found. These values are smaller or equal to those of the lags found between the same indices and fire activity [Chen et al., 2011], proving that hydrological stress precedes ecological stress and showing the forecasting capability of SSTs regarding TWSAs.
[68] Among the different Pacific pools, we found that the Central Pacific pool (the so-called Niño 4 region) was significantly (p < 0.01) more correlated to GRACE-TWAs, which is corroborated by the presence of a quasi-biennial oscillation in the data.
[69] Because the GRACE record is only a decade-long, correlation between GRACE-TWSAs and the Pacific decadal oscillation may highlight climate decadal variations as well as long-term trends of anthropogenic origins. Because their effects in GRACE are consistent with independent evidence coming from in situ observations and model simulations, we suspect that anthropogenic effects are likely in the following regions: the Andes' Cordillera Oriental that is currently experiencing glacier and ice cap melting, and the southern Amazon regions of the Bolivian El Beni department and the Brazilian Rondônia State where deforestation and land use change may have increased river discharge. Decadal variability of the North Pacific climate is likely to control TWSAs decadal change in the northeastern Amazon and northeastern drainages of South America where an overall mass gain has been observed during the past decade due to strong recent La Niña events in contrast to more frequent El Niño events at the beginning of the GRACE period.
[70] As the TWS measurements will be interrupted between the upcoming end of the current GRACE mission and the launch of the follow-on mission planned in 2017, Pacific and Atlantic SSTs may be used as controlling variables in an empirical model designed to forecast TWSAs in tropical South America.
[78] Oscillations are characterized by a pair of modes. Then the EOFs of a same pair are in phase quadrature. The period of an oscillation is estimated as the period of the oscillations of the PCs' autocorrelation function (a pair is considered as an oscillation if the autocorrelations of the PCs have at least two side lobes higher than 0.5).
[79] Each mode can be reconstructed in the initial space by combining each EOF with its corresponding PC: with M t , L t , U t depend on t and are given in Ghil et al. [2002, equation (12)].
[80] As discussed by Rangelova et al. [2012], applying MSSA in the spatial domain may be less efficient than in the spectral domain of the spherical harmonics. MSSA consists of an eigendecomposition of a LM × LM matrix. With the limited size and resolution of our study area, it is however more efficient to work in the spatial domain because LM = 119 × 60 = 7,140, against LM = 2600 × 60 = 156,000 in the spectral domain with a maximum harmonic degree of 50. Our code uses one CPU on a supercomputer, but it does not run on a dual core PC. Increasing the domain size, the spatial resolution and/or the lag-window length M would dramatically increase CPU time.