Introduction

Nuñoa, Peru located in the Department of Puno, is colloquially known as the alpaca capital of the world. There, the livelihood of the predominantly Quechua-speaking population is largely based on herding and small-scale agriculture. In a semi-arid region with a pronounced dry season and frequent drought, high elevation wetlands, with their year-round vegetation, provide essential foraging grounds. These Andean wetlands, locally termed bofedales, are sources of sustenance and income to local pastoralists. They also support many endemic species of high ecological value (ALT-PNUD 2001) and have the potential to sequester large amounts of carbon (Earle et al. 2003). Similar to wetlands elsewhere, they play an important role in storing water and regulating flow. Bofedales are sustained by a combination of precipitation and snow and glacier melt contributions via upslope drainage or overland flow. They are therefore extremely responsive to changes in the water balance (Squeo et al. 2006). Heavily vegetated, the growth and species composition in bofedales are a function of water availability, soil, and climate. Thus, variations in climate as well as differing land use and herding management strategies can result in significant changes in water quantity, water quality, plant diversity and plant cover. Environmental change and its determinants is therefore an issue of consequence for the numerous herders who rely on bofedales in the Nuñoa watershed and elsewhere in the Central Andes.

Mountain regions are particularly sensitive to changes in climate (Diaz et al. 2003). Temperature in the Andean region has increased by approximately 0.1 °C per decade for nearly 70 years (Vuille et al. 2008). Andean tropical glaciers are receding at accelerating rates (Buffen et al. 2009; Jordan et al. 2005) and some are disappearing. Furthermore, recent research suggests shifts in the precipitation cycle that include later onset of the rainy season (Thibeault et al. 2010) and therefore a longer dry season in some areas. In addition to changes in climate, which may directly alter vegetation and hydrological resources, a myriad of social and economic changes are occurring in the Andes. These changes affect the governance and management of natural resources, including bofedales.

The objective of this study is to investigate multi-decadal (1985–2010) vegetation trends during the dry season (June through August) in the Nuñoa watershed (Fig. 1). Based on the size characteristics of the wetland systems of interest and the scale of study we used Landsat derived normalized difference vegetation index (NDVI). Landsat derived NDVI has been used elsewhere in the Andes to investigate changes in vegetation cover both over time (Yager 2009) and in response to drought (Washington-Allen et al. 1998), and specifically to map the distribution of bofedales (Otto et al. 2011; Izquierdo et al. 2015). Bofedales have the highest vegetation cover relative to the surrounding montane ecosystem, particularly during the dry season. This density of photosynthetic vegetation has high reflectance in the near infra-red making these wetlands easily distinguishable using the NDVI.

Fig. 1
figure 1

Elevation map of the Nuñoa Watershed and showing the three districts that form the political boundaries in the watershed. The inset displays the location of the watershed within the Department of Puno, Peru

Full size image

The specific research goals of this study are to:

  1. 1.

    Delineate bofedales and calculate the area of this natural resource in the Nuñoa watershed. The calculation is based on the 1985–2010 dry-season mean NDVI for each pixel in the watershed.

  2. 2.

    Determine the trends in NDVI in the watershed from 1985 to 2010 particularly for bofedales. To determine NDVI trends we perform a multiple linear regression for each pixel in the watershed (3,070,160 regressions) with cell specific annual dry-season NDVI as the response variable (n = 20) and time (year) as one predictor variable. The Altiplano region experiences a lengthy dry season and high interannual variability in precipitation and temperature. As a result we include three other predictor (nuisance) variables in the model. Julian day accounts for intraseasonality. A regional precipitation index and regional temperature index account for interannual variability in climate.

Forkel et al. (2013) show that interannual variability in a time series has the ability to decrease the performance of trend slope estimates. Indeed, trend analyses of time series must include methods that capture the trend of interest and separate out other components in the time series (such as seasonality and noise) (Forkel et al. 2013; Eastman et al. 2009). Owing to the limited number of cloud-free Landsat images available for the geographic region of interest, the current study utilizes 20 scenes that span the time period from 1985 through 2010. The NDVI series does not contain seasonal and cyclic components present in more complete NDVI datasets however, the trend of interest may be influenced by the date of image acquisition and the extreme year to year variability in climate. Studies using more complete time series are able to estimate and remove these effects and the trend analysis is performed on the decomposed time series (Brandt et al. 2014; Eckert et al. 2015). These methods are not applicable where data are limited, such as the current study (only monthly data was available to us). In these cases, the multiple regression approach applied here is able to account for other effects (e.g., intraseasonal and interannual effects) on the longer term trend in NDVI. Similar approaches have been used in other studies investigating trends in NDVI over time and spatially within a watershed. For example, Alatorre et al. (2011) and Alatorre et al. (2015) perform regression analyses with NDVI as the response variable and a set of covariates akin to those used in this study, including climate indices, Julian day, and year of image acquisition. The authors discuss the importance of including such covariates in order to isolate the trend of interest. Furthermore, the addition of year as a covariate allows for the attribution of vegetation trends to non-climatic factors (Alatorre et al. 2011, 2015; Zhang and Wang 2015).

