Statistical prediction of terrestrial water storage changes in the Amazon Basin using tropical Pacific and North Atlantic sea surface temperature anomalies

. Floods and droughts frequently affect the Amazon River basin, impacting transportation, agriculture, and ecosystem processes within several South American coun-tries. Here we examine how sea surface temperature (SST) anomalies inﬂuence interannual variability of terrestrial water storage anomalies (TWSAs) in different regions within the Amazon Basin and propose a statistical modeling frame-work for TWSA prediction on seasonal timescales. Three simple semi-empirical models forced by a linear combination of lagged spatial averages of central Paciﬁc and tropical North Atlantic climate indices (Niño 4 and TNAI) were calibrated against a decade-long record of 3 ◦ , monthly TWSAs observed by the Gravity Recovery And Climate Experiment (GRACE) satellite mission. Niño 4 was the primary external forcing in the northeastern region of the Amazon Basin, whereas TNAI was dominant in central and western regions. A combined model using the two indices improved the ﬁt signiﬁcantly ( p < 0.05) for at least 64 % of the grid cells within the basin, compared to models forced solely with Niño 4 or TNAI. The combined model explained 66 % of the observed variance in the northeastern region,


Introduction
The Amazon River basin has experienced several severe droughts during the last decade (Marengo et al., , 2011Chen et al., 2009;Espinoza et al., 2011;Lewis et al., 2011;Frappart et al., 2012), with extreme events in 2005 and 2010 increasing forest mortality (Phillips et al., 2009;Lewis et al., 2011) and the number of satellite-detected fires (Chen et al., 2013b). More frequent extreme wet events also have been observed in the Amazon in the last 20 years (Chen et al., 2010;Boening et al., 2012;Espinoza et al., 2013;Gloor et al., 2013). Droughts and floods lead to important economic losses by affecting land and river transportation, agriculture, fisheries , and hydropower generation (Lima and Lall, 2010). They also influence seasonal and interannual variability of the regional carbon budget by modifying photosynthesis, tree mortality, fires, rates of autotrophic and heterotrophic respiration, and river-dissolved organic carbon fluxes (e.g., Richey et al., 2002;Baker et al., 2008;Phillips et al., 2009;Lewis et al., 2011;Miller et al., 2011;Chen et al., 2011;Davidson et al., 2012). Extreme hydrological events may become more frequent in a changing Published by Copernicus Publications on behalf of the European Geosciences Union.
climate due to increases in precipitation extremes (Kitoh et al., 2013) or a tendency toward longer and more intense dry seasons (Joetzjer et al., 2013). To limit economic and ecosystem impacts, an important goal is to provide managers with early warning information about drought and flood risk based on observations of key climate system predictors on seasonal timescales. In this context, a key intermediate step is developing a predictive understanding of variations in terrestrial water storage -the integral of surface water, soil moisture, and groundwater.
While the influence of the El Niño-Southern Oscillation on the Amazon's hydroclimatology is well established (e.g., Ropelewski and Halpert, 1987;Richey et al., 1989;Enfield, 1996;Nogués-Paegle and Mo, 1997;Zeng, 1999;Dettinger et al., 2000;Liebmann and Marengo, 2001;Ronchail et al., 2002;Grimm, 2003;Dai et al., 2009;Molinier et al., 2009;Lima and Lall, 2010;van der Ent and Savenije, 2013), several recent studies have identified ocean-atmosphere coupling in the tropical Atlantic as another key regulator of hydrological variability within the basin. Studies provide evidence that the severe 2005 and 2010 droughts in the central and western Amazon were the result of warmer tropical North Atlantic sea surface temperatures (SSTs) that induced a northward shift of the Intertropical Convergence Zone (ITCZ) and reduced easterlies and moisture advection from the tropical North Atlantic to southern and western regions of the Amazon during austral summer Zeng et al., 2008;Molinier et al., 2009;Yoon and Zeng, 2010;Espinoza et al., 2011). These changes in atmospheric circulation indicate that ocean-atmosphere-land dynamics in the region, including the South American monsoon (Nogués-Paegle and Mo, 1997;Zhou and Lau, 1998), interact with (and are modified by) SST variability.
Evidence for a dual external ocean forcing of the Amazon's hydroclimate also comes from studies that show both equatorial Pacific and North Atlantic SST anomalies are correlated with many variables that influence terrestrial ecosystem drought stress, including precipitation, evapotranspiration, surface relative humidity, and the number of satellitedetected active fires (Chen et al., 2011(Chen et al., , 2013a. These variables, while crucial for assessing climate impacts on terrestrial ecosystem function, are not direct indicators of variations in regional water budgets. Specifically, precipitation is only one of several variables contributing to the surface water mass balance of soils and aquifers, and fires occur as a consequence of a complex set of feedbacks between humans, ecosystem processes and climate. Additional information is needed to understand more directly how tropical climate modes influence regional variations in terrestrial water storage within the Amazon Basin. Developing accurate forecasting systems for terrestrial water storage variations and other indicators of flood and drought risk on seasonal timescales within the Amazon is challenging and likely requires a combination of dynamical and statistical approaches. Seasonal precipitation forecasts from climate models often have low to moderate skill in the region starting from a 1-month lead time, as a consequence of uncertainties in initial conditions and incomplete representation of convection and other atmospheric processes (e.g., Lavers et al., 2009). In turn, global hydrology and land surface models often have difficulties representing soil moisture in deeper soil layers (Zeng, 1999) and surface water dynamics in rivers, lakes, and floodplains (Swenson and Milly, 2006;Han et al., 2010;Werth and Güntner, 2010). Important uncertainties also remain in our understanding of biosphereatmosphere interactions, including the role of deeply rooted trees in modulating seasonal and interannual variability in regional climate (Lee et al., 2005;Golding and Betts, 2008) and the impacts of fire-emitted aerosols on radiation and circulation (Tosca et al., 2013). In this context, investments in both statistical and dynamical forecasting approaches have the potential to advance the field. For statistical models, this may occur through extraction of key mechanisms and teleconnections regulating variability in hydrological cycle processes (and also by serving as a benchmark for forecasts from other, more complex models). Dynamical model forecasts, in contrast, are likely to become increasingly robust with advances in climate modeling, increases in resolution, and improved estimates of initial conditions (Barnston et al., 2012;van Dijk et al., 2013). These models also are more likely to maintain their fidelity during periods of rapid Earth system change since they simulate the complete evolving state of ocean-atmosphere conditions.
Here we developed several statistical models to predict terrestrial water storage anomalies (TWSAs) using observations from NASA's Gravity Recovery And Climate Experiment (GRACE) mission. GRACE satellite observations have been used along with precipitation data to predict floods worldwide with a 1-month lead time, by estimating the saturation level of the land surface (Reager and Famiglietti, 2009). In the Amazon Basin, statistically significant teleconnections between GRACE TWSAs and either Pacific or Atlantic SST anomalies were identified by de Linage et al. (2013), suggesting that SST anomalies may be used as a proxy for the meteorological processes influencing the balance between precipitation and evapotranspiration.
We build on this work in the current study by developing a series of semi-empirical models that are driven by SST anomalies and used to explain variability in monthly TWSAs across different river basins and regions within tropical South America (encompassing the region between 23 • S-14 • N and 46 • W-83 • W). These models were forced by tropical Pacific or North Atlantic climate indices with variable lead times relative to the predicted TWSA time series. This allowed us to examine the potential of using SST anomalies one to two seasons ahead to predict TWSAs, taking advantage of both SST and TWSA persistence on these timescales. In Sect. 2 we describe the source and processing steps used to prepare the GRACE and SST anomaly observations used in our analysis. In Sect. 3 we discuss how we developed our statistical

