Multi-century lake area changes in the Southern Altiplano : a tree-ring-based reconstruction

Size fluctuations in endorheic lakes in northwestern Argentina (NWA) and southwestern Bolivia (SWB) are very sensitive to basin hydrological balances, and consequently, very vulnerable to deleterious effects from climatic changes. The management of these water resources and their biodiversity requires a comprehensive knowledge of their natural variability over multiple timescales. In this study, we present a multi-century reconstruction of past lake-area fluctuations in NWA and SWB. The evidence used to develop and validate this reconstruction includes satellite images and a century-long tree-ring record from P. tarapacana. Interannual fluctuations in lake area of nine lakes were quantified based on Landsat satellite images over the period 1975 to 2009. A regional P. tarapacana tree-ring chronology, composite from two sampling sites, was used as predictors in a regression model to reconstruct the mean annual (January– December) lake area from the nine lakes. The reconstruction model captures 62 % of the total variance in lake-area fluctuations and shows adequate levels of cross-validation. This high-resolution reconstruction covers the past 601 years and characterizes the occurrence of annual to multi-decadal lake area fluctuations and its main oscillation modes of variability. Our reconstruction points out that the late 20th century decrease in lake area was exceptional over the period 1407–2007; a persistent negative trend in lake area is clear in the reconstruction and consistent with glacier retreat and other climate proxies from the Altiplano and the tropical Andes. Since the mid 1970s, the Vilama-Coruto lake system recorded an accelerated decrease in area consistent with an increasing recurrence of extremely small lake-area events. Throughout the 601 years, the reconstruction provides valuable information about spatial and temporal stabilities of the relationships between changes in lake area, ENSO, and PDO, highlighting the Pacific influence over most modes of lake area variability. Global and regional climate models for the Altiplano project a marked reduction in precipitation to the end of the 21st century, exacerbating presently dry conditions. These results provide a baseline for the historical range of variability in lake fluctuations and thus should be considered for the management of biodiversity and water resources in the Central Andes during the next decades.


Introduction
The southern Altiplano or dry Puna, is located in the southern subtropical Andes.A prominent feature of the landscape between 21.5 and 24 • S is the shallow salty lakes in the lowest point of endorheic watersheds at very high elevation (above 4400 m a.s.l.).
Productivity in the dry Puna is low and mainly concentrated in the wetlands (lakes and peat bogs), which consequently play a critical role in sustaining a unique diversity of rare and endemic biota (Squeo et al., 2006).Bird and large vertebrate species depend upon the wetlands for grazing, nesting and water.This is an extreme environment with a prolonged dry season, strong winds, high daily temperature amplitude with frequent frosts, intense solar radiation, and hypoxia.These conditions prevent agriculture development and permanent grazing, resulting in the absence of human settlement above 4000 m.During summer, however, the zone is used as grazing area for livestock, which is an important component of the local indigenous economy (Nielsen, 2003).The ecological dynamics of these wetlands are sensi-tive to the watershed's hydrological balance, related mainly to precipitation and temperature (Carilla et al., 2013); and lake ecosystems act as direct and indirect indicators of climate variability.
During the 20th century, the tropical and subtropical Andes have experienced a persistent warming trend (Vuille and Bradley, 2000).Climate scenarios for the 21st century predict a decrease in precipitation and an increase in temperature (Urrutia and Vuille, 2009), posing a major threat for water resources in this arid region.Changes in lake area, reflecting hydrological dynamics and vegetation productivity (Carilla et al., 2013), allow the assessment of climate change effects in areas and periods with poor or non-existing instrumental records.Lake hydrological dynamics are usually measured as fluctuations in lake levels recorded at gauging stations but this region lacks such instrumental data.However, interannual changes in lake area between 1985 and 2009 were quantified based on Landsat satellite images and their association with tree-ring and instrumental climatic records (adjacent to our study area) was analyzed (Carilla et al., 2013).
P. tarapacana are small trees that grow between 4200 and 5200 m a.s.l. on the slopes of the volcanoes near the lacustrine areas studied.This extremely moisture-sensitive woody species has allowed the development of a network of multicentury tree-ring width chronologies along the South American Altiplano from 16 to 23 • S (Morales et al., 2012).As the yearly growth of Polylepis trees varies depending on hydrological balance (Solíz et al., 2009;Christie et al., 2009), annual rings of this species may therefore be used to reconstruct past lake area fluctuations.
The main goal of our study was to use P. tarapacana ring-width records and Landsat images to develop a highresolution multi-century reconstruction of past fluctuations in lake area from the Vilama-Coruto region in northwestern Argentina (NWA) and southwestern Bolivia (SWB).We describe temporal fluctuations of lake area and put recent changes and trends in the context of the past 6 centuries.We compare this reconstruction with palaeo-hydroclimatic records from the subtropical and tropical Andes.In addition, we analyzed the recurrence of extreme drought events, the presence of persistent periodicities and dominant oscillation modes of variation at different timescales.Finally, in order to identify the major climatic factors influencing variations in lake area, we compared the Vilama-Coruto lake reconstruction with indices of El Niño-Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO).