Site description

The study area is the 2763 km2 Nuñoa watershed in the Peruvian Altiplano (Fig. 1). The climate is semi-arid with a seasonal cycle of precipitation characterized by dry months (June through August) when there is often no precipitation to the wet months of December through March when precipitation may be 200 mm per month. Annual precipitation for the region averages 760 mm but year to year variability is high. Relatively frequent drought and heavy rainfall periods are common and are associated with El Niño and La Niña events, respectively (Garreaud and Aceituno 2001). The dry season months of June and July are typically the coldest of the year averaging 4 °C while the months of December through March are warmer with temperature averages of 10 °C (ATDR-Ramis 2008). However, the region experiences greater variation in diurnal temperatures than in annual temperatures and this is most pronounced in the dry season when the daily temperature often varies between −12 and 12 °C.

The watershed ranges in elevation from 3726 to 5520-m above sea level (Fig. 1) and ecologically is characterized as predominately humid puna. Humid puna is one of four phytogeographic regions distinguished in the Central Andes (according to the Rivas-Martinez model: Josse et al. 2009) and occurs at elevations between 3800 and 5000 m. Today, only sparse forests of Polylepis remain and the dominant vegetation is a tussock or bunch grass known as Stipa (e.g., Stipa ichu). The landscape includes rolling hills in the puna ecosystem and lichen and moss covered rock in anticline formations at higher elevations. There are several snow covered mountains and along the northern border of the watershed.

Bofedales are common in the humid puna and range in shape and size according to geomorphic setting (bordering streams, at the toe of a slope, or on relatively flat ground) and can cover several square kilometers (Warner et al. 2008). These wetlands are characterized by saturated soils, a high percent of vegetation cover, and significant wet and dry biomass. As a result, they provide essential forage during the dry season and in drought and are used more intensely during these times. In the Department of Puno, vegetation coverage in wetlands may be as high as 85–98 % (ALT-PNUD 2001). Vegetation includes grasses (Poaceae) such as Deyeuxia rigescens and Festuca d., sedges (Cyperaceae) including the species Carex sp. and Scirpus rigidus, flowers in the Asteraceae family such as Hypochoeris, and cushion plants such as Distichia muscoides (Juncaceae family). A large percentage of the vegetation in bofedales is consumed by alpacas.

Data and methods

Landsat 5 thematic mapper (TM) image acquisition and preparation

Landsat 5-TM (30 m resolution) satellite imagery was obtained through USGS/EROS (Earth Resources Observation and Science Center) and the Brazilian National Institute for Space Research (INPE). Images available for the months of June through August from 1985 through 2010 were reviewed for quality, including cloud coverage. Relatively cloud free, dry season-imagery exists for 20 of the 26 years of interest. The images that were processed and included in the NDVI analysis appear in Table 1 along with image properties and associated climate variables.

Table 1 List of images used in the analysis along with image properties and the associated temperature and precipitation values
Full size table

Bands 1–5 and 7 were stacked using ENVI 4.7 software. Georeferencing was performed on images in their original datum using ENVI software package (Version 4.5, ITT Visual Information Solution) and an ETM scene as the base image. In order to obtain the smallest root square mean error (RMS) possible (with a goal of 0.5 pixel error or less) images were first warped using a linear stretch. If, after warping and upon inspection, a location was off by more than one pixel then the image was warped using a second degree polynomial method.

Atmospheric correction (haze reduction) in the form of dark object subtraction (DOS) was performed on each scene using open lake water as the dark object (Ahern et al. 1977). Images were converted from raw digital numbers to represent true reflectance using rescaled gain and bias factors provided in Chander et al. (2009). Converting from DN to radiance, then to reflectance, and performing dark pixel subtraction was done following Eq. (1) below and calculated in ENVI using the Spectral Math tool (Version 4.5, ITT Visual Information Solution).

