STUDY AREA

The study was conducted in the extreme northwest of Argentina (Rinconada, Jujuy) and extreme southwest of Bolivia (Sud Lipez, Potosí), between 21°41′-23°17′S and 65°10′-67°26′W (Fig. 1). All the area is located above 4400 m, and covers approximately 9300 km², including ca. 30 lakes ranging in size from 1 to 40 km². Vegetation is dominated by Festuca spp. grasslands, several small shrubby species (Parastrephia, Adesmia), and cushion plants (Azorella compacta). The only treelet species that occurs in the area, Polylepis tarapacana, can be as tall as 2.5 m (Kessler, 1995). The most productive areas in the landscape are the peatbogs or “vegas,” wet zones dominated by cespitous Juncaceae and Cyperaceae, and Poaceae of the genus Deyeuxia.

There are no permanent human settlements, but temporary “puestos” serve as logistic bases for summer livestock grazing (llamas and sheep). Since the area represents important habitats for vertebrate species with high conservation value such as Andean flamingos (Phoenicoparrus jamesi, P. andinus), vicuñas (Vicugna vicugna), Andean rheas (Pterocnemia pennata), and high Andes wildcats (Leopardus pajeros, L. jacobitus), and there is an incipient tourist development, both the Argentinean and Bolivian sectors are under protected status: Reserva Provincial Altoandina de la Chinchilla in Argentina and Reserva Nacional Eduardo Avaroa in Bolivia.

CLIMATE