Data
Our study domain encompassed the Amazon River basin and several neighboring regions: 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 Titicaca and Poopó lakes' system to the southwest (Fig. 1).

GRACE terrestrial water storage anomalies
We used the GRACE-Release 2 solutions computed by the Groupe de Recherche en Géodésie Spatiale (GRGS) (Bruinsma et al., 2010). We chose to use the GRGS solutions rather than other publicly available GRACE solutions because of their ability to preserve small-wavelength signals, due to a constrained inversion scheme that adapts to the raw data noise structure and thus improves the signal-to-noise ratio (Bruinsma et al., 2010). No post-processing was therefore performed. These solutions consist of water storage anomalies expressed in equivalent water height, that are the sum of surface water, snow, soil moisture and groundwater, at a spatial resolution of 400 km at the equator and with an uncertainty comprised between 20 and 30 mm on average over our study domain, depending on the estimation method (i.e., by computing the root-mean-square of interannual variations over desert land areas, as in Bruinsma et al. (2010), or over the open-ocean domain within the same latitude range as our study area, following Chen et al., 2009). TWSAs are the deviations from a reference value during the study period, as absolute levels cannot be measured. The 10-day fields were linearly interpolated during gaps and decimated to a monthly time step for consistency with the SST anomaly time series. The available GRACE record spanned the period from August 2002 to July 2012. The TWSAs were computed on a 3 • × 3 • grid (333 km × 333 km at the equator) from the Stokes coefficients to minimize redundant information and to avoid undersampling. In each grid cell, a mean monthly climatology constructed from the 10 year period was subtracted from the data to remove the mean annual cycle and isolate variations associated with hydrologic extremes. The residual interannual TWSA amplitudes ranged from −280 to 360 mm over the 112 grid cells and their standard deviation varied between 20 and 135 mm. The largest standard deviation occurred in the northeastern part of the basin, including the Rio Branco watershed, the lower Amazon downstream of Manaus, and the coastal basins of the Essequibo, Maroni and Courantyne rivers in Guyana, Suriname and French Guyana (Fig. 1). In regions near the mouth of the Amazon, large standard deviations of terrestrial water storage on annual timescales also are observed (Fig. 1a of de Linage et al., 2013). This spatial pattern is correlated with areas that have a large surface storage capacity as a consequence of deep river beds and extensive floodplains that often experience high fractions of inundated area (e.g., Prigent et al., 2007).