$$Ref\,L = \left[ {\pi \times d^{2} /\cos \left( {zenith\,\theta } \right)} \right] \times \left[ {\left( {DN{-}DkPixel} \right) \times G_{rescale} + B_{rescale} } \right] /Esun$$
(1)

where RefL = Top of the atmosphere reflectance (unitless); π = Mathematical constant (unitless; ~3.14159); d = Earth–Sun distance (astronomical units) (Table 6, Chander et al. 2009); zenith θ = Solar zenith angle (degrees); DN = pixel value (unitless; 0–255); DkPixel = minimum dark pixel value; G rescale  = Band-specific rescaling gain factor [(W/(m2sr µm))/DN] (Table 3, Chander et al. 2009); B rescale  = Band-specific rescaling bias factor [(W/(m2sr µm))] (Table 3, Chander et al. 2009); Esun = Mean exoatmospheric solar irradiance (W/m2 µm) (Table 3, Chander et al. 2009).

NDVI calculation

The NDVI is one of several remote sensing vegetation indices used to represent vegetation characteristics (e.g., leaf area index, plant cover) in an area (pixel) (Rouse et al. 1973; Tucker et al. 1985). The index is a single value derived from the ratio of the difference in spectral bands that are sensitive to the spectral characteristics of leaf tissue (infrared reflectance is sensitive to plant cells and water content while reflectance in the visible red wavelengths is sensitive to chlorophyll) to the sum of these bands. Specifically NDVI is calculated as (near infrared − red)/(near infrared + red). Using Landsat TM images, this results in the following:

$$NDVI = \left( {Band\,4 - Band\,3} \right) /\left( {Band\,4 + Band\,3} \right)$$

Typically, NDVI values range between −1 and 1. For this study, in order to remove erroneous values and convert to unsigned integers with precision to the hundredths place, values were transformed to an integer scale of 0–20,000 by adding one to the original NDVI value, multiplying by 10,000, and rounding. Original regression model output and calculations use these transformed NDVI units. In these instances, the presentation of the transformed NDVI units will be denoted as NDVIT. Subsequent presentation of NDVI trends and related discussion use conventional NDVI units and will be denoted simply as NDVI.

Post-processing and classification of bofedales

The 20 transformed NDVI Images were then stacked and cropped to a grid containing the Nuñoa watershed. The stacked and cropped image was exported as a band interleaved by line (BIL) file for compatibility with R to perform the regression analyses. Additional analyses and the classification of bofedales were performed using ArcMap 10 (ESRI 2012). For the classification map of wetlands, the average NDVI for the 20 scenes (1985 through 2010) was calculated for each pixel. Based on field visits, GPS data, and Landsat images, three known wetland systems were delineated as regions of interest (ROI). Landsat images were used as a guide to visually inspect mean NDVI values in the three ROI and the surrounding land cover. It was determined that a threshold of 0.3 for mean NDVI clearly distinguishes wetlands from other land cover, including puna vegetation. This is because in the dry season there is high contrast between wetland vegetation and dry grassland vegetation. This conclusion and the threshold value are supported by studies in similar ecosystems in the Andes (Otto et al. 2011; Yager 2009).

An accuracy assessment was performed on the resulting classification map consisting of two categories: wetland and non-wetland. Reference data for the wetlands include 29 GPS ground truth locations collected as part of a Peruvian Government study (ALT-PNUD 2001). Additional reference points for both wetland and non-wetland were obtained from stratified random samples from a digitized land cover classification map [from the 1:100,000 National Map (Instituto Geografico Nacional 1972)]. The sampling scheme follows the recommendation of Congalton (1991) who also suggests a minimum of 50 reference points per class. The reference points were cross-checked in Google Earth (accessed in March 2016; image dates ranged from 7/2010 to 11/2013) and resulted in the identification of 65 wetlands and 79 non-wetland locations for a total of 144 reference points.

Covariate data

Although we were primarily interested in the change in NDVI through time (Year), we had to include as covariates regional precipitation and temperature to account for high interannual variability in the Altiplano climate associated with El Niño and La Niña events. We also included the date of image acquisition (day of the year or Julian day) as a nuisance variable.