Study system
The study area is located on the border between Argentina (Rinconada, Jujuy) and Bolivia (Sud Lipez, Potosí) from 22 • 10 to 22 • 45 S and from 66 • 25 to 67 • 13 W (Fig. 1).This area shows the influence of prolonged volcanic activ- ity during the Upper Cenozoic (Coira et al., 1996) with the presence of a large concentration of shallow lakes that form an interrupted chain of endorheic watersheds above 4400 m.The study area belongs to the broader ecological region called dry Puna, located in the high-elevation plateaus of the southern subtropical Andes (Cabrera et al., 1976).Vegetation cover does not exceed 20 % and is dominated by Poaceae (Festuca, Deyeuxia), several small shrubby species (Parastrephia, Adesmia, Acantholipia), and cushion plants (Azorella, Pycnophyllum) (Cabrera et al., 1976).P. tarapacana is a small tree, ca 2-3 m tall, and is the largest woody species that grows on the slope of the high volcanoes between 4200-5200 m.Extreme aridity (< 150 mm annual rainfall) and low temperatures (mean annual < 5 • C) characterize this region (http://www.worldclim.org).Precipitation is of convective type and mostly concentrated during summer.More than 80 % of the total annual rainfall falls between November and March (Vuille et al., 2003;Vuille and Keimig, 2004).Summer precipitation variability is related to local orography and with changes of the upper troposphere zonal wind, which in turn favor or block the ingression of wet easterly flow transporting humid air masses from the Amazon basin (Vuille et al., 2000;Garreaud and Aceituno, 2001;Bradley et al., 2003;Garreaud et al., 2003).

Lake area records
Between 1975 and 2009, inter-annual fluctuations in lake area from nine lakes were quantified using Landsat images.The nine lakes are located between 4400 and 4600 m in the Vilama-Coruto region from NWA-SWB (Fig. 1).Lake area variations were derived from Landsat MSS (80 m pixel resolution, path 249 row 076) for the period 1975-1982, and Landsat TM images (30 m pixel resolution path 232 and row 076) for the period 1984-2009.We calculated lake area based on the number of pixels and their sizes.Pixel size for Landsat MSS is 80 × 80 m (0.0064 km 2 ) and for Landsat TM is 30 × 30 m (0.0009 km 2 ).For example in November 1978 Cerro Negro lake encompassed 316 pixels of Landsat MSS corresponding to 2.02 km 2 (316 × 0.0064 km 2 ).In October 1988 the same lake, encompassed 2551 pixels of Landsat TM image, which is the same as saying 2.3 km 2 (2551 × 0.0009 km 2 ).The area of the lakes in each particular image was determined using the non-parametric method support vector machine (SVM; Hsu et al., 2007).Good quality available images for the different years of the period 1975-2009, ranged from one to ten (dates of images are in the Supplement S1).Despite the variable acquisition dates and number of the images, relative low intra-annual area variability was recorded for the nine lakes.A comparison between the coefficients of variation for each lake showed that interannual variability represented more than double the intraannual variability (see the Supplement S2.1-9).Therefore, we assume that the monthly or seasonal variation was negligible in relation to the higher inter-annual lake area variation, and an area value from a particular image could well represent the mean lake area condition of the year.To further validate this assumption, autocorrelation analyses were performed to determine persistence in the monthly and annual time series averaged for the nine lakes from 1975-2009 period (see the Supplement S3).This analysis showed strong 19-month persistence in the monthly series, while 2-year persistence was found in the annual series.AnClim program (Štěpánek, 2008) was used to perform this analysis.
The annual (January-December) lake area record for each of the nine lakes was estimated by averaging every available image (from one to ten) for a particular year.
Based on the nine individual records, we developed a regional record of annual lake area fluctuations.To minimize the influences of lakes with different sizes on the regional mean, annual lake area records were standardized as Z score.Finally, the standardized lake area records from the nine lakes were averaged to develop a regional mean covering the period 1975-2009.There was no image available in the entire region for the year 1983.The model used for the calibration and verification of the past lake area reconstruction requires continuous time series without missing years.Therefore, a linear regression model between the regional lake area and a regional precipitation index from the Altiplano was applied to estimate the 1983 lake area value.A list of the meteorological stations used for this analysis is reported in the Supplement (S5).

Tree-ring data collection and standardization
Our study area encompasses the southern latitudinal limit of P. tarapacana distribution and tree stands are mostly scattered on wetter north and east slopes of high volcanoes but not on the south and west slopes.Sample sites include the P. tarapacana populations located in relatively close proximity to the study lakes.Cross-sections and wedges were obtained from dead and living trees.From a total of seven sampled sites, we selected those that were highly replicated (> 60 tree-ring series), contained at least 300 years of records, and were highly correlated with lake area fluctuations (r > 0.5; Table 1).Two sites met these criteria: the tree-ring records from the Granada and Uturunco volcanoes (Table 1).The Schulman's convention (1956) for the Southern Hemisphere was used to assign to each ring the date of the calendar year in which radial growth starts.We measured total ring width to ±0.001 mm and cross-dated samples using standard procedures (Stokes and Smiley, 1968;Fritts, 1976).The program COFECHA (Holmes et al., 1983) was used to detect measurement and cross-dating errors.The ring-width chronologies from Granada and Uturunco shared large percentage of common signal in tree-ring variations (r > 0.62; 1620-2006 common period), therefore, 155 tree-ring series were merged into a regional chronology.Ring-width measurements were standardized to remove the age-related growth trends and minimize the growth variations not related to climate (Fritts, 1976).A negative exponential curve or straight line was used to standardized ring-width series.When this standardized growth trend is removed, some of the variance related to climatic forcing signal can also be removed, leading to a trend distortion in resulting index series (Melvin, 2004).To prevent trend distortion the signal free detrending procedure (Melvin, 2004;Melvin and Briffa, 2008) was applied to the tree-ring measurements.The regional tree-ring chronology was produced with the RCSigFree program (Tree Ring Lab-LDEO, Columbia University http://www.ldeo.columbia.edu/tree-ring-laboratory/resources/software).
The quality of the tree-ring chronology was tested based on standard chronology diagnostics such as the RBar (Briffa, 1995) and the Expressed Population Signal (EPS; Wigley et al., 1984;Table 1;S7).The RBar is the mean correlation coefficient from all possible pairings among individual cores in the chronology (Briffa, 1995).EPS is a measurement of the strength of the common signal in the chronology used to identify the most reliable interval in the chronology valid for climate reconstructions.To calculate the RBar and EPS, we used a 50-year window with an overlap of 25 years between adjacent windows.

Lake area reconstruction
First, we identified the best relationship between variations of regional lake area and inter-annual tree growth.For that purpose, we computed correlation coefficients between the regional standard P. tarapacana chronology and annual mean lake area records.The highest correlation coefficient with tree-ring variations was found using for comparison the regional mean annual (January-December) lake area (r = 0.70; p < 0.0001).We reconstructed, using a linear regression, the mean January-December lake area from the Vilama- Coruto lacustrine region over the period 1407-2007.To capture possible climate-related persistence in the tree-ring series and enhance the common climate signal, we included as predictors in the regression the regional chronology in the temporal lags significantly correlated (p < 0.05) to annual lake area during the calibration period 1975-2007.The chronology at year t and a lag (t + 1) were considered as candidate predictors of annual lake area in a stepwise multiple regression (Weisberg, 1985).The "leave-one-out" crossvalidation procedure (Michaelsen, 1987;Meko and Baisan, 2001) was used to validate the model (Fig. 2).This method allowed the use of all available data to calibrate and validate the model, and it is very useful when observed records are short temporal series, such as our lake area records.In this approach, the reconstruction model is calibrated using the full available overlap period.Then, the method sequentially withheld one observation and calibrates additional models with the remaining observations.Each model is used to estimate the withheld observation.Finally, all these estimated values are assembled in a time series used as the validation data set.The quality of the model during the calibration period was measured by the proportion of variance explained by the regression (R 2 adj = 0.62).The model accuracy was evaluated using the reduction of error (RE) statistic, which measures the model skill in estimating values of the predictand in the cross-validation process (Wilkes, 1995;Fritts, 1976).We also computed the root mean square error (RMSE) statistic as a measure of inherent uncertainties in the reconstruction (Weisberg, 1985).

Identification of extreme high/low events in the lake area reconstruction
Significant shifts in the mean state conditions of the lake area reconstruction were identified using a simple regime shift detection technique (Rodionov, 2004) with a cut-off segment length l = 25 years, a target probability level p = 0.1 and a Hubert's weight parameter h = 1 (Fig. 3).In order to determine the intensity and duration of drought (pluvial) events in the lake area from the Vilama-Coruto region, we computed the lowest (highest) reconstructed n-year running means for n = 1, 5, 25 and 50 years (Table 2).In addition, we analyzed the occurrence rate of extreme reduction in lake area within the reconstruction (Fig. 4).For that purpose, we employed a kernel estimation technique (Mudelsee et al., 2004;Girardin et al., 2006).This method allows the detection of nonmonotonic trends and imposes no parametric restrictions.A Gaussian kernel function was applied to weigh reconstructed years of extreme decrease in lake area and calculate the occurrence rate (probability per time unit).A Kernel bandwidth of 50 years was used and pseudo-data generation rule "reflection" to reduce boundary bias (Mudelsee et al., 2004).Confidence bands at the 90 % level were obtained by using 2000 bootstrap simulations (Cowling et al., 1996;Mudelsee et al., 2004).The dates of extreme events were defined by a percentile threshold of 20 %.    281 (1908-1933) 1.164 (1908-1957) the reconstructed annual (January-December) lake area and the annual averaged (January-December; 2.5 × 2.5 gridded cell) sea surface temperatures from the NCEP reanalysis global data set (Kistler et al., 2001: http://www.cpc.ncep.noaa.gov/data/indices/;Fig. 5).Additionally, power-spectral and cross-spectral analyses (Jenkins and Watts, 1968) were used to determine the frequency domains of the reconstructed lake area variations in the Vilama-Coruto region and the annual averaged Pacific sea surface temperatures from the N3.4 Pacific sector (N3.4-SSTs;Fig. 6).Blackman-Tukey (BT) spectral analysis (Jenkins and Watts, 1968) was computed for the periods 1407-2007 and 1872-2007 for the reconstruction and the Pacific N3.4-SSTs, respectively.The number of lags used for the auto-and crosscovariance functions was equal to 30 % of the series length.Spectral functions were smoothed with the Hamming filter.The 95 % confidence level of the spectrum was estimated from a "red noise" first-order Markov null continuum based on the lag-1 autocorrelation of the time series (Mitchell et al., 1966).A coherency spectral analysis was computed over the common period 1872-2007 between the reconstructed lake area and N3.4-SSTs records.The coherency spectrum measures the variance in common between reconstructed lake area and N3.4-SSTs as a function of frequency.Both, BT and coherency analyses were performed using the AnClim program (Štěpánek, 2008).To detect the dominant, non-stationary oscillatory modes of variability in the 601 years series of reconstructed lake area we used singular spectral analysis (SSA; Vautard and Ghil, 1989; Fig. 7).SSA is a non-parametric method related to empirical orthogonal function analysis that samples lagged copies of time series at equal time intervals and calculates eigenvalues and eigenvectors of the auto-covariance matrix (Vautard and Ghil, 1989;Vautard, 1995).SSA enables the user to evaluate the changing periodic behavior in a single time series by extracting "reconstructed components" or "waveforms", which represent the dominant periodic modes of the time series.This method detects changes in amplitude and phasing of the dominant cycles over time because waveforms retain their temporal periodic behaviour (Vautard and Ghil, 1989;Vautard, 1995;Speer et al., 2001).In addition, SSA was computed for the N3.4-SSTs and PDO indices to assess the relationship in the time frequency domain between the reconstructed lake area and Pacific climate forcing.The SSA analysis was performed by using the SSA program (Boninsegna and Holmes, unpublished operating manual [on file at Laboratory of Tree-Ring Research, University of Arizona]).

Tree-ring chronology and reconstruction model
In this study, 155 tree-ring samples from dead and living trees of P. tarapacana were combined to develop a robust composite tree-ring chronology (S6) for the Vilama-Coruto region in NWA and SWB.The regional chronology covers the period AD 1242-2008, but is best replicated since AD 1407 (N > 10 ring-width series).As previously demon- strated (Morales et al., 2004(Morales et al., , 2012;;Carilla et al., 2013), P. tarapacana trees are drought sensitive and show relativehigh variability in tree-ring widths (mean sensitivity = 0.30).The RBar and EPS statistics indicate that the chronology is of good quality with a strong common signal in ring width between individuals (Table 1).The chronology statistics are similar to those previously documented for P. tarapacana tree-ring records from the Bolivian Altiplano and the Argentine Puna (Morales et al., 2004;Solíz et al., 2009;Christie et al., 2009;Morales et al., 2012;Carilla et al., 2013).The EPS remains above the threshold of 0.85 over the period 1399-2008.The mean EPS for the period 1407-2008 is 0.96 (Table 1).The "signal free" standardization method, using negative exponential curves, allows to capture appreciable amounts of both low-(multi-decadal) and high-frequency (inter-annual) variability in lake area fluctuations.
We present here a 601 year high-resolution tree-ring reconstruction of annual area fluctuations from nine lakes in the Vilama-Coruto region.Our reconstruction, calibrated and validated on satellite image January-December lake area ; Fig. 2), accounts for 62 % of the total variability in lake area changes (R 2 adj adjusted for loss of degreesof-freedom; Draper and Smith, 1981).The consistence between the observed and estimated lake area indicates good quality in the calibration model for reproducing past lake area fluctuations.The high predictive power of the calibration model is consistent with a highly significant F value (F = 25.4;P < 0.001) and a largely positive Reduction of Error (RE = 0.54). Aalyses of regression residuals indicated no violations of regression assumptions.The residuals of the regression models are normally distributed, not significantly autocorrelated and without a significant trend according to Durbin-Watson tests (Fig. 2b).

Temporal evolution of the lake area reconstruction
The lake area reconstruction for the Vilama-Coruto region is characterized by inter-annual variations embedded within decadal to multi-decadal fluctuations (Fig. 3).Significant persistent changes in the mean state of the reconstruction were detected with the Rodionov analysis (Rodionov, 2004;Fig. 3).This analysis shows that lake-area conditions around the long-term mean for the period 1407-1440 were followed by a lake area increase lasting to the beginning of the 16th century (1503).The 15th century high lake level stand was characterized by high inter-annual variance, with extreme high-and low-lake area years.The longest period of comparatively low-lake area spans most of the 16th century (1504-1595).This long-term dry period was interrupted by 2 decade-scale periods showing contrasting high and low lake areas centered on the years 1600 and 1620, respectively.After that, significant high lake area conditions prevailed for more than 160 years (1627-1790).Two of the five years with the largest lake areas were recorded during this period (1734 and 1744; Table 2).The 35-year smoothing spline shows a period of gradual reduction in lake area between 1750 and 1790).However, no significant shifts in mean lake-area extent were detected over this interval.For the period 1797-1818 low lake area conditions were identified.A marked shift from low-to high-lake areas was identified around 1820.The period 1820-1859 encompasses the highest lake areas reconstructed during the past 601 years.After this long-term increase in lake area, a gradual decrease in lake areas was reconstructed to the first decades of the 20th century.After that, a steady negative trend in lake area begins in the 1930s and persists to the present.The most pronounced negative decline for the entire reconstruction is recorded for the period 1976-2007.During this 30 year period 1987 was the only year to exceed the long-term historical mean area.Three of the five years with the lowest lake areas occurred during the last 3 decades (Table 2).
A different picture of the long-term history of changes in lake area periods is given by averaging non-overlapping intervals of 5-, 25-, and 50-years (Table 2).Based on a 5year average, three of the five periods of lowest lake area occurred since 1930 (i.e. 1938-1942, 1998-2002, 2003-2007), whereas the highest lake-area periods in the reconstruction were centered around the 1740s and 1840s (Table 2).In terms of 25-year averages, the period 1983-2007 has the smallest 25 year mean and is substantially lower than the other four lowest events.The standardized anomalies of lake area during this interval were substantially lower than during 1608-1632, which ranks second among the 25-year long lowest www.clim-past.net/11/1139/2015/periods (Table 2).The marked lake area increase from 1833-1857 is also a prominent feature revealed by the 25-year averages.The periods 1958-2007 and 1808-1857 are similarly the lowest and highest 50 year periods in this 600 year record (Table 2).These results indicate that lake-area during the second half of the 20th century has been comparatively the lowest, whereas during the 19th century have been the higher in the context of the past 601 years.
Following these results, the reconstruction suggests an increase in the occurrence rate of extreme events in lakearea reduction during the 20th century (Fig. 4).To assess changes in the occurrence rate of extreme lake-area along the reconstruction, we employed the recurrence analysis by Mudelsee et al. (2004).The recurrence interval of extreme events (< 20 % of the mean area) was between 4-7 years in the 15th to 18th centuries and between 7-15 years in the 19th century.It has steadily increased from ca. 6 in the 1930s to 2 in recent years.