Climate indices
The El Niño-Southern Oscillation (ENSO) is a globalscale climate oscillation involving a coupling in the oceanatmosphere system in the Pacific, with complex teleconnections worldwide. Large-scale subsidence and reduced precipitation over the northeastern regions of tropical South Figure 2. Central Pacific and tropical North Atlantic climate index regions (Niño 4 and TNAI; shaded gray areas). Area-weighted sea surface temperatures anomalies from these regions were found to be correlated with interannual terrestrial water storage anomalies in tropical South America.
America occur during El Niño, because of a reduction in uplift associated with an eastward shift of the Walker circulation (e.g., Saravanan and Chang, 2000). Concurrently, some subtropical regions in South America benefit from enhanced rainfall as a consequence of strong northerly lowlevel winds that advect moisture from the Amazon Basin and are associated with a weakening of the South America Convergence Zone (e.g., Nogués-Paegle and Mo, 1997;Grimm, 2003;Carvalho et al., 2004).
Here we used the Niño 4 index from NOAA's Climate Prediction Center (www.esrl.noaa.gov/psd/data/climateindices/ list/) as a measure of ENSO influence on TWSAs in South America (Fig. 2). We chose Niño 4 because in many regions of South America GRACE TWSAs are correlated more significantly with central Pacific SST anomalies (represented by the Niño 4 index) than eastern Pacific SST anomalies (represented by the Niño 3 index) over the 10 years of our study, as shown by de Linage et al. (2013). During this period, central Pacific events were more frequent than eastern Pacific events, which might explain the stronger correlation with Niño 4, compared to Niño 3 or Niño 3.4. Analysis by Kao and Yu (2009) indicates different ENSO types have different spatial-temporal dynamics and teleconnections with precipitation.
To represent the influence of tropical North Atlantic SST anomalies on water storage changes in our study domain, we used the Tropical North Atlantic Index (TNAI), also from NOAA's Climate Prediction Center. TNAI is the areaweighted SST anomaly (after removing a long-term monthly climatology) averaged across the region indicated in Fig. 2 (Enfield et al., 1999). TNAI is the oceanic northern component of a meridional sea surface temperature gradient found in the coupled atmosphere-ocean system in the tropical Atlantic Ocean called the Atlantic Meridional Mode. Stronger convergence of surface winds and more precipitation occurs in the warmer hemisphere due to a meridional shift of the ITCZ (Enfield, 1996;Chiang and Vimont, 2004). When the tropical North Atlantic Ocean becomes warmer, the ITCZ shifts northward and northeasterly trade winds weaken along with the northwesterly low-level jet that usually transports moisture to the southernmost Amazon, thus cutting off an im-portant moisture source to the region Zeng et al., 2008). Although the Atlantic Meridional Mode is equivalently well correlated to GRACE TWSAs compared to TNAI, we preferred to use TNAI as a predictor because it remains significantly correlated with TWSAs when the lead times between TNAI and TWSAs are extended (de Linage et al., 2013). For each climate index, the mean over the GRACE period was subtracted from the data.

Model description
We developed three linear models in which TWSA is a solution of the following first-order differential equation: The first term in the right-hand side of the equation represents the evolution of the system in the absence of external forcing (F (x, t) = 0), in other words, the relaxation from the TWSA initial value to zero with a time constant τ , one of the model parameters. We included this relaxation term to reflect the tendency of the plant-soil-aquifer system to retain water during periods of terrestrial water storage deficit as a consequence of negative feedbacks on evapotranspiration and subsurface runoff (and a reversal of this tendency during periods with above-average terrestrial water storage). The second term, F (x, t), represents the forcing (i.e., the impact of the SST anomalies on the balance between precipitation and evapotranspiration at each location, x, within the study domain), which prevents the system from returning to a steady state. For the first two models we considered, the forcing consisted in a single climate index (Niño 4 or TNAI), while that of the third model was a linear combination of both indices (Table 1). We assumed the two indices were independent enough to be used simultaneously as predictors. A lead-lag correlation analysis revealed that Niño 4 explained at most 17 % of TNAI's variance (p < 0.01) from 1950 to 2012, when Niño 4 led TNAI by 4 months. We also assumed that existing teleconnections between the climate indices and TWSAs were stationary, that is, we assumed the parameters a 0 and b 0 in the equations in Table 1 were constant over the calibration and validation periods considered here. In the context of developing a forecasting system, it is important to note that these parameters may change as a consequence of decadal variability in the Pacific (Kosaka and Xie, 2013), including changes in the relative frequency or intensity of central Pacific and eastern Pacific El Niño (e.g., Lee and McPhaden, 2010). Model calibration was done using the 10-year-long time series of monthly GRACE TWSAs. For each model and each parameter set, the solution of the model equation was found numerically using an ordinary differential equation solver, with TWSA (t 0 ) as the initial condition. Among the different parameter combinations that we tested, the one leading to the lowest root mean square error (RMSE) between the observations and the model was selected. We ensured that the minimum RMSE was reached by investigating large ranges for parameters a 0 and b 0 (Table S1 in the Supplement). The range of the relaxation constant τ and that of the predictors' lead times α and β were limited by a priori constraints: an upper limit of 12 months was chosen for τ , while α and β were allowed to vary between 3 and 8 months in our primary set of model analyses presented in Sects. 4.1-4.3. A minimum value of 3 months for α and β was chosen to place a reasonable minimum bound on model prediction lead times. Shorter lead times improved model performance, as shown in Sect. 4.4, but time delays associated with SST processing and data availability make it unlikely this information could be used effectively for operational forecasts. A maximum value of 8 months was selected because it allowed for the investigation of a wide range of possible lead times (i.e., over oneto-two seasons) and because model performance was found to be considerably degraded at longer intervals. Model 3 (the combined model) used lead time values that had been optimized for models 1 and 2 (for the individual Niño 4 and TNAI models; Fig. S1 in the Supplement) to reduce the number of iterations during parameter optimization.