The regional precipitation value (Precip) associated with each image represents the sum of precipitation for the preceding months of December, January, February, and March (Table 1; Fig. 2a). These are the wet season months in the Altiplano when collectively, approximately 71 % of precipitation falls (based on data in this study). For the images acquired from 1985 to 2007, the wet-season value is an average across nine meteorological stations. The meteorological stations are currently maintained and operated by Servicio Nacional de Meteorologia e Hidrologia del Peru (SENAMHI-Peru) in the Department of Puno, southern Peru. They are located at the northern extent of the Altiplano from 14°41′ to 15°29′ south and 69°45′ to 70°42′ west and range in elevation from 3820 to 4400 m. Seven monthly values were missing in total from two stations corresponding to three wet seasons. Using backward stepwise regression and the Akaike information criterion (AIC), approximately 2 % of the wet season totals across all nine stations for the 17 years were imputed. The consistency of station records was reviewed with double mass analysis (Searcy and Hardison 1960). Precipitation values for the last 3 years of image analysis (2008 through 2010) were available from one station only (also used in the 1985–2007 analysis). This station is located at 3971 m. The correlation result between this station and the nine station average for the years 1985 through 2007 is r = 0.52, p < 0.05.

Fig. 2
figure 2

Graphs of the explanatory variables: a precipitation (mm); b temperature (°C) and c JDay, in their original values

Full size image

Regional temperature data (Temp) for all years was obtained from the meteorological station located at 3971 m and maintained by SENAMHI-Peru. The temperature value associated with each image is the average for the month preceding the image acquisition date if the scene was acquired within the first 15 days of the month. If the image was acquired after day 15, then the temperature average for the month of acquisition is used (Table 1; Fig. 2b).

Julian day (JDay, included as a nuisance variable), is the date of image acquisition [as the number of days since the end of the previous year (midnight on 12/31)]. For the 20 images used in the analyses it ranges from day 176 (June 25) to day 230 (August 18) (Table 1; Fig. 2c).

The covariate Year is simply the year of image acquisition (Table 1).

Pixel by pixel NDVI regression analysis (1985–2010)

The following ordinary least squares regression model was applied to each pixel (30 m × 30 m) in the watershed using the BIL file containing the 20 stacked NDVI processed images:

$$\begin{aligned} NDVI_{cell} & = \alpha + \, (\beta_{0} \times Year) + (\beta_{1} \times Precip) \\ & \quad + (\beta_{2} \times Temp) + (\beta_{3} \times JDay) + \varepsilon \\ \end{aligned}$$
(2)

where the response variable is NDVI for a given cell, α is the intercept, β i is the slope coefficient for the ith predictor variable, and ε is the residual error. Each predictor variable was centered (C x  = (x − mean(x))) to reduce high correlations between slopes (β 0, β 1, β 2, and β 3) and the intercept (α). In addition, a log base 10 transformation was applied to Temp and Precip values to improve homoscedasticity of residuals. The regression model was individually fitted to each of the 3,070,160 pixels in the watershed and uses 20 years of images and associated NDVI cell (each pixel regression has n = 20). As a result, each cell is assumed to be independent of its neighbors. The model was run in R using the package ‘raster’ (Hijmans and van Etten 2011) which allows for the following output in the format of raster images:

  1. 1.

    Intercept values (α);

  2. 2.

    Slope coefficient values for each of the covariates: Year, Precip, Temp, and JDay;

  3. 3.

    p value (F-statistic) values for each of the covariates;

  4. 4.

    R 2 values for each pixel’s model; and

  5. 5.

    p value (F-statistic) for each pixel’s model.

Utilizing the parameter estimates for each pixel, predicted NDVI cell was calculated for each cell in a given analysis year as were residuals. Boxplots for observed NDVI by year, predicted NDVI by year, and residuals by year are included in “Appendix” (Fig. 8, respectively). Model fit and assumptions were checked by reviewing plots of residuals versus predicted values for patterns (Figure S.1 through S.20 in Supplementary Information). In addition, absolute percent error (100 × |Predicted − Observed|/Observed) was calculated for each cell in the watershed. The maximum absolute percent error and mean absolute percent error for the 20 years of analysis at each cell are presented below together with a brief discussion of the R 2 values for each pixel’s model.

For all analyses, a mask was applied to remove cells within lakes. Lake perimeters were digitized from the National Map covering the region at a scale of 1:100,000 (Instituto Geografico Nacional 1972). Furthermore, with the exception of the reporting of R 2 values and maximum absolute percent error, all other analyses associated with the regression analysis are limited to significant cells. In these instances, cells that were not significant in the regression model at p < 0.05 have been excluded from the analysis. In total 1,588,491 cells in the watershed were not significant at the 5 % level. This represents slightly more than half of the total cell count in the watershed.