Spectral properties and Pacific influence on lake area reconstruction
The spatial correlation pattern between the reconstructed lake area and the global gridded sea surface temperatures (SSTs) were determined for the period 1948-2006 (Fig. 5).
The annual lake area is significantly negatively correlated with SSTs from the tropics to the mid-latitude Pacific coasts of the Americas, resembling an El Niño Southern Oscillation (ENSO) -Pacific Decadal Oscillation (PDO) like pattern (Fig. 5).
Figure 6 shows the BT spectra for the annual lake-area reconstruction over the Vilama-Coruto region and the N3.4 SSTs.Peaks that exceed the 95 % confidence limit in the lake area reconstruction are prominent at 12.9, 5.9-6.3,3.2 and 3 years.Except for the decadal peak of 12.9, all significant peaks coincided with the N3.4 SSTs peaks recorded in Fig. 6b (5.9, 3.6-3.7).At an inter-annual timescale, the lake reconstruction contained similar periodicities to N3.4 SSTs, falling within the preferred frequency band of ENSO (2-7 years).
The coherency spectrum was significant at 2.8, 3.4-3.9years cycles, revealing that the lake area reconstruction captures the inter-annual variability in SST records (Fig. 6c).These results suggest that part of the variance in the lake-area reconstruction from the Vilama-Coruto region reproduces the timing and duration of high frequency SST anomalies.
We used singular spectral analysis (SSA), to isolate the main oscillatory modes of the lake area reconstruction variability (Fig. 7).Several important waveforms were recorded at multi-decadal (26-33 years), decadal (11-15 years) and inter-annual (5-8 and 2-4 years) timescales (Fig. 7).The multi-decadal waveform amplitudes around 1600s were the largest in the series at a time when two severe decade-scale reductions in lake area (1583-1595; 1614-1627) were interrupted by a decade of increased lake area (1596-1613).The decadal oscillation mode showed high amplitude variability during the 15th, 18th and the beginning of the 19th centuries, and was reduced during the 16th and the second half of 19th centuries.The reconstruction contains an important 5-8 year cycle waveform that explains a high percentage (27 %) of the lake-area reconstruction variability.The amplitude of this oscillatory mode was relatively large over the entire reconstruction but decreased from the 1920s to present.The interannual oscillation mode (2-4 years) shows high amplitudes during the 15th, first half of the 16th, second half of the 19th and in particular during the second half of the 20th centuries.
A comparative analysis between the dominant oscillatory modes in the reconstruction with those in the Pacific Decadal Oscillation index (PDO) and in the N3.4 SSTs, highlighted the occurrence of common waveforms between them at different frequencies (Fig. 8).The lake-area reconstruction shares similar oscillations with PDO at multidecadal modes of variability, while by decadal and interannual modes, the reconstruction waveforms are closely related to those recorded at N3.4 SSTs.High correspondence in the inter-annual modes are evident in the amplitudes of the lake reconstruction and N3.4 SSTs waveforms showing that both records have large and small amplitudes during 1870-1890/1960-2008and 1930-1960, respectively;(Fig. 8c).