Model evaluation
We computed the coefficient of determination R 2 and the normalized root mean square error (RMSE) of the linear regression to evaluate model performance in each grid cell. Normalized RMSE was computed as RMSE divided by the standard deviation of GRACE TWSAs in each grid cell. For each model and in each grid cell, we performed an F test (null hypothesis F = 0, p < 0.05) to check whether or not the model was statistically better than the temporal mean of the observations. The statistics of the regression are provided in Table 2. Finally, to estimate whether the bivariate (combined) model (forced by Niño 4 and TNAI) was significantly better than each of the two univariate models (forced either by Niño 4 or TNAI), we used another F test (null hypothesis F = 0, p < 0.05) on the residual sum of squares differ-ence, accounting for the difference in the models' degrees of freedom (given in Table 1).

Initial condition and forecasting uncertainties
We estimated the sensitivity of our model predictions to uncertainties in the GRACE TWSA initial condition by propagating the estimated GRACE data uncertainty (±30 mm) in Eq. (1). For each grid cell, we defined the model uncertainty as the mean monthly standard deviation of a set of model estimates in which random errors with a Gaussian distribution were added to the initial condition.
As the available GRACE record is relatively short, we could not both train and evaluate our model with equally long data records. However, to provide an estimate of forecast uncertainties, we computed an ensemble of 10 models, each of which was trained using monthly observations from 9 out of the 10 available years, following the least-square approach described above for the model calibrated with the full 10year time series. Models developed using the 9 years of data were then used to estimate TWSAs for the left-out year. We then computed a composite 10-year time series comprised in each year of the simulation from the ensemble member that was trained (calibrated) over all years but the left-out year. The total model RMSE (or NRMSE) was defined from difference between this composite time series and the observed GRACE time series. We then estimated the forecast uncertainty as the difference between the total model RMSE (or NRMSE) and the RMSE (or NRMSE) from the model calibrated with the full 10-year time series. In the results, we first present the model calibrated using the 10-year time series because of its simplicity (its parameters remain constant over the study period). We then explicitly compare the different sources of model uncertainty, including the forecast uncertainty derived from the composite time series, in Sect. 4.5.

Spatial patterns of climate index influence on terrestrial water storage anomalies
The influence of Niño 4 and TNAI on terrestrial water storage interannual variability varied considerably across different regions in tropical South America (Fig. 3). Niño 4 explained the largest amount of TWSA variability in the northeastern Amazon River basin (accounting for 61 % of the variance in region B as shown in Table 2). In contrast, TNAI was the primary external forcing across central and western regions of the basin (accounting for 35 % of the variance in region A). In northern South America, including the Orinoco River basin, TNAI explained more than half of the observed variance (54 % in region C). Combining Niño 4 and TNAI forcing terms significantly (p < 0.05) improved model performance for many regions within South America (Tables 2 and 3   . Coefficient of determination (R 2 ) and normalized root mean square error (NRMSE) for the three models described in Table 1 and computed for the period from August 2002 through July 2012. The RMSE was normalized at each location using the standard deviation of GRACE interannual terrestrial water storage anomalies shown in Fig. 1b. The prescribed minimum lead time for the three models was 3 months. Cells where the linear fit was significant (p < 0.05) are marked with a black dot. Supplement), with 39, 66 and 63 % of variance explained in regions A, B and C, respectively. In each region, the combined model significantly outperformed (F test, p < 0.05) the best of the two univariate models in at least 64 % of the cells (Table 3). The improvement from the addition of Niño 4 was higher and more widespread in northern South America than in the central and western Amazon, and more significant than the addition of TNAI in the northeastern Amazon. Given the 3-month minimum lead time between the climate indices and model predictions, the relatively high amount of variance explained by the combined model in many regions indicates considerable potential to develop statistical models for TWSA forecasting. Depending on the sign and magnitude of the secondary climate index, model estimates of positive or negative TWSAs were either enhanced or damped in specific regions. Dry con-ditions in the central and western Amazon were triggered by warm SST anomalies either in the tropical North Atlantic or in the central Pacific (Niño 4 and TNAI were in-phase) ( Fig. S2a and b in the Supplement). In contrast, in the northeastern Amazon and in northern South America, dry conditions were caused by either anomalously warm central Pacific or cold tropical North Atlantic SSTs (Niño 4 and TNAI were in anti-phase). These relationships are consistent with previous studies documenting relationships between SST anomalies and precipitation (Ropelewski and Halpert, 1987;Zeng, 1999;Liebmann and Marengo, 2001;Ronchail et al., 2002;Gloor et al., 2013;Grimm, 2003;Molinier et al., 2009;Yoon and Zeng, 2010;Espinoza et al., 2011) or SSTs and streamflow (Richey et al., 1989;Dettinger et al., 2000;Dai et al., 2009;Molinier et al., 2009;Espinoza et al., 2011). We refer to de Linage et al. (2013) Fig. 1a, for the three models described in Table 1. Regional averages of model performances (Table 2) are also provided for each model. The prescribed minimum lead time between the climate indices and TWSAs was equal to 3 months.
the intrinsic lags between climate indices and precipitation, TWSA, and streamflow observations. The 2010 drought had a considerable impact on TWSAs in every region and was preceded by very large positive Niño 4 and TNAI indices (Fig. 4). The 2005 drought was visible in the central and western Amazon, where it was responsible for a wide trough in the time series. The 2005 drought also was associated with positive Niño 4 and TNAI indices, although the magnitudes were smaller as compared to 2010. The impact of the 2011-2012 La Niña on TWSAs may have been damped because of the considerable water deficits remaining after the 2010 drought.