Results

Delineation of bofedales

Figure 3 shows the classification map of wetlands in the Nuñoa watershed. The error matrix with accuracy percentages appears in Table 2. The total classification accuracy is 90 % and the κ statistic is 0.79. The results show good agreement and a high level of accuracy. Using the average NDVI value of 0.3 as a threshold, approximately 451 km2 are classified as bofedales. This represents approximately 16 % of the watershed area. The majority of wetland systems are riparian. Their features are narrow and dendritic, paralleling river and streambeds (Fig. 3b). There are a few wetland systems, namely in the northwest and in the south of the watershed that are of the broader type found in shallow depressions (Fig. 3c). NDVI change in the wetlands is discussed following the presentation of the regression results below.

Fig. 3
figure 3

a Classification map displaying wetlands in the Nuñoa watershed delineated using the threshold value of 0.3 mean NDVI. Two examples of wetland types in the watershed include; b a riparian wetland associated with a stream or river and c a broader wetland, in this case associated with an alluvial fan

Full size image
Table 2 Error matrix for wetland and non-wetland categories
Full size table

NDVI regression analysis

In this section we present the results of the regression model (Eq. 2) applied to each cell in the watershed (3,070,160 regressions). First we review the model output for the coefficient of determination and error terms. This is followed by a review of the spatial distribution of each of the covariates and finally, the results of the NDVI change analysis for bofedales.

Model review

A review of the spatial distribution of the coefficient of determination (R 2) values (Fig. 4) reveals many cells in the north and northeast of the watershed with very low values. This region of the watershed contains the highest elevations (see Fig. 1) therefore some low R 2 values may be the result of the mountainous terrain and related slope, aspect, and shadow effects. Snow and ice cover may also reduce the explanatory power of our covariates. A Pearson’s product-moment correlation analysis between elevation and R 2 indicates a low inverse association between the two (Pearson’s r = −0.31, p < 0.05). Note that elevation cannot be put directly into our model because each regression for a single pixel has a single elevation. Very narrow linear features, such as streambeds and smaller wetlands, may also exhibit low R 2 values due to the 30-m pixel resolution.

Fig. 4
figure 4

Coefficient of determination (R 2) values for each cell modeled in the watershed, but with lake areas removed

Full size image

To investigate the magnitude of absolute percent error for all years and the spatial distribution of this error, a map displaying the maximum absolute error for each cell was created (Fig. 5a). While the mean for the watershed is 7.2 % the range of maximum absolute percent error is great and extremely large values appear mostly in the north and northeast of the watershed (similar to R 2). Again, this is presumably due to the topography of the watershed and varying snow and cloud coverage which is difficult to account for. There are also some extreme error values in the southern portion of the watershed (near lake shorelines, for example). These values are associated with the annually changing lake boundaries and also the area of intensive irrigated agricultural use in this region.

Fig. 5
figure 5

Error maps: a maximum absolute error (%) in the 20 years of analysis for each cell (with the exception of cells within lakes); b absolute percent error in the watershed averaged over the 20 scenes and limited to significant cells only and c counts for mean absolute percent error categories

Full size image

Average absolute percent error for cells that are significant at p < 0.05 appears in Fig. 5b. The bar chart in the figure (Fig. 5c) shows that the majority (92 %) of significant cells in the watershed have an average absolute error value of 2 % or less. The watershed average is 1.9 %.

Covariate results

Figure 6 shows the spatial distributions of the slope parameter β i for each of the four explanatory variables. Summary statistics for each covariate and slope parameter along with the corresponding change in NDVI for the given mean and mode appear in Table 3. Minimum and maximum values are the result of truncation of the original data sets to eliminate extreme data values with low counts.

Fig. 6
figure 6

Slope parameters, β, for the four explanatory variables in the regression model: a JDay; b Precip; c Temp and d Year. Only significant cells are included. (Color figure online)

Full size image
Table 3 Summary statistics for the covariates and slope parameters in the regression analysis. The corresponding change in NDVI based on the mean and mode slope parameters is also provided
Full size table