Discussion and concluding remarks
Our study represents the first high-resolution, tree-ring based reconstruction of variations in lake area over the past 6 centuries in the Andes.Through the analysis of 601 years of water balances of the dry Puna lakes, this study provides a novel insight into the climate variability of the subtropical Andes.The evidence used to develop and validate this reconstruction includes satellite images and century-long tree-ring series of P. tarapacana, a climate-sensitive species from the Altiplano.Consistent with the results obtained by Carilla et al. (2013), the use of satellite-derived lake area records as a measure of water balance, resulted in a useful approach to explore hydrological trends in this region devoid of local meteorological records.
The regression model used to develop the reconstruction explains 62 % of the total variance in the January-December lake area over the 1975-2007 calibration period.This value is comparatively higher to that previously reported for a treering based precipitation reconstruction of the Southern Altiplano (R 2 adj = 0.54; Morales et al., 2012).At the sampling sites, precipitation explains high percentage of the total radial growth variance, but also summer temperatures which increase evapotranspiration and reduced soil water supply, affect tree growth variability (Morales et al., 2004).In the same way, the areal extent and depth of water in lakes respond to changes in precipitation and evaporation integrated over the lakes and their catchment basins (Mason et al., 1994).Both tree-growth and lake area integrate changes in climatic parameters including precipitation, temperature, radiation and wind speed.As a consequence, a closer association is expected between tree growth and lake extents than tree growth and singular parameters such as precipitation or temperature.
Extended periods of small lake areas in our reconstruction, such as those recorded in 1504-1595, 1614-1626, 1790-1818 and 1976-2007, highlight the occurrence of decadal to multi-decadal droughts in the region.Similar, relatively longterm humid conditions were reconstructed during the periods 1441-1503, 1596-1613, 1627-1796, 1819-1934.Concurrent dry and wet periods to those registered in our reconstruction have been documented in a tree-ring based precipitation reconstruction for the South American Altiplano (Morales et al., 2012) and in a gridded multiproxy PDSI (Palmer Drought Severity Index) for the South American subtropics (Boucher et al., 2011).
Wet and cold conditions in NWA-SWB, associated with large areas of the lakes during the period 1627-1796, have also been documented for other climate proxies from northern locations in tropical Andes from Peru, such as the δ 18 O isotope records in the ice core from the Quelccaya ice cap (Thompson et al., 2006) and in the lake lake sediment core from Pumacocha (Bird et al., 2011).However, wet and cold conditions in these records occurred during almost 300 years (1500-1800), whereas in our tree-ring reconstruction, these conditions last 169 years from 1627-1796.In contrast to the Quelccaya and Pumacocha records, our reconstruction shows persistent low lake areas in the 16th century.Differences between these records may reflect distinct climate conditions between sites separate by more than 1000 km in north-south direction along the Andes.Discrepancies between proxies are also related to the individualistic nature of records (treerings, ice cores and lake sediments) capturing environmental changes at different frequency and time lags.In addition, small lake areas were recorded in this relatively wet time period (1627-1796) for the years 1661-62, 1724, 1740, and 1753-1755.These dry years were also identified in other high-resolution precipitation reconstruction based on historical documents from the silver mines from Potosi (southern Bolivia) across the period 1585-1807 (Gioda and Prieto, 1999).The annual dry events of the historical series and lake area reductions were consistently related at inter-annual variations during the common period (see statistical analysis in the Supplement S7).
The 19th century showed the highest lake area records of the reconstruction, consistent with the most important wet period identified in the precipitation reconstruction from the Altiplano (Morales et al., 2012).Also, other low-resolution environmental proxies such as the peat-accumulating wetland from Tuzgle mountain in the dry Puna, Argentina (Schittek et al., 2015) and packrat middens in the Andean precordillera from the northern Atacama Desert, Chile (Mujica et al., 2015), highlight the large pluvial period of the 19th century.
Analyses of satellite-derived lake area and vegetation index in the Vilama-Coruto region by Carilla et al. (2013), revealed sustained negative trends in water balance and ecosystem productivity for the periods 1985-2009 and 2000-2010, respectively.These environmental changes have been associated with recent climate variability in NWA-SWB (Carilla et al., 2013).Our lake-area reconstruction indicates that the recorded dry condition of the late 20th century is exceptional over the 1407-2007 period.Three of the five 5-year periods with smallest lake areas occurred from 1930 to present, and the driest 25-and 50-year periods occur in the second half of 20th century.
Recurrence of extremely small lake-area events has been a common feature in the Vilama-Coruto region over the past 6 centuries.However, the 20th century shows the highest rate of extreme events, that changes from a recurrence period of 15-years in the second half of the 19th century to every 2years in the late 20th century.
Since the mid 1970s, the Vilama-Coruto lake system recorded an accelerated decrease in area consistent with the driest condition of the past 6 centuries.This shift to arid condition has also been identified in δ 18 O isotope records in the ice core from Quelccaya (Thompson et al., 2006) and sediment core from Pumacocha (Bird et al., 2011).In both records, this pervasive dry period was also unprecedented in the context of the last 600 years.The dry conditions seem to be a regional pattern affecting the southern tropical An-des, and a major cause of the steady shrinking of glaciers from 1930 to 2005 in southern Perú and Bolivia (Ramirez et al., 2001;Francou et al., 2003;Vuille et al., 2008).Glacier retreat is also related to the increasing temperatures during the last decades across the region (Vuille and Bradley, 2000;Urrutia and Vuille, 2009;Rabatel et al., 2013).Annual temperature deviations averaged over the tropical Andes (1 • N-23 • S) indicate a significant warming trend for the period 1939-1998.However, this trend has intensified since the mid 1970s (Vuille and Bradley, 2000) and was highly consistent with the faster rate of global temperatures increase (Jones et al., 1999).
Significant changes in precipitation patterns in tropical and subtropical South America have been identified in instrumental records during the 20th century.A trend analysis  of atmospheric vertical motion and winds across a north-south transect at 65 • W in tropical-subtropical South America shows a significant increase in precipitation and cloudiness in the inner tropics whereas the southern tropics and subtropics become drier and less cloudy (Vuille et al., 2008).These changes in precipitation and cloudiness have been associated with an intensification of the regional Hadley circulation, with a more vigorous vertical ascent, favorable for convective activity, in the tropics, balanced by enhanced subsidence and less cloudiness in the subtropics (see Fig. 6 in Vuille et al., 2008).
The main moisture source for the Altiplano is the easterly influx from the Amazon Basin (Lenters and Cook, 1997;Garreaud et al., 2009).However, wet or dry episodes are related to changes in the mean zonal wind, largely modulated by Pacific Ocean surface temperature (Vuille et al., 2000;Garreaud and Aceituno, 2001;Bradley et al., 2003).Pacific SSTs properties, such as amplitude, frequency and associated teleconnection patterns, showed considerable natural variability from inter-annual to centennial timescales (Li et al., 2013).In this study spectral analyses reveal periodicities and dominant waveforms in the lake area reconstruction at multi-decadal (26-33 years), decadal (11-15 years) and inter-annual (5-8 and 2-4 years) modes of variability.These oscillations were consistent with the spectral properties from the N3.4 SSTs and PDO indices.The amplitudes of the different modes change throughout the reconstruction.A dominant multi-decadal mode shows large-amplitude oscillations around 1600, whereas the decadal mode was more important in the 1500s and from 1700s to the beginning of 1800s.The inter-annual 5-8 year bandwidth was dominant for almost the entire period analyzed, except after 1950, where amplitudes decreased.The inter-annual 2-4 year bandwidth was strong during the second half of the 19th and 20th centuries.The shift from a dominant decadal to an inter-annual mode was also recorded in a PDO index reconstruction based on western North America and southern Central Andes tree-ring chronologies (D'Arrigo et al., 2001;Villalba et al., 2001).These authors present evidence of a more dominant decadal mode from 1700 to middle 1800s with a shift towards more inter-annual-dominated forcing around that time (D'Arrigo et al., 2001;Villalba et al., 2001).
The main oscillatory modes between the lake-area reconstruction and the instrumental PDO and N3.4 SSTs indices are similar though negatively related over the 1900-2007 and 1872-2007 periods, respectively.Decrease in reconstructed lake areas and Pacific-SSTs amplitudes at interannual modes was recorded from 1930-1960, which was coherent with a reported low ENSO activity during the same period (Aceituno and Montecinos, 1993;Torrence and Webster, 1999;Sutton and Hodson, 2003).These results were highly consistent with the reported relationships between the main oscillation modes of the precipitation reconstruction from the Altiplano and the N3.4 SSTs (Morales et al., 2012).
The BT spectral density analysis suggests that the lakearea reconstruction and N3.4 SSTs contain similar interannual periodicities within the preferred frequency band of ENSO (2-7 years).In particular, the strong spectral coherency between both series within the 2-4 year bandwidth indicates that the timing and magnitude of high-frequency anomalies in Pacific SSTs modulate the interannual lakearea variability in the Vilama-Coruto system.During the past 600 years, 16 % of the variance in lake area is associated with the 2-4 year bandwidth (Fig. 7d); however, the largest contribution to the total variance (27 %) arises from cycles at the 5-8 year bandwidths (Fig. 7c).The amplitude of the 5-8 year waveforms decreased from the 1930s to the present.In contrast, an increase in amplitudes of the 2-4 year oscillation is recorded over the same period.
In contrast to high-frequency oscillations associated with ENSO, decadal to multi-decadal modes of variability in the lake-area reconstruction appear to be associated with changes in the PDO index.Previous studies have documented a long-lived pattern of Pacific climate variability in palaeoclimate records from western North America and the Central Andes (D'Arrigo et al., 2001;Villalba et al., 2001;Masiokas et al., 2010).The persistent trend to arid conditions recorded from the mid 1970s to the present in the Vilama-Coruto lake area reconstruction was coincident with a documented shift from the negative to the positive phase of the PDO in 1976 (Mantua et al., 1997).This PDO positive phase, which lasted to the end of the 20th century, coincided with a marked increases in global mean surface temperatures (Trenberth et al., 2013(Trenberth et al., , 2014) ) and has significantly impacted climate and ecosystems across the western coasts of North and South America (D'Arrigo et al., 2001;Villalba et al., 2001;Masiokas et al., 2010).During this period of positive PDO phase, an increased in ENSO was recorded in both instrumental and reconstructed data (Li et al., 2013).
Interactions between PDO and ENSO have been documented for both North and South America (Andreoli and Kayano, 2005).More intense El Niño (La Niña) rainfall anomalies occur when both ENSO and PDO share positive (negative) phases.In contrast weakened anomalies occur when they are in opposite phases.It has been suggested that the PDO switched to a negative phase (La Niña-like pattern) around 1999, coincident with a pause in the sustained increase in global mean surface temperatures (Trenberth et al., 2014).However, we note that the downward trend in Coruto-Vilama lake areas has persisted since 1999.This suggests that, in addition to PDO, different forcings may have contributed to the persistence of dry conditions in our study region for the past few decades.A dominance of the upper-tropospheric westerlies during the second half of the 20th and early 21st centuries, leading to decreased precipitation by blocking the moisture easterly flow from the Amazon basin (Neukom et al., 2015), together with a continuous warming trend of the upper troposphere over the Central Andes (Vuille et al., 2015) may thus explain much of the recent dramatic lake area decrease.Furthermore, an anomalous ENSO activity in the late 20th century (Li et al., 2013) would even strengthen the negative relationship between westerlies and precipitation, with the consequent intensification of drier conditions.
In the present study, common cycles and periodicities were identified between lake area and ENSO (at interannual/decadal oscillation modes) and PDO (at multidecadal oscillation mode).These results allow to identify ENSO and PDO frequency, amplitudes and associated teleconnection patterns.These findings were consistent with those recorded for the ENSO-related precipitation reconstruction from the Altiplano (Morales et al., 2012), while no previous analyses were made about the PDO influence on the reconstructed precipitation.Since most palaeoreconstructions of ENSO are devoid of tropical/subtropical records (Li et al., 2013), our lake-area reconstruction provides valuable information about ENSO and also PDO evolution during the past 600 years and should be considered a key input in further high-resolution reconstructions of these hemispheric forcing.Given that ENSO and PDO variability provided major factors affecting hydrological patterns in the high-elevation lake systems from NWA-SWB, the longer time frames provide here may also complement those studies discussing the regional manifestations and long-term variability of these large-scale ocean-atmosphere features (Stahle et al., 1998;Evans et al., 2001;D'Arrigo et al., 2001D'Arrigo et al., , 2005;;Villalba et al., 2001;Cobb et al., 2003;Gergis and Fowler, 2009;Li et al., 2013).
At local and regional scales, the lake-area reconstruction should be used as input in water system models to refine the reliability of water-balance assessments and ecosystem productivity in this arid region with persistent hydrological deficits.This study is especially useful in understanding the temporal variability of water availability in the Vilama-Coruto lake system, where the water resource is fundamental for biodiversity conservation and socioeconomic activities such as grazing, mining and tourism.Global and regional climate models for the Altiplano and Puna (Urrutia and Vuille, 2009;Minvielle and Garreaud, 2011;IPCC, 2013;Neukom et al., 2015) project a marked reduction in precipitation to the end of the 21st century, exacerbating presently dry conditions.If the negative trend in precipitation documented in recent decades continues, it is likely that the ranges of natural variability of these high-altitude ecosystems will be exceeded.Therefore, quantitative information on past hydrological changes is required to provide reliable information to better adaptation to future drier conditions in the Altiplano.Further improvements in the existing lake area and precipitation reconstructions could be achieved by expanding the temporal and spatial coverage of the tree-ring network in the Puna and Altiplano from Argentina, Bolivia, Chile and Peru.Similar reconstruction exercises could be performed for neighboring high-altitude lake systems from the Altiplano and northernmost tropical Andes.
The Supplement related to this article is available online at doi:10.5194/cp-11-1139-2015-supplement.