Model lead times and relaxation parameters
The optimal lead time for a given climate index was usually shorter in regions where TWSAs were dominantly forced by this index (Fig. S1a in the Supplement). Lead times for Niño 4 were 3 months in most grid cells, except for several small contiguous areas in the central and western part of the basin (where they ranged up to approximately 6 months). Lead times for TNAI were uniformly 3 months in the central and western Amazon, and were much longer elsewhere (from 5 months in region C to 7-8 months in region B).
The shortest lead times obtained from the optimization were often equal to the prescribed minimum value and generally decreased when the prescribed minimum lead time was reduced to zero (Fig. S1b in the Supplement), indicating that in many areas there was no intrinsic phase difference between the timing of the primary index and its impacts on TWSAs for the models considered here. Only areas near the Amazon River main stem and river mouth retained an irreducible 3-month lead time for Niño 4, which was consistent with expected time delays introduced by routing of runoff through Amazon sub-basins, and the importance of surface water in this region driving month-to-month variations in terrestrial water storage. The impact of changing the minimum lead time on models skill is discussed in Sect. 4.4.
A shorter relaxation time τ indicates a shorter memory and thus a stronger sensitivity of the system to rapid changes of the climate indices. In either the central and western Amazon or the northeastern Amazon, mean relaxation times were shorter for the model forced by the dominant index, while they were longer for the secondary index ( Fig. S2c in the Supplement). In northern South America, however, relaxation times were larger for the dominant index, TNAI, than for Niño 4. For the combined model, the mean relaxation time in the central and western Amazon and the northeastern Amazon was approximately one-month longer than those obtained from the univariate model forced by the dominant index. Larger τ values (6/12 months) were found in the downstream parts of the main rivers, floodplains and wetlands, where the land surface memory increases due to longer residence times of surface water storage and time delays associated with stream and river transport. These parameter values are broadly consistent with results from Bonnet et al. (2008), who estimated a 3 month residence time for the water within a floodplain lake, and an estimated 1 month delay for every 900 km of water transport through the river and floodplain network based on an effective routing velocity (Miller et al., 1994).  Fig. 1a, for the three models described in Table 1 (with a 3-month prescribed minimum lead time). Dashed lines represent the R 2 computed for all months. The gray shaded areas represent the 3 months with the lowest average TWSA, computed from the GRACE TWSA mean monthly climatology.

Seasonal variations in model performance
Each climate index had considerable month-to-month differences in their standard deviations, with annual peaks occurring between December and February for Niño 4 and March through May for TNAI (Fig. 5). Likely as a consequence of the different seasonal dynamics of the climate indices and their teleconnections with TWSAs, model performance varied considerably over an annual cycle (Fig. 5). Seasonal variation in the R 2 metric for the combined model was larger for the central and western Amazon than for other regions. Within each region, the monthly R 2 of the best two models had a similar phase. In the central and western Amazon, the combined model had higher levels of performance between May and October, before and during the time terrestrial water storage was near its annual minimum (i.e., during the precipitation dry season, see Ronchail et al., 2002). In the northeastern Amazon, the combined model performance was highest from January to June. This period spans the time when water storage increases in floodplain and wetland areas along the Amazon main stem downstream of Manaus (Prigent et al., 2007;Paiva et al., 2013). In northern South America, we found that the highest R 2 values occurred between March and June, during and after the annual minimum in terrestrial water storage (corresponding to the end of the dry season and beginning of the wet season for precipitation).