The slope parameter for JDay (Fig. 6a) is the only one to have a negative mean and mode (Table 3). This indicates that as the dry season progresses NDVI values decrease. This relationship is logical if falling NDVI values are the result of plant water stress and one assumes that there is generally a decrease in soil moisture as the dry season progresses. The southern half of the watershed, characterized by flatter terrain and lower elevations, contains the largest percentage of cells with a positive relationship between JDay and NDVI. This may be due to several factors including processes such as infiltration and subsurface flow in the region and irrigation activities.

The slope parameter values for Precip (Fig. 6b) are mostly positive indicative of the relationship between increased precipitation and vegetation. A few cells such as on mountain tops and in the depths of the bofedales exhibit an inverse relationship between the two. The reasoning for this may be that as precipitation falls on mountain tops in the form of snow, NDVI decreases. In the bofedales, increased precipitation may raise the water level to a critical point that may cover some of the vegetation, thereby reducing NDVI.

The slope parameter for Temp (Fig. 6c) exhibits considerable variation in magnitude and sign throughout the watershed as demonstrated by the mode and standard deviation of the slope parameter (Table 3). This may be due to the difference in vegetation response to temperature at various elevations, slopes, and aspects and, to various landforms and vegetation species. Furthermore, the use of a regional temperature variable may not be the most appropriate choice ecologically. However, a review of model performance at the cell level at various locations indicates that while the covariate itself may not be significant in many instances, model performance is enhanced when it is included.

The slope parameter for Year (Fig. 6d) is mostly positive with both mean and mode slightly above zero (Table 3). Since the beta coefficients represent the associated change in NDVIT with respect to a unit increase in Year (holding the other covariates constant) we can calculate the total change in NDVIT over the time period for a given beta value. For example, the mode slope parameter value for the covariate Year is 13 NDVIT/year (Table 3). The change in NDVIT that is associated with this beta value is calculated as 325 (13 NDVIT/year × 25 years). Since NDVIT was scaled by a factor of 10,000, this equates to an increase of 0.03 units on the untransformed NDVI scale. In other words, the change in NDVI from 1985 to 2010 most encountered in the watershed is the modest increase of 0.03 units. We continue to use this method to analyze the increasing and decreasing NDVI trends in the watershed.

NDVI trends

The Nuñoa watershed is 2763 km2. After removing pixels that are well within lake perimeters as well as pixels that are not significant in the model at p < 0.05, the area summarized for change in NDVI over time is 1329 km2 or approximately 48 % of the watershed area. The spatial patterns of positive and negative NDVI trends are shown in Fig. 7a, b, respectively. The trends are in total NDVI units of change from 1985 to 2010. Of the summarized area, 81 % exhibits an increase in NDVI from 1985 through 2010 (Fig. 7a). Another 16 % of the cells exhibit a decrease in NDVI (Fig. 7b) over this time and approximately 3 % of the cells display no change in NDVI (i.e., these cells have slopes <0.0025). Frequencies for the changes in NDVI appear in the histograms (Fig. 7c, d). Cells with an increase in NDVI of 0.1 or greater appear mostly in three locations in the watershed: in the north/northwest along the watershed boundary where there are snow- and few, small ice-capped mountains; in the central section of the watershed in association with the wetland below the alluvial fan; and in the southern portion of the watershed. Cells with decreasing NDVI trends appear as linear features in the north and northeast of the watershed area of the watershed. This is in association with the stream network. However, the majority of cells exhibiting a large magnitude of decreasing NDVI appear in the south of the watershed.

Fig. 7
figure 7

Maps displaying significant a increasing NDVI and b decreasing NDVI from 1985 to 2010. Many, but not all of the wetland systems are marked. Histograms for cells with c increasing NDVI and d decreasing NDVI appear below the maps. (Color figure online)

Full size image

NDVI change in bofedales

Of the 451 km2 classified as bofedales approximately half that area (47 %) was included in the NDVI change analysis. If we calculate the number of cells with a given NDVI trend (limited to model significance of p < 0.05) within bofedales, we find that 68 % (144 km2) of the analyzed wetland area exhibits an increase in NDVI. An additional 29 % (61 km2) of the analyzed wetland area exhibits a decrease in NDVI while 3 % (7 km2) exhibits no change in NDVI values for the years from 1985 to 2010 (Table 4). In bofedales, there are over twice as many cells with increasing NDVI compared to cells with decreasing NDVI.

Table 4 Area (km2) and percent of total area (in parentheses) of increasing, decreasing, and no change in NDVI within bofedales and otherwise in the landscape
Full size table