Figure 1 .
Figure 1.Map showing the Vilama-Coruto lake system in the northwestern Argentina and southwestern Bolivia.Green stars represent the selected tree-ring chronologies used in this study.

Figure 2 .
Figure 2. (a) Observed (green line) and tree-ring predicted (blue line) annual lake area (January-December) variations in the NWA-SWB.Annual lake area expressed as a Z score of the satellitederived 1975-2007 lake area mean.Calibration and verification statistics: the Pearson correlation coefficient (r) between observed and reconstructed values, the coefficient of determination adjusted for the degrees of freedom in the model (R 2 adj ) over the calibration period, the F value of regression (the ratio between the explained and the unexplained variability in the model), and the reduction of error (RE).(b) The linear trend of regression residuals (slope: blue line) and the Durbin-Watson (D-W) statistics used to test for firstorder autocorrelation of the regression residuals are shown.

2. 6
Spatial and spectral analyses of the reconstructed lake area and climate forcing The influence of Pacific sea surface temperatures on lake area fluctuations in the Vilama-Coruto area was estimated by examining the spatial correlation coefficient pattern between Clim.Past, 11, 1139-1152, 2015 www.clim-past.net/11/1139/2015/

Figure 3 .
Figure 3. Annual (January-December) Vilama-Coruto lake area reconstruction for the period AD 1407-2007.Annual lake area expressed as Z score of the 1975-2007 lake area mean.The 35-year smoothing-spline curve highlights the multi-decadal variability.Significant (95 % cl) regime shifts (blue arrows) and the mean of periods (orange horizontal line) detected by the Rodionov (2004) method (window length = 25 years).Dates of the regime shifts are shown in the figure.Uncertainties of the reconstruction are shown by the light green band (±1 RMSE).