Changes in model performance with varying lead times
To study the degradation of model performance with increasing minimum lead times, in another set of simulations we varied the lower bound of the climate index lead times from 1 to 8 months. Overall, model performance degraded monotonically (R 2 decreased and RMSE increased) with increasing minimum lead times. In most areas, the model trajectories usually did not cross each other (Figs. 6 and S4 in the Supplement), so that the model ranking for any given lead time was the same as that found for the nominal 3-month minimum lead time described above. In the central and western Amazon, however, model 1 (Niño 4) provided better extended forecasts than model 2 (TNAI), while the opposite was true for medium-and short-range forecasts. Also, the combined model was not significantly (p > 0.05) better than model 1 for extended forecasts with 8-month lead times.
In the central and western Amazon the evolution of the model skill was mainly influenced by TNAI, while in the northeastern Amazon, it was influenced by Niño 4, so that the trajectories of the combined model and the best individual univariate model tended to be parallel in these regions. In northern South America, both indices had important, complementary roles as predictors: Niño 4 influenced the short-to medium-range forecasts, while TNAI influenced the longer-range forecasts. As a result, in this region the combined model diverged from the univariate model forced by TNAI for short and medium lead times. The time-varying Hydrol. Earth Syst. Sci., 18, 2089-2102, 2014 www.hydrol-earth-syst-sci.net/18/2089/2014/ Table 4. Initial condition and forecasting uncertainties averaged over each region delineated in Fig. 1.
In the Amazon Basin, the combined model was able to explain 31 % of the observed variance for an 8-month lead time compared to 45 % for a 1-month lead time (Fig. S4 in the Supplement). Model degradation with increasing lead times was smaller in the northeastern Amazon (where 49 % of the variance was explained at 8 months versus 69 % at 1 month) than in the central and western Amazon (where 22 % of the variance was explained at 8 months versus 43 % at 1 month). In northern South America, model degradation was even smaller, with still 52 % of the variance explained at 8 months versus 65 % at 1 month.

Initial condition and forecasting uncertainties
Normalized monthly mean errors arising from uncertainties in initial conditions in the different regions ranged between 2 and 10 % ( Table 4). As shown in Fig. S5 in the Supplement, the errors introduced from the initial conditions decreased rapidly within the first year. The spatial pattern of normalized model error for initial condition uncertainty is shown in Fig. 7a. The forecasting error, defined in Sect. 3.3, represents the additional error introduced by forecasting TWSAs using a model developed for a different time period. These errors had NRMSEs ranging from 3 to 12 % (or RMSEs from 3 to 11 mm), depending on the region and model, with the largest uncertainties found for the combined model (forced by both Niño 4 and TNAI) (Fig. 7b, Table 4). The forecast error associated with the combined model was approximately 10 % of the total error (Fig. 3, Table 2) in the Amazon Basin and up to 20 % in the northern South America. As expected, the forecasting error increased with the prescribed minimum lead time (Fig. 8).  Fig. 4). The total error is the sum of the forecasting errors (b) and the prediction error (Fig. 3).

Potential for flood and drought forecasting in the Amazon
Floods hinder ground transportation and damage crops and infrastructure in many regions in South America. The northeastern Amazon is characterized by an extensive network of wetlands and floodplains (Prigent et al., 2007;Paiva et al., 2013) that are inundated seasonally. A strong wetting trend was observed by GRACE in these regions (Fig. 4), and may be related to an equatorial Pacific surface cooling trend (de Linage et al., 2013;Kosaka and Xie, 2013) or with warming of the tropical North Atlantic over the last 20 years (Gloor et al., 2013). Both of these pathways are consistent with the sign of our coefficients linking Niño 4 and TNAI with predicted TWSAs (Fig. S2 in the Supplement). Large hydroelectric reservoirs are also found in these regions (like the Balbina Reservoir near Manaus) that are used to generate hydroelectricity. Our 3-month combined model explained a considerable amount of the variance in this region (66 % in region B) and had a level of performance higher than average during the precipitation wet season (i.e., when floodplains and wetlands begin to fill with excess water runoff). An important next step in this context is to relate the TWSA observations and model predictions analyzed here to floodplain area, river height, and reservoir height time series using additional satellite observations (see Paiva et al., 2013). Although floods have significant economic impacts, farmers and fishermen may be more vulnerable to droughts, because droughts reduce crop yields, hinder river transportation, and harm fish communities . Also,  Fig. 1a. The difference between the two RMSEs is the forecasting error (shaded area). The vertical dotted line represents the minimum lead time that was used to display the forecasting error spatial distribution shown in Fig. 7b. droughts impact larger areas than floods: the 2005 and 2010 droughts respectively affected 1.9 and 3 × 10 6 km 2 of land, based on satellite precipitation data (Lewis et al., 2011). Severe hydrological droughts like the 2005 and 2010 events may increase tree mortality in intact forests and promote understory fires that contribute to the conversion of tropical forests to savannas and croplands. Regions that are critically affected by fires include the Amazon's central and southern regions (the so-called "arc of deforestation", see Chen et al., 2013b). In these regions, the amount of soil moisture recharge during the wet season critically affects transpiration, surface humidity, fire weather, and fire spread rates during the following dry season (Chen et al., 2013a). In the central and western Amazon our 3-month model predictions had above-average seasonal performance during the precipitation dry season (May-August), suggesting that TWSA may be another useful variable for integration into an early warning system for Amazon fires.