How do these percentages compare with NDVI change in vegetation outside of bofedales? We find that 83 % of the analyzed non-wetland cells in the watershed have increasing NDVI values while 14 % have a decreasing trend and 3 % display no change in NDVI (Table 4). There is an overwhelming percentage of vegetation in the watershed that exhibits an increase in NDVI values.

The question remains whether the decreasing vegetation trends in wetlands (29 vs. 14 % in other areas) constitutes a significant clustering of decreasing cells. It should be noted that the percentage of decreasing cells quantified in bofedales is conservative given the methods used in the regression analysis and the resulting significance filter. Furthermore, there may be areas that were once wetland but have undergone rapid and extreme changes in land cover so that they were not identified as wetlands in this study. An example of this is a portion of the large patch of decreasing cells in the south of the watershed (see Fig. 7b). Much of the area was classified as wetland on earlier maps (e.g., the 1:1000,000 National Map 1972) but due to intensive use for agriculture and cultivated pastures in more recent years it no longer meets the mean NDVI threshold to be classified as a wetland for this study.

Finally, we would like to comment on the spatial configuration of increasing and decreasing NDVI cells in the bofedales in the watershed. There appears to be no consistent pattern throughout the watershed. That is, there are bofedales that contain a majority of cells with decreasing NDVI values such as the wetland in the northwest and several of the wetland systems in the south, while others contain a majority of cells with increasing NDVI values, and still others, such as the large wetland just south of the town of Nuñoa that contain a good number of both. These spatial relationships and the processes behind them warrant further attention.

Discussion

We have used the NDVI to map wetland areas and review dry season vegetation conditions in the Nuñoa watershed between 1985 and 2010. The spatial resolution and spectral properties of Landsat 5-TM were more than adequate for the purposes and scale of this study. We were able to obtain a good number of images because relatively cloud free scenes exist for the region during the dry season (June–August). This timing is favorable since the contrast between bofedales and background vegetation is most distinct at this time. Additionally, wetlands are easily visible due to the absence of dense forest cover in the watershed (except for small areas of Polylepis forest comprising <0.2 % of the watershed). Therefore, the methods employed here in classifying wetlands would be suitable in other high elevation watersheds and regions with a semi-arid climate.

Our objective was to analyze NDVI trends in the watershed from 1985 to 2010 at a resolution of 30-m since this information is useful to local herders as well as regional, national, and international agencies with an interest in natural resources management and conservation in the humid puna. The results from the regression analysis show that there is a slight increase in NDVI with time for a majority of the watershed (81 % of all analyzed cells and in 68 % of the analyzed wetland area). A decreasing trend in NDVI occurs in 16 % of the watershed area and in 29 % of the area classified as wetland. Thus, there is a concentration of decreasing vegetation in wetland areas despite an overall greening of the watershed.

Van Leeuwen et al. (2013) conducted an NDVI analysis for a similar time period (1982–2011) and for a wider geographic region along the Andes. Vegetation trends were calculated using advanced very high resolution radiometer (AVHRR) composite derived (8 km) annual mean NDVI. They found increasing vegetation trends along the Andes of southern Peru, including the region presented here. These results suggest that the greening of the Nuñoa watershed may be part of a more regional phenomena occurring in parts of the Andes.

Many studies have related NDVI to temperature and precipitation at global and regional scales (e.g., Tucker et al. 1991; Liu et al. 2013; Piao et al. 2011; Sun et al. 2013). Productivity may be limited by either variable or more sensitive to one over the other. One limitation of the cell by cell regression method used here is that it is difficult to establish the strength of the relationship between the predictors and NDVI. To do so, Mazzarino (2014) first calculated the spatially averaged NDVI for the watershed for each of the 20 dry-season scenes (n = 20). To understand the relationship between watershed averaged dry season NDVI and precipitation during the preceding wet season, a linear regression was performed [NDVImean = α + (β0 × Precip) + ε]. Roughly half (R 2 = 0.56, p < 0.05) of the variability in mean dry-season NDVI for the watershed can be attributed to precipitation during the preceding wet season (if one assumes a stationary relationship). A time series analysis of the precipitation index used in this study however reveals no trend. Therefore, while a good percentage of the variability in NDVI is attributed to precipitation during the previous wet season, NDVI trends observed in the watershed are not attributed to a linear trend in precipitation.