Figure 4 .
Figure 4. Changing probability of extreme low values (< 20th percentile; blue line) of the Vilama-Coruto lake area reconstruction from 1407-2007.A kernel smoothing method was used with a bandwidth of 50 years.The shaded represent 95 % confidence interval based on 1000 bootstrap simulations.

Figure 5 .
Figure 5. Spatial correlation patterns over the period 1948-2006 between the 2.5 × 2.5 gridded monthly averaged January-December sea surface temperatures (SSTs) and the reconstructed January-December lake area for the Vilama-Coruto region.Spatial correlations were obtained from the National Oceanic and Atmopheric Administration website (http://www.esrl.noaa.gov/psd/data/correlation/).The reconstructed lake area region is indicated by the red square.

Figure 6 .
Figure 6.Blackman-Tukey power spectra analysis of the (a) reconstructed lake area (1407-2007) and (b) the sea surface temperatures from the N3.4 Pacific sector (N3.4-SSTs; 1872-2007).(c) Coherency spectrum of N3.4-SSTs and reconstructed lake area in the Vilama-Coruto region, estimated over the common period 1872-2007.Short dashed and dotted lines represent the 0.05 probability levels and the red noise band, respectively.Dominant periodicities are indicated in each panels.

Figure 7 .
Figure 7. (a) Multi-decadal, (b) decadal and (c, d) inter-annual main oscillation modes of variability extracted by Singular Spectral Analysis (SSA) of reconstructed lake area from the 1407-2007 period.The frequencies for each SSA are indicated in years with the corresponding explained variance in percentage (%).

Figure 8 .
Figure 8. Comparisons between the waveforms of the lake area reconstruction (red line), PDO (blue line) and the N3.4-SSTs (green line) extracted by Singular Spectrum Analysis (SSA) over the (a) 1900-2007 (PDO) and (b, c) 1872-2007 (N3.4-SSTs) intervals.The periodicities and percentage of variance explained by each frequency are indicated in parentheses.The correlation coefficients between the two series are shown at the right corner.The PDO and N3.4-SSTs waveforms are shown inverse to facilitate the comparison between records (right axis).

Table 1 .
Geographical location and statistics of the Polylepis tarapacana ring width chronologies.Mean RBar and mean EPS were averaged over chronology periods including > 10 tree-ring series.

Table 2 .
Lowest and highest non-overlapping averages of the reconstructed (1407-2007) lake area fluctuations for the Vilama-Coruto lake system.Annual lake area expressed as Z score of the 1975-2007 lake area mean.Ranks 1-5 correspond to the five most extreme reconstructed lake area years or set of consecutive years.Rank 1 represents the most important minimum (maximum) lake area record.