Ocean-atmosphere dynamics and implications for TWSA prediction
Complex teleconnections exist between tropical Pacific and North Atlantic ocean-atmosphere systems. Linkages between the eastern Pacific and the tropical North Atlantic are well known, with El Niño events leading to positive Atlantic SST anomalies with delays of approximately 2-4 months (Enfield, 1996;Giannini et al., 2000;Saravanan and Chang, 2000). This is consistent with TNAI lagging Niño 4 by 4 months (and both indices being positively correlated) and explains the severity of the 2005 and 2010 droughts when both pools were warmer than usual (Fig. 4). On the other hand, positive SST anomalies in the tropical North Atlantic contribute to La Niña conditions in the central Pacific three seasons later (Ham et al., 2013), which may partly explain the 2006 and 2011 La Niñas. In terms of the statistical models evaluated here, it suggests that the two forcing terms in the combined model are not completely independent, which may lead to a slight ambiguity in identifying the forcing source of the observed TWSA variations. In future work, spatial optimization of the forcing regions within the Pacific and Atlantic by means of empirical orthogonal function analysis  or other multivariate methods (Westra et al., 2008) may lead to increases in model skill. Regional indicators of the atmospheric circulation (Yin et al., 2014) also may be useful in model refinement. Combined, these steps may reduce possible correlations among the predictors and improve our understanding of the modes that drive hydrological cycle variability in South America.
In some regions of South America, our analysis provides evidence that predictability of TWSAs may be possible for relatively long lead times. For example, for the northeastern part of the Amazon, performance of the combined model was reduced by only 30 % when lead times were extended from 3 months (R 2 = 0.66) to 8 months (R 2 = 0.49). In this region, Niño 4 was the dominant index and some of the success of the model for the longer lead times may be attributable to intra-seasonal persistence (and characteristic seasonal behavior) of the climate index. During the building phase for ENSO (from JJA to DJF), for example, SST anomalies vary monotonically toward their peak value, sustained by positive feedbacks between surface winds and SSTs (Rasmusson and Carpenter, 1982).

Importance of regional and local controls on TWS variability
The underlying teleconnections between climate variability in the ocean regions used to drive our model and Amazon hydrological processes are highly non-linear and involve complex interactions with regional atmospheric circulation and the local land surface. During austral summer, the South American monsoon (SAM) serves as a likely vector for these teleconnections (Zhou and Lau, 1998). However, depending on the timing of ocean-atmosphere variability in Pacific and Atlantic regions with respect to the onset of the SAM, either forced or internal variability dominates in the atmosphere over South America, inducing confusion on the causality of the observed precipitation anomalies (e.g., Yin et al., 2014). The SAM starts in austral spring with a low-level convergence over the Amazon and stronger trade winds (easterlies) flowing from the tropical North Atlantic (Sahara high) toward the northern Andes. Since the SAM onset is concurrent Hydrol. Earth Syst. Sci., 18, 2089-2102, 2014 www.hydrol-earth-syst-sci.net/18/2089/2014/ with the ENSO building phase, ENSO is likely to interact with the initiation of the SAM (its tropical branch) and may impact its propagation. A higher model skill was indeed found at this time, in other words, during the local wet season in the Amazon northeastern regions (Fig. 5). The SAM mature phase is characterized by convection and rainfall over its subtropical branch in austral summer and is called the South American Convergence Zone or SACZ (Nogués-Paeble and Mo, 1997). The SACZ extent is modulated by ENSO on interannual timescales (Nogués-Paeble and Mo, 1997;Carvahlo et al., 2004). During El Niño the SACZ is shifted offshore in the western South Atlantic, hence suppressing precipitation over eastern tropical South America. Precipitation is concurrently enhanced in subtropical South America due to an increased meridional moisture transport from tropical regions to subtropical regions. Future work should investigate the feasibility of developing statistical TWSA prediction models for these subtropical regions.
As part of the SAM, the equatorial northeasterly winds are re-directed by the northern Andes and become northwesterlies traveling along the eastern flank of the Andes to the southeast (e.g., Zhou and Lau, 1998). Thus the western and southern Amazon regions provide moisture to the downwind subtropical regions, and may be prone to droughts whenever their own moisture source (provided by the equatorial trades from the Atlantic) is reduced because of El Niño or tropical North Atlantic warm events. These complex dynamics may partly explain why our model skill was lower in central, western, and southern regions of our domain as compared to the northeastern region, which is more directly under the influence of ENSO.
The impact of interannual SST anomalies in the tropical North Atlantic on the SAM is less studied relative to that of ENSO. Warmer tropical North Atlantic SSTs weaken the equatorial trade winds and shift the ITCZ to the north, which in turn results in the weakening of the northwesterly lowlevel jet that usually brings moisture from the north of the continent to the southern Amazon regions, and reduces convection over these regions . The maximum of TNAI variability occurs in austral fall, at the end of the SAM, so it may not directly affect the onset of the SAM. Further work is needed to better understand how the tropical North Atlantic circulation interacts with the regional dynamics in tropical South America.