The study area (Fig. 1, part B) is cold and dry, with mean annual temperatures lower than 5 °C and an annual precipitation lower than 100 mm ( http://www.worldclim.org). Precipitation over the southern tropical Andes is characterized by a marked seasonality. December to February precipitation accounts for 70–90% of the total annual precipitation (Garreaud et al., 2003; Vuille and Keimig, 2004). Across the Altiplano, precipitation decreases from 600–1000 mm in the northeastern Cordillera to less than 100 mm in the southwestern mountains adjacent to the Atacama Desert (Vuille et al., 2003; Vuille and Keimig, 2004).

Summer precipitation is related to local convection, which in turn is associated with variations in the mean zonal wind in the upper troposphere, and with an easterly zonal flow favoring wet conditions (Garreaud et al., 2003). Rainfall variability appears to be related to the position and intensity of the Bolivian High, an upper-level pressure cell that develops over the central Andes during summer due to the intense convection over the Amazon basin (Lenters and Cook, 1997). Wet episodes are related to a southward displacement of the Bolivian High, the consequent reduction in meridional baroclinicity at subtropical latitudes, and the weakening of the westerly flow over the central Andes promoting the ingression of wet air masses from the east (Garreaud and Aceituno, 2001; Vuille and Keimig, 2004).

TABLE 1

Lakes measured in this study (with codes, see Fig. 1); minimum and maximum lake sizes; coefficient of determination (r²); slope of size trend between 1985 and 2009; % of total lake size change in relation to initial value (2009–1985 size); and Pearson's correlation (r) with regional chronology ring width raw data.

Interannual changes in the mean zonal flow over the southern tropical Andes are modulated by sea surface temperatures (SST) over the tropical Pacific (Vuille et al., 2000; Garreaud and Aceituno, 2001; Bradley et al., 2003). Positive anomalies in SST over the tropical Pacific of the South American coast (El Niño events) are associated with an intensification of the westerly flow and below-average precipitation over the southern tropical Andes. In contrast, rainfall is more abundant during negative SST in the tropical Pacific (La Niña events), concurrent with enhanced easterly flow bringing humid air masses from the east over the region (Aceituno, 1988; Lenters and Cook, 1999; Vuille, 1999; Vuille et al., 2000; Garreaud et al., 2003).

Instrumental records from La Quiaca, the closest instrumental weather station, shows a slight increase in precipitation during the 20th century (r²= 0.1, p < 0.05) which has become more abrupt since 1970 (r²= 0.4, p < 0.001). However, La Quiaca is more than 1000 m lower in elevation than the study area, and 150 km east, thus is likely to be partially influenced by different general atmospheric circulation patterns. The same trend (r²= 0.1, p < 0.01) is evidenced in the average of 11 weather stations of the region during the common period (1903–2007).

LAKE FLUCTUATIONS AND RELATIONSHIPS WITH REGIONAL CLIMATE

To quantify interannual fluctuations in lake sizes (objective 1), we used Landsat TM images (30 m pixel resolution), path 232 and row 076, from 1985 to 2009. We selected nine lakes located between 4400 and 4600 m a.s.l. (Fig. 1, part B; Tables 1 and 2). To describe lake fluctuations we report area (km²), relative size (annual area/maximum area during the period 1985–2009), and area differences between two consecutive years relative to the initial area ([A_x - A_x-1]/A_x-1). To determine the area of the lake in each particular image, we classified images using the non-parametric method Support Vector Machine (SVM; Hsu et al., 2007;  see online Appendix (Lake Fluctuations.PDF)).

TABLE 2

Pearson correlations between average of relative lake size for 1985–2009 period. For explanation of lake codes, see Table 1.

We assigned one size data for each of the 24 years between 1985 and 2009, the only exception being 1986 for which there was no good quality image available. For years during which we had more than one good quality available image, we used the average of the different lake size measures. Available images for each year varied with regard to their acquisition date, but we found that intra-annual variation was on the order of 15% of relative range in comparison to more than 100% of interannual variation, so we assumed that the seasonal variation was negligible. To further validate this assumption, intra-annual fluctuation was not significantly different (repeated measures ANOVA: F_(2,90) = 2.21; p = 0.11) between lake size measures of different months. Dates of images used for this analysis are reported in the online Appendix.

TABLE 3

Spearman correlation between interannual differences in average annual lake size and interannual variations in precipitation records (R) for 11 meteorological stations located along Bolivian and Chilean Altiplano (see Fig. 1, part D, for specific locations), regional precipitation average (R), and de Martonne's Aridity Index (AI performed for La Quiaca weather station). Countries codes: Bo—Bolivia, Ch—Chile, Ar—Argentina. For explanation of lake codes, see Table 1.

To assess the relationship between instrumental records of precipitation variability and changes in lake size (objective 1), we compiled the available high-altitude precipitation records from meteorological stations relatively close to our study area or at similar elevation although located as far as 300 km north (Fig. 1, part C). Monthly precipitation records were obtained from the Servicio Meteorológico Nacional in Argentina (SMN), the Servicio Nacional de Meteorología e Hidrología in Bolivia (SENAMHI), and the Dirección General de Aguas in Chile (DGA). The 11 weather stations used are located between 17°S and 22°S, and between 3500 and 4600 m a.s.l. (Fig. 1, part C; Table 3). We combined the 11 records into one regional precipitation record. To avoid overweighting weather stations with higher rainfall on the regional mean, each precipitation record was standardized ([year value — interval average]/standard deviation) with respect to the common interval 1983–2000 period, and then precipitation anomalies were averaged. We considered a 12-month cycle of lake size and considered annual precipitation from August to July which is more consistent with biological cycles in the region than the calendar year.

To assess the relationship between lake levels and water balance, we calculated de Martonne's aridity index, defined as precipitation (mm)/temperature (°C) + 10(deMartonne, 1926). This aridity index, based on precipitation and temperature records, was only calculated for La Quiaca weather station (ca. 150 km to the east, at 1000 m lower elevation), the only one with a reliable long-term temperature record. We also assessed the relationship between La Quiaca temperature and changes in lake size.

SATELLITE-DERIVED INDICES OF PLANT PRODUCTIVITY

We used Moderate Resolution Imaging Spectroradiometer (MODIS) images from 2000 to 2010, georeferenced and reprojected to UTM 19 south (DATUM WGS84), in order to describe temporal patterns of vegetation productivity using EVI (enhanced vegetation index; Huete et al., 2002). These patterns were based on a 16-day composite (23 images/yr) of 250 × 250 m spatial resolution with rescaled values ranging from 0 to 1. EVI enhances the vegetation signal with improved sensitivity through a decoupling of the background signal and a reduction of atmosphere influences (Huete et al., 2002). EVI temporal series obtained from MODIS were analyzed using TIMESAT software (Jönsson and Eklundh, 2004), which quantifies phenological signals found in temporal series of satellite image data. TIMESAT adjusts local functions for each temporal series point and combines these functions in a model of phenological patterns. Based on these functions, TIMESAT provides statistical descriptors of the seasonal pattern of the analyzed variable (EVI in this case) through the year. For this study, we selected four phenological variables calculated the following for the 9300 km² study area (Fig. A1 in online Appendix). (1) Growing season duration (GS; in days between GS starting and ending). For the first year in the series, GS starts in mid-2000, and ends in mid-2001. (2) Maximum EVI value (MV). (3) Growing season amplitude (AMP, the difference between MV and base value of EVI that is the average of the EVI minimum values; beginning and ending of GS). (4) Seasonal total integral (TI), an index of absorbed photosynthetically active energy accumulated in one growing season (Running et al., 2004).

To explore associations between lake fluctuations and plant productivity (objective 2a), we used Pearson's correlation coefficients.

DENDROCHRONOLOGICAL ANALYSIS OF POLYLEPIS TARAPACANA

In this study, three tree-ring chronologies from P. tarapacana were developed by incorporating two new chronologies from Ramada and Cerro Negro sites as well as updating and extending back in time the previously published chronology from Granada Volcano (1659–1999) (Morales et al., 2004) (Fig. 1, part A; Table A1 in online Appendix). All three sites are in the study area, between 4400 and 4800 m a.s.l., and share similar climatic characteristics (Fig. 1). Cross sections were chainsaw extracted from branches of living and dead trees, which remain well preserved for several centuries due to the dry and cold climate. Samples were processed following standard dendrochronological techniques (Stokes and Smiley, 1968). For dating purposes, we followed Schulman's convention (1956), which assigns to each tree ring the date of the year in which radial growth started. Annual rings were measured under a binocular stereoscopic microscope with a Bannister dendrometer measuring machine (0.001 mm precision) connected to a computer (Robinson and Evans, 1980). The quality of the tree-ring chronology was tested by mean sensitivity (MS), mean correlation between series (RBAR), and expected population signal (EPS) statistics (described in online Appendix) generated by the computer program COFECHA (Holmes, 1983). Chronology statistics ranges met generally accepted reliability standards; mean intercorrelations: 0.51–0.57; MS: 0.26–0.31; RBAR: 0.26–0.32; and mean EPS: 0.71–0.87 (Morales et al., 2012).

The three P. tarapacana ring-width chronologies showed a similar interannual growth pattern across the study area, which allowed composition of a regional tree-ring chronology, merging in a single record the 93 tree-ring series from the three sites. A relative high series intercorrelation (r = 0.56) indicates the strong internal coherence between the three site chronologies (Morales et al., 2012) computed over the common period 1850–1999 (>6 samples in each site).

The ring-width measurements of the regional chronology were standardized to conserve the low frequency signal and to remove components of the variability likely not related to climate. Details on standardization methods and the statistics of the tree-ring chronologies are given in the online Appendix. The resulting regional standardized tree-ring chronology as well as the raw chronology (non-standardized ring width measurements) were used to explore correlations with lake size and with productivity index fluctuations during the recent past (objective 2) in order to evaluate the potential of P. tarapacana growth as a proxy of both factors. In addition, the 258-yr chronology was used to contextualize recent ecosystems changes into long-term ranges of variability (objective 3).

In general, when normality and homoscedasticity assumptions were not met we used non-parametric Spearman correlation. We report independent “p” values for each pair of variables analyzed without adjusting for multiple comparisons (e.g. Bonferroni correction) or spatial autocorrelations. Although we acknowledge this fails to meet assumptions of statistical independence and degrees of freedom, we preferred it to provide readers with more transparent quantitative information for each independent analysis.

LAKE SIZE FLUCTUATIONS AND RELATIONSHIP WITH INSTRUMENTAL RECORDS

Interannual size fluctuations of all lakes for the entire study area were highly positively correlated, suggesting a strong regional climatic control. Only Coruto lake (Bolivia) showed comparatively low correlation with some Argentinean lakes (between 0.3 and 0.4), and 22 cross-correlations were larger than 0.6 (Table 2). Lake sizes showed a negative trend between 1985 and 2009 (Fig. 2, parts A–C; Table 1), suggesting consistently drying conditions. Particularly dry years include 1996, 2004, 2007, 2008, and 2009. In 1996 lakes consistently decreased in size; Cerro Negro lake shrank to 2.4% of its maximum size, and in 2009 five lakes showed a significant decrease in size, with Isla Grande lake becoming completely dry. In contrast, in 1987 six of the nine lakes recorded their maximumsize. The larger lakes, Vilama and Coruto, also showed a decreasing trend (Fig. 2, part D), but since the interannual fluctuations were higher (variances of 15.7 and 5.1 km², respectively), the correlation coefficients were significantly lower (Table 1).

FIGURE 2.

Trends in lake size, relative to maximum, for the 1985–2009 period. (A, B) Different groups of smaller lakes, (C) average of the six smaller lakes, and (D) larger lakes. Gray lines represent the linear regression.

Fluctuations in lake size were correlated with the instrumental precipitation records from 11 high-elevation weather stations, and with de Martonne's aridity index calculated from La Quiaca instrumental records; values of Spearman's correlation coefficient ranged between 0.34 and 0.7 (Table 3). In particular, strong correlations were found between the annual differences in the mean of the six smaller lakes size and interannual variations in precipitation of the previous hydrological year, with higher coefficients (0.51–0.66) in stations closer to the study area. Weather stations located below 4000 m a.s.l. presented weaker correlations with lakes than stations above this elevation (Fig. 1, part D; Table 3). With only two exceptions, Chojllas and Coruto, lake sizes were correlated with the regional average rainfall of the 11 weather stations with r values of 0.41–0.60 (Table 3). The different lakes were correlated with de Martonne's aridity index with r values between 0.34 and 0.57 (Table 3). The joint fluctuation of the differences in the nine lakes was correlated with precipitation of all meteorological records, with regional precipitation, and with de Martonne's aridity index of La Quiaca (Fig. 3, parts A and B). Lake size fluctuations were also negatively correlated with La Quiaca temperatures of the previous year (for six lakes with r²ranging from 0.14 to 0.38, and the average of the six smaller lakes; Fig. 3, part C) and of the current year (for eight lakes with r²ranging from 0.16 to 0.32, and the average of the nine lakes, r²= 0.33, p = 0.003).

LAKES AND VEGETATION PRODUCTIVITY PATTERNS

Between 2000 and 2010, TIMESAT output for mean value (MV) and amplitude (AMP) was lowest in 2008 and highest in 2000 (Fig. 4, part A). Total integral (TI) presented the lowest values in 2000 and the highest in 2002, and GS duration ranged from 188.4 days (in 2000) to 229.8 days (in 2002) (Fig. 4, part B). The four EVI variables showed a decreasing trend in the 10 years of analysis. The decreases in MV and AMP were consistent throughout the period, while the TI and GS durations were not as a consequence of very low values in the first year of the series.

Fluctuations in the size of the six smallest lakes were associated with EVI variables; in particular, average size of the six smaller lakes were positively associated to MV (Fig. 4, part C) and with AMP (Fig. 4, part D) of the current and previous years (r²= 0.89 and 0.85, respectively; p < 0.01). Instead, Vilama, Coruto (the two largest lakes), and Chojllas lakes were not significantly related to variables of EVI seasonality. Regional precipitation was correlated to MV (r² = 0.75, p < 0.01) and to AMP (r²= 0.56, p < 0.05) of the current year, and to the aridity index (MV: r² = 0.73, p < 0.01; and AMP: r² = 0.60, p < 0.05).

Between 1985 and 2008, average annual lake size significantly correlated with P. tarapacana annual radial growth. Regional raw tree ring width chronology presented correlation coefficients higher than 0.44 with the size of the six smaller lakes (Table 1) and with the average sizes of these lakes (Fig. 5, part A), as both variables showed a decreasing trend during the last decades. The standardized tree ring chronology was correlated with interannual differences in the average lake sizes (Fig. 5, part B). Despite the short period with common data (n = 10 years), tree radial growth (raw and standard chronologies) also showed positive correlations with TI (r = 0.7, p < 0.05; Fig. 5, part C) and GS duration (r = 0.64, p = 0.06) of the current year.

ECOSYSTEM PATTERNS AND LONG-TERM VARIABILITY

A record of 258 years of growth variability was provided by the regional P. tarapacana ring width chronology (Fig. 6). This chronology covers the period 1553–2008, and a reliable number of tree ring series (N > 15 trees) have been recorded since 1750. The regional chronology showed a clear pattern of below-average growth during the first 70 years (1750–1820), followed by a peak of above-average growth between 1825 and 1850. A second period of relatively high growth occurred between 1945 and 1960, followed by a decrease during the last five decades of the series (1960–2008). Multidecadal decreasing trend in radial growth of P. tarapacana during the last decades is coincident with lake area reductions and decreasing water balance reflected in lake size. According to the tree ring chronology, these relative low-growth conditions associated with a dry environment with low productivity in the current period have not occurred in the area since 1825, and are becoming similar to the 1750–1825 period.

FIGURE 3.

Relationship between lake size and meteorological data of the previous year. (A) Regional precipitation (averaged anomalies from the 11 meteorological stations with respect to the common period 1985–2009), (B) de Martonne's Aridity Index (based on La Quiaca station, 1985–2001 period) (size average of all lakes), and (C) mean temperature of La Quiaca (size average of the six smaller lakes).

FIGURE 4.

Temporal trends of TI MESAT parameters derived from the enhanced vegetation index (EVI) product from MODIS; (A) mean value (MV) and amplitude (AMP) (p < 0.05), and (B) total integral (TI) and growing season (GS) duration. Relationship between lake sizes (average of the six smaller lakes) and (C) MV and (D) AMP for the 2000–2009 period (p < 0.01).

FIGURE 5.

Relationship between P. tarapacana radial growth and indices of ecosystem function. (A) raw tree ring chronology vs. average lake size, (B) standardized tree ring chronology vs. differences in average of six smaller lakes sizes, and (C) raw tree ring chronology vs. EVI total integral.