A millennial multi-proxy reconstruction of summer PDSI for Southern South America

. We present the ﬁrst spatially explicit ﬁeld reconstruction of the summer (DJF) Palmer Drought Severity Index (PDSI) for the Southern Hemisphere. Our multi-proxy reconstruction focuses on Southern South America (SSA, south of 20 ◦ S) and is based on a novel spectral analogue method that aims at reconstructing low PDSI frequencies independently from higher frequencies. The analysis of past regimes and trends in extreme wet spells and droughts reveals considerable geographical and temporal variations over the last millennium in SSA. Although recent changes are in some cases notorious, most were not exceptional at the scale of the last thousand years. Our reconstruction highlights that low frequency water availability ﬂuctuations in Patag-onia were generally in antiphase with the rest of the subcontinent. Providing the fact that modern patterns of changes are transferable to the past, we show that such antiphases within SSA’s hydroclimate could be attributed to the spatially contrasted response of summer PDSI to the Antarctic Oscillation (AAO). However, El Ni˜no Southern Oscillation (ENSO) and Paciﬁc Decadal Oscillation (PDO) signals are also embedded within the PDSI series during the 20th century. All these ocean-atmospheric forcings acted syn-ergically, but the dominant inﬂuence appeared highly compartmentalized through space, highlighting clear AAO-(e


Introduction
Numerous evidences exist to support an increasing trend in the severity and intensity of both extremely dry and wet spells over the last century in Southern South America (SSA) (Trenberth et al., 2007;Magrin et al., 2007). These episodes have among the costliest impacts on the economy, society, and natural environment in the area. However, the analysis of present-day climate has shown that an important geographical variability also exists regarding these trends. Indeed, some parts of SSA like northeastern Argentina, southern Brazil, Paraguay, and Uruguay are actually experiencing a wetter climate by comparison to the conditions that prevailed at the beginning of the 20th century (Dai et al., 2004;Magrin et al., 2007). By contrast, southern Chile, south-west Argentina, and southern parts of Peru recorded a clear decline in the amount of precipitation. A general rise in summer temperatures is also superimposed in most regions, except in central Argentina where a cooling of maximum temperatures has been observed during the second part of the 20th century (Magrin et al., 2007). These contrasted climatological trends all impact water resources. Moreover, a large part of the variability associated with these changes originates from the influence of tropical (El Niño Southern Oscillation, ENSO) and high-latitude (Antarctic Oscillation, AAO) ocean-atmospheric climate forcings over the area (Garreaud, 2007). From a geographical point of view, the influence of these forcings is modulated by the presence of the Andes that create strong W-E climatic gradients over the continent, resulting in considerable climatological inhomogeneities through time and space in SSA. Given the number of factors involved and the heterogeneous spatio-temporal trends, it becomes crucial to place the recent hydro-climatic fluctuations in a broader perspective in order to better identify the long term and large scale patterns associated with these changes.
Published by Copernicus Publications on behalf of the European Geosciences Union.
958É. Boucher et al.: A millennial multi-proxy reconstruction of summer PDSI Whereas instrumental records are usually too short and scarce to perform such an analysis, highly resolved multiproxy reconstructions can be useful to place recent climatic trends in a larger spatio-temporal perspective. Most notably, they allow comparison of recent trends with well-known climatic periods such as the Medieval Warm Period (MWP) or the Little Ice Age (LIA) and placement of local patterns of change into a continental, hemispheric or even global perspective. Yet, most efforts have focused on reconstructing past temperatures in the Northern Hemisphere (Bradley and Jones, 1993;Jones et al., 1998;Mann et al., 1998Mann et al., , 1999Briffa et al., 2002;Moberg et al., 2005), owing to the number of available proxies and the density of high quality instrumental data. While local reconstructions have existed since the 1970s in SSA (see Boninsegna et al., 2009, for a review), spatial reconstructions are relatively new. The first highly resolved multi-proxy climate field reconstructions were performed recently (Neukom et al., 2010a,b). The summer temperature reconstruction showed that the MWP (9th to 14th century) was generally warmer than the last century's average. Complementarily, the precipitation reconstruction covered the last 500 yr and showed that summer and winter precipitation behaved in opposite directions over the last centuries (summer precipitations increased since the LIA, while winter precipitations decreased). However, despite these findings, it remains clear that the continental hydrological balance itself (i.e. the linkages between water inputs, storage, and outputs) is affected by all these parameters synergically. Thus, the reconstruction of an individual parameter (either temperature or precipitation) does not yield a full representation of past water availability fluctuations. Moreover, moisture conditions depend not only on precipitation and temperatures occurring during the season of interest, but also on previous moisture conditions. Thus, to better describe water availability fluctuations, one has to take into account both water supply and demand at the earth's surface and evaluate the changes considering prior conditions. The Palmer Drought Severity Index (PDSI, Palmer, 1965) is a prominent index that incorporates antecedent moisture conditions and has been extensively used to monitor past droughts and wet spells in the US and elsewhere (Apipattanavis et al., 2009). It has also commonly been reconstructed from dendrochronology mostly in North-America (Cook et al., 1999(Cook et al., , 2004Woodhouse and Overpeck, 1998;Woodhouse and Brown, 2001) and in the Mediterranean basin (Nicault et al., 2008). In SSA, Christie et al. (2009) have successfully reconstructed high frequency variations of the PDSI in the Temperate Mediterranean Transition in the Andes (TMT) back to the 14th century. However, as for temperature and precipitation, the PDSI also needs to be reconstructed at larger spatial (i.e. sub-continental) and temporal (i.e. millennial) scales in order to place the recent changes in water availability into a broader perspective. Our work therefore focuses on reconstructing austral summer (DJF) PDSI in SSA from a multi-proxy approach. We particularly aim at reconstructing the full spectra of PDSI variations and at analyzing the possible role of ocean-atmosphere forcings (AAO, ENSO, PDO) on large-scale and long-term PDSI variations in SSA.

PDSI
Our aim is to produce a millennial gridded reconstruction of the Palmer Drought Severity Index (PDSI, Palmer, 1965) in SSA (south of 20 • S, the northern limit of most tree ring records in the area). The PDSI is a climatic metric that measures the departure from abnormally dry or wet conditions (relative to average local conditions). Three input variables are required to calculate the PDSI: monthly precipitation, monthly temperatures, and the available water content of soils. Within the calculation of the PDSI, moisture supply and demand are approximated in a simple hydrological accounting model that evaluates evapotranspiration, soil recharge, runoff, and moisture losses from the surface layer (see Palmer, 1965 for original equations). The PDSI produces monthly indexes of meteorological droughts that account for both contemporary and antecedent climatic conditions. The PDSI fluctuates between −10 (extremely dry) and +10 (extremely wet), with 0 corresponding to the local average conditions. The PDSI should be regarded as an approximation for dryness or wetness conditions rather than as an absolute value, as many factors like solar irradiation or aerosol concentrations are not considered explicitly. In the present paper, we use the gridded (2.5 • × 2.5 • ) PDSI dataset computed by Dai et al. (2004) (Fig. 1a) and extended it to the last millennium in SSA. This dataset covers the 1870-2004 period and was calculated from long temperature (CRUTEM2, Jones and Moberg, 2003) and precipitation (Dai et al., 1997;Chen et al., 2002) series. We provide reconstructions of summer (DJF) PDSI for each pixel in SSA. Summer PDSI better integrates the effects of summer temperatures and precipitation along with lagged effects from previous seasons. We finally provide regional reconstructions for four regions similar to those identified by Neukom et al. (2010a) (Fig. 1a): Patagonia (PG), the Pampas (PM), Sub-Tropical SSA (ST), and the Andes (ANDES). Although climatic conditions and trends are not always uniform within these regions (e.g. PG, ANDES), the coarse source grid used here does not allow subdivision into smaller subsets.

Proxies
To achieve the PDSI reconstructions, we adopted a frequency-dependent approach similar to the one used by Moberg et al. (2005) where each proxy is used only to reconstruct the periodicities for which it provides reliable information (see Sect was compiled from high-and low-resolution proxies (Table 1, Fig. 1b,c). High-resolution proxies are mostly tree rings, most of which are available on the International Tree Ring Data Bank (ITRDB). We kept only the mean tree ring chronologies that were longer than 250 yr, giving a total of 82 tree ring predictors. In order to circumvent the well known problems associated with the standardization of tree ring series and the problematic interpretation of medium to low frequency variations (Esper et al., 2002(Esper et al., , 2004, we retained only the high frequency variations contained within these series. Low frequencies were removed in each tree ring chronology using a low-pass filter (based on the Fast Fourier Transform procedure). Basically, we removed all periodicities below a 12-yr threshold (f ≤ 0.08). In order to avoid the over-representation of tree ring series in the multi-proxy reconstruction, we selected (after infilling with the AM method, see next section) the first 7 PC explaining about 30 % of the variance of the high frequency tree ring series. Low-resolution proxies included marine sediments, lake sediments, and ice cores from the Andes and Antarctica (Table 1). The Frías sequence (Chapron, unpublished) is a new sedimentary record from a varved proglacial lake (Ariztegui et al., 2007) draining the Frías River Valley. The Frías sequence consists of a detailed record of goethite content measured every 0.5 cm by sediment diffuse reflectance on a 187 cm long core retrieved in 2008 within Lago Frías main basin. The chronology (Fig. A1) has been established back to AD 1649 by varve counting on digital core images and by the identification of three mass wasting deposits induced by historical large earthquakes in AD 1960; AD 1751, and AD 1737. Geothite content in these annually laminated fine-grained sediments reflects soil erosion after rainfall.

Analogue method (AM)
The analogue method (Guiot et al., 2005) is commonly used to infill proxy-climate matrix containing missing values. The AM procedure aims at identifying, for each year i in the past where no PDSI value exists, the year k within the instrumental record that has the most "similar" proxy vector. Here, similarity (d 2 ik ) is measured as an Euclidian distance (Guiot et al., 2005): where x ij represents the value of proxy j at year i, x kj corresponds to the value of proxy j at year k, M and m refer to the maximum and minimum value of the proxy j . δ J ik is an index equal to 1 if a value is available for proxy j at year i and k, and equals 0 otherwise.
In the present paper, the AM was first used to infill both proxy and PDSI matrices ( Fig. 2[1]). Since most available tree rings series terminated before ∼1993 (Table 1), we chose not to infill the proxy matrix up to the present in order to reduce the possible bias associated with the use of an infilling method during the calibration period. Additionally, Dai et al. (2004)'s PDSI dataset presented some irregularities (in both mean and variance) before ∼1930 in SSA. These irregularities are most likely associated with the very low number and poor quality of temperature and precipitation records in the early 20th century in that area. We therefore defined the 1930-1993 period for calibration. A total of 101 PDSI series were retrieved from Dai et al. (2004)'s dataset. We selected the first 12 PC explaining 77 % of the variance of the original PDSI series (Fig. 2[2]).

Spectral analogue method (SAM)
The reconstruction of the 12 PC of PDSI is also based on the AM, but it is performed separately for each frequency band using a spectral analogue method (SAM). The SAM (Guiot et al., 2010) is a combination of the AM with a spectral decomposition procedure that aims at achieving the reconstruction separately for each frequency band (HF: high frequency band, LF: low frequency band). We defined the LF band as all periodicities below f = 0.08 (or T <12 yr, fixed by experimentation) and the HF band as all periodicities comprised between 0.08 and 0.5 (12 yr > T > 2 yr). While it is assumed that all frequency bands are present in the PDSI series (0 < f < 0.5), the SAM method allowed choosing, for each frequency band, which proxies are the most reliable predictors. The decomposition of both proxy and PDSI series into their LF and HF component is achieved through a spectral decomposition algorithm (SD, Fig. 2[3]) based on the Fast Fourier Transform (FFT). The FFT first transforms each series into the frequency domain. Then, all frequencies outside the band are set to zero. Finally, the residual spectrum is back-transformed into the time domain using the inverse FFT. LF and HF bands are complementary because their summation restitutes the original series. Thus, for each frequency component, the AM is used to reconstruct the corresponding band of the PDSI PCs over the last millennium ( Fig. 2[4]). The LF and HF bands of each PC are reconstructed independently from one another. Finally, the two bands are restituted into one single band (SR) and back-transformed into the original gridded PDSI data using an inverse PC (PC −1 ) procedure ( Fig. 2[5,6]).

Uncertainty assessment
In the present paper, steps 4 to 6 on Fig. 2 were repeated iteratively (100 times here) and allowed for the implementation of an h-block jackknife procedure (Guiot et al., 2010), which provided independent validation statistics and confidence intervals for the reconstructed values. At each iteration, a randomly selected year t was removed from each frequency band, and a training dataset defined after the exclusion of year t along with the four preceding and following years (h = 9). The AM was then used to predict the left-out year t from the training dataset in each band. This adaptation of the traditional jackknife method has proven to be particularly useful with series (like LF series) that present a high degree of autocorrelation (Guiot et al., 2010). In our case, it forbade that the best analogue for year t could be found between year t − 4 to t+4. An analysis of the autocorrelation structure within summer PDSI series revealed that the years located immediately outside the block had less than 30 % of variance in common with year t. Such an h-block jackknife procedure was performed 100 times with different randomly chosen years and allowed for the calculation of the RE statistic (Reduction of Error); a validation statistic ( Fig. 2[7]) commonly used in dendroclimatology (Fritts, 1976). A check of the autocorrelation structure of residuals was also performed using the Durbin-Watson (DW) statistic (Durbin and Watson, 1950). Moreover, the information pro-duced at each iteration allowed for the computation of 95 % confidence intervals around reconstructed values. As a further validation for our reconstructions, we used a published sedimentological record of Mar Chiquita's lake levels fluctuations (Piovano et al., 2002). Mar Chiquita is a large Sub-Tropical lake (about 6000 km 2 in size) that drains a vast area (about 127 000 km 2 ). It is reasonable to assume that the fluctuations of such a vast hydrosystem are controlled by parameters (temperature and precipitation) that are similar to those that also influence summer PDSI in the area. Thus, this record might be adequate to evaluate the quality  of our reconstruction, at least the long-term trends. Finally, despite the fact that a number of proxies are shared between the studies, we compared our PDSI reconstructions with the precipitation and temperature reconstructions of Neukom et al. (2010aNeukom et al. ( , 2010b. This was done to analyze which climatic parameter (P or T ) is the best related, over time and space, to summer PDSI in SSA.

Calibration and validation statistics of the SAM
The SAM yields satisfying results and suggests that the PDSI can be reconstructed in most areas of SSA. A good fit exists between observations and predictions for the full calibration period   (Fig. 3), and this fit is reasonably good in each frequency band (Fig. 4). The mean correlation coefficient (R) for SSA is 0.64 and ranges from 0.52 (minimum) to 0.78 (maximum) ( Table 2) for the full spectra. The average correlation coefficient is also very similar between the studied regions, with the lowest R value found in PM (0.61) and the highest in the PG region (0.67). On average, the SAM reconstructs about 37 % of the variance of the PDSI in SSA ( Table 2). The percentage of variance accounted for by the model (coefficient of determination, R 2 ) varied between 19 % and 58 %. Figure 5 shows the spatial distribution of these calibration statistics. In general, the relationships between observations and predictions tend to be slightly weaker in the southern part of the Andes. A possible reason could be that the observed PDSI values might be less well estimated in this region of high climatological variability and altitudinal contrasts, a situation that would result in locally spurious trends. A second reason could be that the proxies themselves, although numerous in that area, do not cover the full range of variability within these pixels. amount of the original variance within the series if the reference period is sufficiently diversified. In other words, when a regression technique is used, the proportion of reconstructed variance decreases with the R 2 . This attenuation effect becomes even more exaggerated in the past because the number of proxies, and then the R 2 , diminishes with time. Yet spectral analogues are less sensitive to these effects since they are similarity-based and consequently they are able to reproduce all the variance contained in the reference period. To illustrate this effect in SSA, we calculated the ratio of standard deviations between reconstructed and observed PDSI series in each region and observed that the former corresponds to about 85 to 95 % of the latter. Therefore, we conclude that this property facilitates the analysis of extremes. Validation statistics confirm the predictive skills of our model. RE is generally positive over SSA (Table 2), suggesting that the SAM provides better estimates of the PDSI than the climatology. Moreover, the DW statistic is always very close to 2, an indication that there is no significant first-order autocorrelation structure in the residuals (Durbin and Watson, 1950). The mean RE value was, on average, around 0.30 in all sub-regions except the Andes (0.25) but generally oscillated between 0 and 0.55 (taking into account the 95 % confidence intervals over all SSA). These performances are comparable to the RE values obtained with the same method applied in Europe to reconstruct growing season temperatures (Guiot et al., 2010). They are below the average RE (0.73) obtained by Neukom et al. (2010a) in a principal component reconstruction of summer DJF temperatures in SSA, but they are comparable to those found in the precipitation reconstruction (average RE around 0.27 in both seasons, Neukom et al., 2010b). Our RE values in the Mediterranean region of the Andes (around 0.35) are also comparable to the RE calculated by Christie et al. (2009) in a tree-ring reconstruction of PDSI in this area (around 0.45). However, Christie et al. (2009)'s reconstruction focused on high frequencies, while our multi-proxy reconstruction aims at reconstructing the full spectra of PDSI variations. Therefore, the comparison between RE statistics should be interpreted with care, as errors propagate differently when HF and LF bands are reconstructed simultaneously.
The capacity of the SAM to reproduce modern spatial patterns of variations of the PDSI in SSA is another indication of the performance and usefulness of our modeling approach. Figure 6 presents the maps of the average observed and reconstructed PDSI values for two consecutive 30-yr windows (1931-1961 and 1961-1993). The first period  is marked by wet conditions in PG and the ANDES and a drier climate in PA and ST. This pattern is also very clearly  seen in our reconstruction. The second period  has the opposite pattern: a wet climate in PM and ST and drier conditions in PG and the ANDES. This spatial inversion of PDSI trends was well captured in our reconstruction.
As an independent validation for our work, we checked the correspondence between our reconstruction and the 1767-2002 reconstruction of Laguna Mar Chiquita water-levels fluctuation (Piovano et al., 2002). We used the average PDSI within the lake's watershed for a comparison (60-65 • W; 20-25 • S). The agreement between the two series is fairly good (Fig. 7), at least in the mid to low frequency domains, considering that both reconstructions were performed from totally independent proxies. The abrupt rise of Mar Chiquita's lake   . 9. Comparison between the reconstructions made from the full proxy dataset (colored lines, same colors as in Fig. 3) and the reconstructions performed from only the proxies that go back to AD 1000 (black lines). The top and middle rows present the LF and HF components, respectively, while the lower row depicts the comparison over the full spectra. Correlation coefficients (R, N = 994) are shown in the title, for each frequency band and for each region. Not shown here: ANDES (R = 0.28, for the full spectra) and ST (R = 0.57, for the full spectra).
levels during the 1970s also corresponds to a rise in summer PDSI, although the latter has a smaller amplitude. It is worth mentioning that a recent modeling of Mar Chiquita's water levels has shown that the recent rise is attributable to an increase in runoff in the northern sub-basin (Troin et al., 2010), evoking an increasingly dominant tropical influence in the area. Thus, since the PDSI is a value that is standardized to reflect departures from local mean conditions, the average PDSI within the mar Chiquita watershed may not adequately reflect this possible flow redistribution. Before exploring millennial summer PDSI fluctuations, it is important to establish that our reconstruction is reliable over the full period. To do so, we provide a series of nested reconstructions and their corresponding verification (RE) statistics. Nested reconstructions were computed exclusively from subsets of proxies that are older than AD: 1000AD: , 1025AD: , 1200AD: , 1222AD: , 1388AD: , 1505AD: , 1649 respectively corresponding to the calendar dates to which a new LF proxy is added to the reconstruction. We present the evolution of RE statistics in each sub-region on Fig. 8. Even with a limited number of proxies at the beginning of the last millennium, the RE statistics remains positive in all regions, meaning that the spectral analogues have good predictive skills over the last thousand years. As an example, we graphically compare the full reconstruction and the AD 1000 nested reconstruction ( Fig. 9) for two regions: PG and PM. Our analysis shows that LFs are comparable between reconstructions over the last thousand years. However, HFs are less similar and this probably relates to the fact that higher frequencies contain a lot of the local climatic signal (i.e. noise) that cannot be adequately reconstructed from a very limited number of HF proxies. In conclusion, the correlations for the full spectra remain acceptable, suggesting that the long-term trends can be interpreted since AD 1000.

PDSI fluctuations over the last millennium in SSA
Summer PDSI reconstructions (1000-1993) and regime shift detection analysis (extended until AD 2005 using Dai et al. (2004)'s PDSI data) are presented in Figs. 10 and 11. The period between 1000 and ∼1250 clearly appears as a distinct regime in all areas of SSA. Breakpoints corresponding to the end of the latter period were identified in 1242 (PG) 1241 (PM), 1234 (ST), and 1255 (ANDES) (Fig. 11). In most regions except PG, the first part of the millennium was probably slightly wetter that today. This period corresponds to the end of the MWP that extended until 1200-1350 in the Northern Hemisphere (Jansen et al., 2007;Mann et al., 2009;Guiot et al., 2005Guiot et al., , 2010DahlJensen et al., 1998), and until about 1200-1250 in the Southern Hemisphere . In SSA, the MWP seemed to have persisted until about ∼1350 (Neukom et al., 2010a) 100-yr lag between our PDSI reconstruction and Neukom et al. (2010a) summer temperature reconstruction needs further investigation but could possibly be explained by changes in precipitation patterns that are not present in temperature reconstructions. It is interesting to note, however, that from the perspective of PDSI, the first part of the millennium in SSA was probably a period of important geographical contrasts rather than a widespread drought.
The ∼1250-∼1380 period is also clearly distinct in all regions of SSA and is characterized by drier than normal conditions (Figs. 10 and 11), except in PG where wetter conditions prevailed. Again, regime-ending dates cluster in time: 1376 (PG), 1375 (PM), 1386 (ST), 1384 (ANDES), and 1369 (SSA). Interestingly, this well-defined regime corresponds to the Wolf Minimum (Eddy, 1976), the first period of the LIA with an extremely low solar activity. The physical link between solar irradiation and the PDSI in SSA, however, needs to be substantiated by modeling studies in order to clarify the processes and forcings involved. Additionally, the possible climatic impact of volcanic eruptions should be investigated, Regime shift detection is based on a sequential a t-test algorithm.To be detected, a regime needs to be at least 50 yr long. If a regime has less than 50 yr, the method can still allow for its detection; however, the significance required to determine that a date corresponds to a changepoint is adjusted to be inversely related to the length of the regime. The bottom panel corresponds to the 20 yr running proportion of the number of series, recording a changepoint at year t. It gives a general idea of where these changepoints cluster in time. Color code from Fig. 3 is used here.
as a series of major eruptions occurred during that period (Wanner et al., 2008). The physical processes driving these changes are currently being examined through a model-data comparison study.  The 1780-1820 is also an important regime change period (Figs. 10 and 11) and corresponds to the Dalton minimum (Eddy, 1976), though its spatial extent is less marked that the two previously described periods. Wetter than normal conditions were recorded in PG, while drier conditions were reconstructed in PM and ST.
Although important regime changes occurred over the last centuries, most fell within the bounds of the last millennium variability (Fig. 11). The best example is the drying episode that occurred in PG since ∼1880. Similar regime conditions could be found between 1400 and 1750 in the area, and even before ∼1250. We must point out that, while individual years within the regimes are reconstructed by analogy using the SAM, regimes drier or wetter than today can actually be reconstructed if, for example, extreme years (either wet or dry) cluster in a given past period. This was rarely the case in our analysis. Nevertheless, several exceptions exist. The recent wet regime in PM seems rare (if not exceptional) at the scale of the last millennium; so is the dry regime observed in the ANDES since the 1930s that has very few (if no) equivalent in the past.
A striking feature of our reconstruction is an antiphase between PG and the rest of the subcontinent. A comparison of 50-yr-filtered reconstructions (Fig. 12) reveals that SSA's summer PDSI evolution, rather then being uniform over all the studied area, was characterized by an important geographical contrast, especially between the northeast and the southernmost part of the continent. That contrast can be most easily perceived in the low frequency domain (f < 0.03), but it characterizes, typically, all frequency bands below 0.2 (T = 5yr). In order to better interpret that antiphase, we compared our reconstructions of contrasted PG and PM regions to those of Neukom et al. (2010a and2010b) (Fig. 7). In PG, summer PDSI seemed to be more responsive to temperatures than to precipitation. In this area, warm summers and winters were associated with a wet climate (positive relationship). By contrast, in PM, PDSI seemed closely linked to DJF precipitation. Moreover, DJF and JJA temperatures related differently to summer PDSI. Warm summers were associated with a dry climate while warm winters were generally coupled to wet summer conditions. These contrasting dynamics underline the fact that the PDSI might not have responded similarly to precipitation or temperature in every region, but instead that variations in the drought index were probably driven by different parameters whose importance are likely to vary through space.
A frequency analysis was performed to determine whether or not extremely wet (or dry) spells have been more common (or rarer) in the past. For each region, a reference PDSI value was defined. That value corresponded to the 50-yr wet spell and drought identified for the 1930-1993 period. Then, past return periods ( = 1/probability to exceed the reference PDSI value) were retrieved after fitting a log-normal distribution to the 100 yr preceding each year t. The results are presented on Fig. 13  . Past return periods (Ω) (and their 95% bootstrap confidence interval) of events equivalent in magnitude to the reference 50-yr wet spell (left panel) or drought (right panel). The reference 50-yr event (value indicated in parenthesis for each region) was calculated for the 1930-1993 period. All return periods were obtained after fitting a log-normal distribution to the 100 years preceding year . Color codes for each region are the same as in Figure 3. Return period values on the y axis are exponential values ie. 5e years =5*10 years. Thus values of 1,2 and 3 respectively yield 50-yr, 500-yr and 5000-yr return periods. Confidence intervals are incremented at each 5%, with reversed color directions for wet spells and droughts. Color codes for each region are the same as in Fig. 3. Return period values on the y axis are exponential values i.e. 5ek years = 5 × 10 k years. Thus k values of 1, 2 and 3 respectively yield 50-yr, 500-yr and 5000-yr return periods. Confidence intervals are incremented at each 5reversed color directions for wet spells and droughts. 50-yr wet spell generally had a larger between 1250 and 1400 ( > 500 yr), except in PG where their was closer to the modern value ( ∼100 yr). Between 1400 and 1600, a tendency towards shorter (<100 yr) can be observed everywhere except in PG where very rare occurrences are noted ( > 500 yr). From 1600 to the 20th century, wet events progressively became more common, except in PM where they were rarer. The latter tendency reversed after 1900 in some regions: the ANDES (wet spells were rarer after 1970), PM (wet spells were more common after 1970), and PG (wet spells were rarer after 1900).
The frequency analysis of droughts revealed the antiphase as well (Fig. 13). In general, extreme droughts were common before 1400 in most regions ( < 100 yr), except in PG where they were clearly rarer between 1250 and 1400 ( > 500 yr). Between 1450-1700, extremely dry events were rarer than today ( > 500 yr) almost everywhere in SSA, except in PG where an opposite trend can be observed. An interesting point is that although the antiphase between PG and SSA seems symmetrical for both droughts and wet spells, the behaviour of droughts with respect to their wet spell counterparts is not necessarily symmetrical in each region, possibly indicating a relative independency between the two phenomena. A good example of that independency is found in PM after 1400. While wet spells became less frequent in PM during that period (possibly indicating drier conditions in the area), extreme droughts also became less frequent between 1400 and 1600. Afterwards, severe droughts in PM became more common (20 < < 100 yr) until about 1900, while wet spells were still quite rare (200 < < 1000). This asymmetry highlights the fact that, regardless of the long-term trends in the mean of summer PDSI series, extremely dry and/or wet spells can still occur isolatedly in small clusters and independently from the regime conditions. In other words, the driest (wettest) years do not systematically occur during the driest (wettest) regimes.

Links with ocean-atmosphere forcings (AAO, ENSO, PDO)
The climate of SSA is under the influence of high latitude and tropical ocean-atmospheric forcings, but the magnitude and the direction of this influence varies considerably between regions. To explore the teleconnections between summer PDSI and the dominant indices in the area (summer means for AAO, ENSO, and PDO, all retrieved from the Climate Diagnostics Center (NOAA) at http://www.esrl.noaa.gov/psd/ data/), we first used standard Pearson correlation analysis. The AAO index is estimated by the first PC of the 850 hPa geopotential height anomalies south of 20 • S ( Thompson and Wallace, 2000). The ENSO index corresponds to the mean SST anomalies from the N3.4 Pacific region and is a common proxy for the ENSO phenomenon (Trenberth, 1997). The PDO (Zhang et al., 1997) is an ENSO-like phenomenon that exhibits decadal to interdecadal variability (at least 20 to 30 yr). The estimation of the PDO interdecadal variability remains poorly known because very few stations in South America are century-long. In order to evaluate our reconstruction's ability to reproduce climatic patterns driven by the major indices in SSA, reconstructed PDSI values were used as predictands. However, we also tested the relationships with instrumental PDSI values and the results are comparable (Fig. 14). Over the common period , AAO is well correlated to summer PDSI in SSA. In PG and the ANDES, the correlation is negative and reaches R = −0.53, while in PM and ST, the correlation is positive (Fig. 14)   variations (Fig. 12) and trends in the extremes (Fig.13) of summer PDSI in PG have been in antiphase with those on the rest of the sub-continent. It is now possible to argue that these antiphases could have been driven by past lowfrequency variations in the AAO, since this index is the only one that can be associated to such a contrasted PDSI response. As underlined by Garreaud et al. (2009), the AAO is the leading pattern of tropospheric circulation variability south of 20 • S. Its influence on temperature over the recent period is unequivocal over the recent period (Gillett et al., 2006;Garreaud et al., 2009). The most striking aspect is a contrasted response of annual surface temperatures with a clear summer warming south of 40 • S and a cooling elsewhere during positive phase of the AAO (Garreaud et al., 2009). These results highlight the fact that summer PDSI variations in SSA are responsive to AAO-induced temperature fluctuations at the continental scale, a result that is also supported by the comparison with Neukom et al. (2010a)'s temperature reconstruction. Summer PDSI is also related to the ENSO phenomena in SSA . However, the response is far less con-trasted than with AAO and mostly positive throughout the studied area (Fig. 14). The strongest correlations are found in the Mediterranean area of the Andes and along the eastern part of Argentina. Positive (negative) ENSO years typically generate conditions that are wetter (drier) that the normal in SSA. PDO is also related to the fluctuations of summer PDSI . This index has a general positive influence, but the highest correlations (R = 0.5) are found in the Mediterranean Andes region. Much like with ENSO, positive (negative) PDO years generate conditions that are wetter (drier) than the normal almost everywhere in SSA.
Finally, it is clear that all these ocean-atmosphere forcings act synergically to influence summer PDSI in SSA. In order to better understand the interactions between AAO and ENSO, model-data comparison studies will be necessary. Here, we simply underline the fact that SSA's summer PDSI can be compartmentalized according to the dominant ocean-atmosphere forcing. To determine which index has a predominant influence over the 1950-2005 period, we performed a regression analysis using each index as a predictor. For each 2.5 • × 2.5 • pixel in SSA, the equation had the . Red symbols represent a dominant AAO influence while blue and green symbols indicate a dominant ENSO or PDO influence, respectively. Filled (empty) squares denote a positive (negative) influence on summer PDSI. The regression was performed for the longest common period . Spots with no symbols had a negative RE value.
form PDSI = a(AAO) + b(ENSO) + c(PDO) + noise, where a, b, and c are coefficients. For each regression, we retained the largest coefficient (absolute value) and noted its sign (Fig. 15). Interestingly, this figure clearly shows that the various regions of SSA are under the influence of different dominant indexes. In PG and northeastern ST, the dominant process is clearly the negative and spatially coherent influence of AAO. A positive influence of AAO is also dominant west of 70 • W in PM and in lower ST. Both PDO and ENSO are both dominant in the Mediterranean ANDES, in western PM, and in ST. The dominant influences appear to be spatially coherent (points with a common influence tend to cluster in space), suggesting a true physical link between PDSI variations and these atmospheric circulation modes. The question remains whether or not these influences have been stationary through time and how did they interacted through time. Data-model comparison studies over the whole millennium are necessary to clarify these issues.

Conclusions
This study presents the first spatially explicit summer PDSI reconstruction in the Southern Hemisphere. We provide a 2.5 • × 2.5 • gridded reconstruction that extends the Dai et al. (2004) PDSI dataset over the last millennium at the subcontinental scale, south of 20 • S in SSA. Calibration and ver-ification statistics, along with the spatial analysis of modern patterns of changes, show that in most areas of SSA (PG, PM, ST, ANDES), the summer PDSI is well reconstructed using the SAM. Over the last millennium, SSA has experienced considerable PDSI variations, some of which might have been at least as important as those recently observed in the area. The temporally well-defined 1000-1250 period possibly associated with the MWP was characterized by wetter than normal conditions in most studied regions, except in PG where the climate was clearly drier that the normal. This regime terminated very abruptly around ∼1250, as indicated by the regime shift detection and the analysis of extremes. In most regions, the climatic spatial pattern reversed between 1250 and ∼1400. SSA's climate became much drier, except in PG where the climate humidified during the austral summer. Afterwards, SSA slowly humidified to approach conditions somewhat closer to the normal (although with some important discrepancies between regions, e.g. in PM the climate dried up).
The analysis of recent PDSI fluctuations in the context of the last millennium reveals analogous patterns for the modern period over the last thousand years, especially between 1000 and 1400. Our reconstruction also shows that PG's fluctuations were generally in antiphase with the rest of the continent. Such an antiphase could be driven by the Antarctic Oscillation (AAO), providing that modern patterns of response are transferable to the past. We reveal evidences that the AAO has a contrasted effect on the PDSI in SSA over the calibration period. During its positive phase, the climate tends to be humid in the northeastern part of SSA and much drier in the south. During its negative phase, the latter situation reverses. However, AAO has never acted alone to modulate fluctuations in the mean and extremes of PDSI. Nevertheless, in some regions like PG, AAO clearly is the dominant driver (at least during the calibration period), while in others it plays a smaller role.
Our analysis finally shows that past low-frequency PDSI variations can be reconstructed quite successfully in most regions of SSA. We have shown that even with a reduced number of proxies, low frequencies reconstructed using the SAM remain quite similar to the low frequencies reconstructed from the full dataset. However, we have also shown that high frequency PDSI variations are less well reconstructed from a reduced set of proxies. In order to better reconstruct high frequency variations of PDSI in SSA, more highly resolved proxies are needed, as stated by Villalba et al. (2009), especially in proxy-lacking areas such as ST area and eastern PM.
Finally, future analysis should focus on studying the interactions between ocean-atmosphere forcings (e.g. interactions between AAO an ENSO) in order to better identify how they can modulate SSA's climate in both time and space. Such work should preferentially be achieved through climate models.  Fig. A1. Goethite (derivative 445 nm) chronology of Lago Fras varved lake sediments. Stars represent some well-known earthquakes in the area (AD 1960;AD 1751, andAD 1737). Figure A1 represents the Goethite (derivative 445 nm) varved chronology of Lago Frías sediments . This sedimentary proxy reflects the amount of precipitation-induced erosion in the watershed. The dating of these sediments was confirmed by three well-known earthquakes (1737, 1751, and 1960) also shown on the figure.