Model structure and uncertainties
Other studies (Chen et al., 2013a;de Linage et al., 2013) have used climate indices as a proxy for TWSAs or other observables. Our approach was different in essence because (1) we considered predictive models that can inform the development of seasonal forecasts and (2) we modeled TWSA changes (dTWSA/dt) as a function of SST anomalies, or equivalently, that TWSAs responded to the integral of SST anomalies. This form of the equation is consistent with SST anomalies imparting precipitation-evapotranspiration imbalances in South America with relatively short atmospheric transport times by means of the ocean atmosphere mechanisms described above in Sect. 5.3. We also modeled TWSA as experiencing a forced relaxation (as described by Eq. 1), thus accounting for the land surface memory. In the absence of any additional external forcing (P -E anomalies), the relaxation term causes the system to gradually return back to the mean hydrologic state. During times with positive TWSAs (periods of land water excess), this return to the mean would likely occur by means of increased drainage, runoff, and evapotranspiration. During times with negative TWSAs (periods of water deficit), return to the mean would occur by means of reduced plant transpiration, reduced soil evaporation, higher soil water retention (lower conductivity), and reduced surface and subsurface runoff. We used a spatially varying parameterization because the relaxation time is likely to change from region to region as a function of storage capacity and routing. For example, longer relaxation times were expected (and found) in regions near the mouth of Amazon River, where transport times from upper basins introduce time delays and variations in floodplain storage are relatively large. Important next steps for reducing model uncertainties include (1) a more comprehensive evaluation of model forecasting success for time periods that were not used to calibrate model parameters as more GRACE observations become available, (2) improvements to the approach for optimizing parameters for the combined model, and (3) more detailed analysis and study of the mechanisms contributing to TWSA interannual variations in different South American regions. Forecast errors obtained using our approach of parameterizing the model using 9 years of observations and predicting the remaining (left-out year) likely underestimate the true range of uncertainty, given the relatively limited range of Pacific and North Atlantic Ocean SST variability observed over the past decade. Improved parameterizations may include the development of new models that allow for a seasonally varying TWSA sensitivity to the forcing, to increase the inter-seasonal performance of our model. To improve our understanding of the underlying mechanisms, more study is needed of the relative contributions of remote oceanic sources to local and regional precipitation patterns (van der Ent and Savenije, 2013), as well as more detailed analysis of vegetation and deforestation influence on precipitation recycling by means of land-atmosphere couplings (Coe et al., 2009;Lee et al., 2011).

Conclusions
We found that the spatial pattern of interannual TWSAs in tropical South America was significantly influenced by variations in SST anomalies from the tropical Pacific and North Atlantic, as represented by Niño 4 and TNAI, respectively. To predict the spatial and temporal variability of TWSA in tropical South America, we built a series of simple statistical models. Niño 4 was the primary external forcing for the northeastern Amazon (with 61 % of variance explained with a 3-month forecast), whereas TNAI was dominant in the central and western regions of the southern Amazon (35 % of variance explained with a 3-month forecast). Forcing the model with a combination of the two indices improved the fit significantly (p < 0.05) for more than 64 % of the grid cells, compared to models forced solely with Niño 4 or TNAI. The combined model was able to explain 43 % of the variance in the Amazon Basin as a whole, 66 % in the northeastern region, and 39 % in the central and western region.
We studied how the predictive skill of our combined model was degraded with increasing lead times from 1 to 8 months. For the Amazon Basin as a whole, the model was still able to explain 31 % of the observed variance using an 8-month lead time versus 45 % for a 1 month lead time (equivalent to a 31 % degradation). Model degradation with increasing lead times was smaller in the northeastern Amazon (up to 29 %) and larger in the central and western Amazon (up to 49 %).
These statistical models have the potential to provide early warning information about flooding in the northeastern Amazon, where floodplain areas are extensive and the sensitivity of floods to external SST forcing was high (model skill was up to 25 % higher than the annual mean during the wet season and the time of peak of flooding in this area). They also may enable drought and fire season severity prediction, in particular in the central and western Amazon where our models had above-average performance (up to 50 % higher) during the dry season.
An important next step is to improve our understanding of ocean-atmosphere processes that enable these statistical predictions, including for example, how the South American monsoon and internal atmospheric variability influence the spatial pattern of TWSA variability over the continent. Further work also is needed to validate the models developed here against longer time series of GRACE observations and to relate predictions of TWSAs to estimates of flood and drought extent. Refinement of the statistical models may be possible in the future by allowing parameters to vary seasonally and by optimizing the ocean regions used as model drivers.
The Supplement related to this article is available online at doi:10.5194/hess-18-2089-2014-supplement.