Mean annual temperature in the Andean region has increased approximately 0.7 °C since the mid-1900s (Vuille et al. 2008). It might be hypothesized that the increase in temperature is playing a role in the greening of the watershed. We have accounted for interannual variability in climate by including regional temperature and precipitation indices in the regression model. The temperature index in this study was chosen based on the reasoning that vegetation might be impacted by low temperatures and associated frost just prior to scene acquisition. With our methods we cannot address the issue of how vegetation might change with climate change.

Of particular interest are the drivers of decreasing vegetation trends seen in nearly one-third of wetland area, especially since these trends are in contrast to the greening of the watershed. Again, we cannot rule out the attribution of climate factors given our methods but we hypothesize that hydrology and land use are both potential predictors of NDVI change in wetlands. We were not able to include reasonable indices for these predictors in the regression model. In a global analysis of vegetation trends in semi-arid regions, Fensholt et al. (2012) use an “evidence from absence” approach to infer secondary indicators of vegetation change. Based on our knowledge of the watershed we do the same here.

Wetland systems are sensitive to changes in hydrology and many of the bofedales in the watershed are fed by springs in addition to precipitation (ALT-PNUD 2001). It is noted that conversations with herders in the upper reaches of the watershed have elicited statements regarding decreasing volumes of water in springs and streams with significant changes circa 2000 (Mazzarino 2011, unpublished interview data). These observations are supported by studies in the region that show both strong glacier retreat since the mid-1980s (Buffen et al. 2009; Salzmann et al. 2013) and an increase in precipitation rather than snow events at higher elevations for a majority of the year (Bradley et al. 2009). Therefore, it is plausible that the decreasing NDVI trends seen in the many bofedales in the north of the watershed are related to an increasing water deficit in the dry season due to decreased groundwater contributions as a result of changing hydrological pathways. This hypothesis warrants further research.

Lastly, there are many changes occurring in the watershed with respect to land use and management. For example, there has been a trend toward increased cattle production throughout the region (INEI 2012). Unlike alpaca, cattle are quite damaging to soil and vegetation. Due to altitudinal constraints, the greatest increases in cattle numbers have been in the lower half of the watershed (PIEP 2007; Inquilla 2013). The large patch of decreasing NDVI in the south of the watershed (Districts of Orurillo and Asillo), in what was once a wetland connected to a lake system, is likely directly due to increased cattle production and agricultural activities there. However, many cells adjacent to the area display increasing trends in NDVI. The juxtaposition of the opposing trends may in fact be related in that the same districts also report an expansion of cultivated and irrigated pasture for forage (oats, alfalfa, barley) (Inquilla 2013; DIA Puno 2012).

As a result of climate disturbances and the heterogeneity of the environment, herders in the Andes have developed several land management practices that ensure increased access to critical resources. These practices include the use of a multiple resource base supported by a rotation system that includes daily and seasonal movement of herds (Thomas 1979). In recent years, however these practices have been changing. A study in 2000 by the Peruvian government stated that 44 out of 48 bofedales analyzed in the Department of Puno were now used year-round rather than seasonally (ALT-PNUD 2001). A more detailed analysis of land use and management strategies throughout the watershed combined with the results of the present NDVI analysis would be useful in eliciting relationships between these variables.

Conclusion

We demonstrated that high altitude wetland systems in the Andean region can be successfully classified and mapped using the NDVI at 30-m resolution (κ = 0.79).

A vegetation trend analysis was performed for the years 1985–2010 using a multiple regression model applied to each cell in the Nuñoa watershed. The area of study represents a region that historically and currently produces the greatest number of alpaca in Peru. The livelihood of the predominately Quechua speaking pastoralists is tightly connected to the vigor of vegetation, particularly that found in wetland systems (bofedales). The time period for the study correlates to significant changes in land use and production systems in Peru and to observable and accelerating changes in climate. Because bofedales are a natural resource heavily used by herders throughout the watershed, it is difficult with the methods employed here, to untangle the anthropogenic effects of land use and herding strategies from exogenous climate factors on vegetation. However, we have demonstrated that there are concentrations of decreasing vegetation density, as measured by the NDVI, in wetland systems in the Nuñoa watershed. This is occurring within the broader context of a greening of the puna grassland. Processes from the biophysical and social environments are surely interacting to create the pattern in changing NDVI discussed here. To what degree each process is contributing to the trends is unclear. However, the information garnered from the present study is useful in quantifying the area of and vegetation trends in a critical natural resource.