Trends and ENSO/AAO Driven Variability in NDVI Derived Productivity and Phenology alongside the Andes Mountains

Increasing water use and droughts, along with climate variability and land use change, have seriously altered vegetation growth patterns and ecosystem response in several regions alongside the Andes Mountains. Thirty years of the new generation biweekly normalized difference vegetation index (NDVI3g) time series data show significant land cover specific trends and variability in annual productivity and land surface phenological response. Productivity is represented by the growing season mean NDVI values (July to June). Arid and semi-arid and sub humid vegetation types (Atacama desert, Chaco and Patagonia) across Argentina, northern Chile, northwest Uruguay and southeast Bolivia show negative trends in productivity, while some temperate forest and agricultural areas in Chile and sub humid and humid areas in Brazil, Bolivia and Peru show positive trends in productivity. The start (SOS) and length (LOS) of the growing season results show large variability and regional hot spots where later SOS often coincides with reduced productivity. A longer growing season is generally found for some locations in the south of Chile (sub-antarctic forest) and Argentina (Patagonia steppe), while central Argentina (Pampa-mixed grasslands and agriculture) has a shorter LOS. Some of the areas have significant shifts in SOS and LOS of one to several months. The seasonal Multivariate ENSO Indicator (MEI) and the Antarctic Oscillation (AAO) index have a significant impact on vegetation productivity and phenology in southeastern and northeastern Argentina (Patagonia and Pampa), central and southern Chile (mixed shrubland, temperate and sub-antarctic forest), and Paraguay (Chaco).


Introduction
The Millennium Ecosystem Assessment [1] illustrates how South American ecosystem services are governed by high biodiversity, productive land cover, a variety of land use changes, and by desertification that is taking place in some areas.One attribute used to monitor and assess ecosystem services is annual above ground primary production.Vegetation productivity has been estimated by using NDVI time series data [2][3][4].Additionally, long term NDVI times series have been used to analyze trends [5][6][7][8][9] and responses of vegetation to environmental variables such as climate [6,[10][11][12][13] and land use [14][15][16][17][18].
The Intergovernmental Panel on Climate Change (IPCC) indicated that phenology, the timing of recurring biological events in response to the environment, is a relatively simple measure or proxy for climate change and variability [19].Land surface phenology, is defined as "the seasonal pattern of variation in vegetated land surfaces observed from remote sensing" [20], is key to many earth surface processes and impacts carbon and hydrological processes as well as societies and economies [21].Therefore, information about land surface phenological responses is needed to enhance earth systems science and model representation and understanding.Spatially distributed land surface phenology has been derived using NDVI time series data from the Advanced Very High Resolution Radiometer (AVHRR) and Moderate resolution Imaging Spectroradiometer (MODIS) to examine trends and the timing of landscape changes [22][23][24][25][26] and monitor and assess ecosystem responses [27][28][29][30][31], and possibly adapt to asynchronous ecosystem response to environmental variability and change [32].South America's land cover is extremely variable across the continent [33] and is strongly impacted by interannual variability in climate.The topography of the Andes Mountains in combination with climatic cycles such as ENSO results in complex seasonal variability and spatial distribution of precipitation and temperature fields across the continent [34].Interannual climate variability is causing increased temperatures and droughts in some and increased rainfall and decreased solar irradiance in other areas [34].Some of this climate variability is captured by atmospheric indices such as the Multivariate ENSO Index (MEI) and Antarctic Oscillation index (AAO).One of the objectives of this research is to explore how MEI and AAO are affecting vegetation response.A better understanding of vegetation-climate interaction in South America will address efforts associated with adaptation to climate variability and warmer temperatures, which has consequences for ecosystem functioning, and land and water use across the continent in a variety of ways.For example, agricultural, mining and energy related demands for water are putting pressure on optimizing and managing these often limited resources.Furthermore, changes in the timing of biological events, such as the timing of green-up, flowering and bird migrations can impact water demand, rates of carbon sequestration, and overall ecosystem functioning.
The main objective of this research is to characterize and examine the interannual variability, trends and changes of land surface productivity and phenology in response to climate variability and change using the longest and newest generation of continuous normalized difference vegetation index (NDVI3g) time series data set currently available [35,36], in combination with climate data.We focus this research on the southern cone of South America alongside the Andes Mountains, including Argentina, Bolivia, Chile, Paraguay, Uruguay, and parts of Peru and Brazil.In the context of the objective we are asking the following research questions: (1) What are the productivity and phenological characteristics of the diverse land cover types in South America?(2) Are there any interannual trends in the phenology?(3) How does land surface phenology respond to the environment (e.g., climate, land use)? ( 4) What and where are the impacts of climate on local and regional scale land surface responses?
With respect to questions (3) and ( 4), we will mostly focus on examining the functional relationships between interannual climate indices (ENSO and AAO), and vegetation productivity and phonological variables.

Data and Methods
The 30-year record of AVHRR GIMMS NDVI data were investigated for part of South-America alongside the Andes, between the latitudes of −9°S and −60°S and the longitudes of −80°W and −50°W, to analyze changes in inter-and intra-annual trends as a function of several climate indices to infer possible causes for observed inter-annual changes in greenness.The following set of methods were used to answer the questions raised in the introduction: (1) trend analysis of annual mean NDVI; (2) trend analysis of phenological variables-start and length of season climatic variables; (3) correlation analysis between annual mean NDVI and seasonal climate variables (MEI and AAO); and (4) correlation analysis between phenological variables (SOS and LOS) and seasonal climate variables (MEI and AAO).

Study Area-South America
The study area is defined by the outlined extent in Figure 1.The focus is mostly on arid and semi-humid areas of the southern cone of South America, a region that has not been a focus of climate-vegetation response research using vegetation time series data spanning 30 years.Figure 1 provides an overview of the aggregated land cover types based on work published by Eva et al. [33].Eva et al. [33] represents multiple agricultural and forest types, such as the endemic temperate deciduous forests in Chile, better than most other comparable land cover products.The Andes Mountains provide for complex topography and climate patterns along a large latitudinal gradient.From a biogeographic viewpoint, the study area contains four bioclimatic zones (Table 1).They are formed by the presence of the Andean mountain in the west, a precipitation gradient with an arid diagonal zone in the middle of the continent (from the Atacama desert, to La Pampa), and by a north-south temperature gradient associated with latitude [37][38][39].The four bioclimatic zones contained within the study area have been transformed by human activities related with forest extraction for energy, land habilitation for agriculture and forest plantation, and urban expansion during the last two centuries.A summary of the areal extent of land cover types is provided in Supplementary Table S1.The climate of this particular area of the continent is impacted mostly by the ENSO and AAO patterns.This region presents a wide variety of climates as a consequence of atmospheric circulation patterns (represented mainly by the Intertropical Convergence Zone and easterly winds in the tropics and the predominance of Westerlies, prevailing winds at southern mid-latitudes), the influence of oceanic currents that transport heat (for instance the Peru, Cape Horn, and Brazil currents) and the presence of major geographic features, like the Andes, that contribute to create a mosaic of south America's climate types.
The influence of the Andes becomes particularly important for the explanation of a semi-arid diagonal strip of land that covers the central north part of Chile and the Patagonian region of Argentina (Figure 1).A rain shadow effect is created by this mountain chain.While between 15°S and 27°S latitudes precipitation is dominated by convective activity in the Amazonia, the Westerlies bring moisture and frontal systems to the southern latitudes.The presence of the Andes precludes rainfall to occur on the western side of the Andes at low latitudes and on the eastern side at high latitudes [40].
The relevance of agriculture and forestry in this continent is quite significant.In this region, a small subsistence agriculture coexists with a more sophisticated export oriented industry, representing a very dynamic sector for the economy (on average, this is a sector that represents approximately 4% of the gross domestic product [41]) and a relevant source of employment.

GIMMS-NOAA-AVHRR NDVI Time Series Data
Time-series of biweekly (15 days) composited 8-km resolution Normalized Difference Vegetation Index (NDVI) values from July of 1981 to December of 2011 provided by the Global Inventory Modeling and Mapping Studies (GIMMS) project were acquired by the Advanced Very High Resolution Radiometer (AVHRR) sensors 7, 9, 11, 14, 16 and 17 on board the National Oceanic and Atmospheric Administration (NOAA) satellite platforms [35,36].The NDVI is the ratio of the difference of the near-infrared (NIR) and red reflectance (ρ) values (ρ NIR − ρ red ) divided by the sum of the red and NIR reflectance values (ρ NIR + ρ red ).The NDVI time series data were calibrated and corrected for view geometry, volcanic aerosols, and other miscellaneous issues, unrelated to vegetation response, were minimized [35,36].These NDVI time series data are extensively used for trend and land surface response analysis [2,6,31,42,43] and currently spans 30.5 years, 24 images per year, resulting in 732 gridded images with observations of vegetation dynamics.The continuity of multi-sensor NDVI data sets remains a challenge due to differences among sensors with respect to spectral response functions, calibration, sensor drift, and atmospheric correction causing small changes in NDVI dynamic range among sensors [35,[44][45][46][47][48][49].The most recent comparison between AVHRR and MODIS based NDVI values showed good correlation for phenologically distinct areas and larger differences for tropical forest areas [42].We will explore the continuity of the 30-year NDVI time series using the Atacama Desert as an invariable target.
NDVI data are also very helpful to examine the capacity of carbon assimilation and water use of different vegetation types during their growing season [50,51].NDVI values integrated during the growing season are frequently used as a proxy measurement for seasonal vegetation productivity [52].In this research the annual mean NDVI is used as an NDVI-based productivity metric.The annual mean NDVI value for each year (t) is computed by averaging all 24 observations for each pixel for each year, referred to as NDVI mean,t .Since vegetation growth tends to occur in spring and summer in the Southern hemisphere it was necessary to redefine the period of integration as the one between July and June, spanning two calendar years.This method has been applied to other similar research done in this region [18].This metric is expected to represent different and unique spatio-temporal vegetation responses to global and regional drivers like climate and land use.

Land Surface Phenology-NDVI Time Series Data
The NDVI time series data provide seasonal and interannual information about ecosystem responses to the climate.These phenological trajectories represent ecosystem changes and variability in response to environmental dynamics including changes in land use.Furthermore, the seasonal and interannual trends provide a qualitative indication of the continuity of the NDVI data across AVHRR sensors.
TIMESAT time-series analysis software was used to extract metrics of vegetation dynamics (phenological metrics) for each year (1981−2011).An adaptive Savitzky-Golay smoothing filter [53] was used because it maintains distinctive vegetation time-series trajectories and minimizes various atmospheric effects [29,54].For the purpose of this study, two timing metrics were used: (1) the start of growing season, the time when the NDVI value has increased by 20% [54] of the distance between the pre-season minimum and the seasonal maximum; (2) the length of growing season, the difference between the start and end dates of the growing season.The end dates of each growing season were determined based on the same 20% threshold.The timing based pheno-metrics are used to characterize seasonal photosynthetic activity and patterns among different vegetation types as they respond to different environmental conditions along latitudinal and topographical gradients.

Climate Indices-Monthly Multivariate ENSO Index (MEI) and Antarctic Oscillation (AAO) Data
ENSO has a very strong impact on the interannual variability of precipitation and temperature patterns in many regions of South America [34].The Andes Mountains also play a pivotal role in these climate patterns.Global scale ENSO-related rainfall anomalies show that the El Niño periods are usually correlated with anomalously wet conditions in central Chile and the southeastern portion of the continent and with below normal rainfall and warmer than normal conditions in the northern part of South America.La Niña events are typically associated with opposite rainfall anomalies in both regions [34].The Multivariate ENSO Index (MEI) represents multiple variables (e.g., sea level pressure, sea surface temperature, Niño regions), related to El Niño and La Niña phenomena [55].A positive MEI is associated with El Niño and a negative MEI with La Niña periods.Decadal and interdecadal variability has been associated with the Pacific Decadal Oscillation (PDO) anomalies of precipitation and temperature over South America, which have a similar spatial structure to the ENSO anomalies, but less strong [34].Finally, positive phases of the Antarctic Oscillation (AAO) typically cause warmer surface air temperature to the south of 40°S.Precipitation anomalies related to the AAO are largest at 40°S, and substantial in southern Chile and along the subtropical east coast of South America [34].
In this research we examine the relationship between MEI and annual mean productivity and phenological metrics.We also examine the impact of the AAO on the same variables.Since the impact of the PDO is relatively weak, we will not investigate this further within the scope of this research.Monthly MEI [55,56] and AAO [57,58] data were used to calculate mean tri-monthly MEI s and AAO s values to represent the growing seasons (s): Winter (JAS), Spring (OND), Summer (JFM), and Fall (AMJ).

Ecosystem Characterization: Variability and Trends in Productivity and Phenology During 30 years
The thirty year mean/median and interannual variability of ecosystem productivity and phenology are represented by the mean and coefficient of variation (COV) of NDVI mean,t and LOS t data, and the median and standard deviation (STD) of SOS t data.The thirty year average annual NDVI was computed based on the average of the thirty yearly mean values of NDVI mean,t .The COV(y) is computed for the thirty years of NDVI mean,t and LOS t data using the following equation: where STD is the standard deviation of the thirty observations per pixel.Interannual trends in the thirty year data record were examined for the NDVI mean,t , SOS t and LOS t using a first order linear regression analysis.Each of the three dependent variables are analyzed as a function of time represented by thirty consecutive years (t = 1 to 30).The following equation is an example of the functional relationship establishing the trend in the interannual NDVI thirty year time series data: where a is the slope, b is the intercept, and e t the error term of the least squares solution.Only significant trends for p ≤ 0.05 are reported.

Ecosystem Response to MEI and AAO
One of the main causes of the variability in ecosystem response is climate variability.We examine how the thirty years of NDVI mean,t , SOS t and LOS t data are related to the MEI s and AAO s indices by using a first order linear regression analysis.Antecedent MEI s and AAO s conditions, preceding the annual growing season, were included in the analysis for the summer and fall seasons.Coincident MEI s and AAO s conditions were included for the winter, spring, summer and fall seasons.Each of the three dependent variables is analyzed as a function of these seasonal and interannual varying climate indices.The following equation is an example of the functional relationship between the annual NDVI and the MEI climate indicator: where a is the slope and b is the intercept and e t the error term of the least squares solution.Significant trends for NDVI mean,t are reported for p ≤ 0.05.The levels of significance for SOS t and LOS t results are set to p ≤ 0.1.

Land Surface Response and Phenological Characteristics and Interannual Variability
The characteristic seasonal responses and interannual variability (standard deviation based on 30 years per composite period) for a range of vegetation land cover samples are shown in Figure 2. The samples are a subset of the samples identified in Figure 1.The NDVI time series data show how each of South America's land cover types have diverse seasonal patterns.Vegetation cover types, such as tropical forest, grasslands, agriculture, savannah, shrublands and deserts, all have land cover type specific seasonal NDVI curves with differences in seasonal dynamic amplitudes and ranges as well as differences in the onset, peak and end of the growing seasons.The timing of the maximum NDVI values varies considerably throughout the year: tropical forest peaks mid-winter, while sub-humid, temperate and arid vegetation types mostly peak during the October thru February time frame.The timing of the minimum NDVI values varies considerably throughout the year as well: tropical forest have much lower NDVI values during the rainy season, while sub-humid, temperate and arid vegetation types mostly have minima values during the fall and winter (April-September).It should be noted that the NDVI time series for the desert areas are very flat for these multi-sensor time series NDVI data.
The magnitude and variability of the standard deviation during these seasonal NDVI patterns are harder to interpret as they represent the interannual variability among all the 30 years of NDVI observations within each of the composite periods.This interannual variability is often caused by year to year changes in the timing and magnitude of the seasonal growth patterns, thus providing us with information about ecosystem response to the environment.This is especially relevant to, for example, a water-limited and possibly temperature-limited shrubland, which will not start greening up until there is soil moisture available and temperatures allow vegetation growth.Hence the variability in the timing of precipitation is often coincident with variability in the timing of NDVI increase.On the other hand, the timing of the onset of green-up and senescence in a temperature-limited Montane forest (>1,000 m) shows the highest variability during periods of water and temperature limitations, while showing very low variability during the winter.However, the interpretation of the interannual variability based on the standard deviation is often more complicated in tropical or temperate vegetation types (e.g., Closed evergreen tropical forest, closed deciduous forest, and temperate mixed evergreen broadleaf forest) where cloud cover and/or smoke aerosols also will cause some of this variability.Since cloud cover and aerosols will reduce the NDVI independent of the vegetation cover beneath it, some of the larger standard deviations that we observe are most likely due to cloud/aerosol impacts on the NDVI data.For the closed evergreen tropical forest seasonal NDVI profile, the large difference between summer and winter NDVI and standard deviation values is arguably mostly due to clouds.The seasonal and interannual variability in the NDVI for a few sample points are further illustrated in Figure 3.The closed deciduous forest point represents a decrease in NDVI of more than 5% while the barren point represents a decrease of less than 2%, significant at p ≤ 0.05.
The 30 years of biweekly NDVI data show temporal consistency for the desert areas (Figure 3), despite the fact that the data comes from multiple sensors.The barren soil 30-year trajectory seems to make a slight drop in 2004.This drop coincides with the time when AVHRR-NOAA-16 transitioned to AVHRR-NOAA 17.When we do a simple linear regression on the annual mean of the barren area, there is a small decline in the NDVI of about 0.01 actual NDVI units over 30 years (p ≤ 0.05).The desert area seems to have a slightly higher NDVI response in the last few years, reportedly due to some sporadic rainfall events.El Niño was strong in 1997/1998, followed by La Niña years in 1999 and 2000.The precipitation limited shrubland and grassland samples show an above normal NDVI response coinciding with increased MEI values and a lower NDVI due to low MEI values.The 30 years of NDVI observations for the closed deciduous forest sample show an overall decrease over time, but also illustrate the volatile responses within short time periods.These are most likely due to cloud cover.Hence, it is best to use the upper envelope, or higher NDVI values, of these NDVI trajectories to minimize these adverse effects.
The 30-year annual mean NDVI for the selected vegetation types shows a large range of values (Figure 4).The lowest land NDVI values are for barren soil and desert land cover types, while the highest NDVI values are for temperate mixed forest.One can assume that the NDVI values for closed evergreen tropical forest are suppressed by frequent cloud cover.One sample site classified as a water body showed higher NDVI values, likely caused by a mixture of water and vegetation mixed within the pixel.For each of the sample sites a summary with the interannual means, COV and trends are included in Supplementary Table S2.
Figure 4. Thirty-year mean annual NDVI and ±STD for representative vegetation cover types samples (see Figure 1).
The spatial distribution of the 30-year annual mean NDVI and the associated COV are shown in Figure 5.The low mean NDVI values or low productivity areas are diagonally situated across the Andes Mountains.Strong increasing productivity gradients are observed from the desert coast across the Andes mountains into tropical forests in the northern latitudes (North of −32°S latitude), while we observe the opposite pattern in productivity in the more temperate southern latitudes (~−32°S latitude), temperate forest on the coast and shrubland across the Andes mountains to the east.Tropical and temperate forested areas show the regions with the highest productivity, followed by agriculturally intensive areas.Many of the patterns seen in the mean NDVI in Figure 5 follow along well established precipitation and temperature gradients classified by the recently updated Köppen-Geiger system map [59].Most of the higher interannual variability represented by higher COV of the annual NDVI values are associated with areas on and along the Andes Mountains, complemented by some of the shrubland regions in the central-northern parts of Chile and southwestern parts of Argentina.These shrubland regions are often drought deciduous/tolerant and will have strong leaf-out and growth of herbaceous plants in response to rainfall.The median of 30 years of SOS values and the STD of the 30-year mean SOS are shown in Figure 6.The SOS values respond variably to temperature, moisture, and radiation limitations [2,27] or a combination thereof.Most of the SOS values seem what we would expect.There is a strong latitudinal and elevational gradient in the SOS along Chile and the Andes Mountains.It should be noted that the SOS values for the tropical regions are likely less reliable due to the strong interference by aerosols and cloud cover that impacts the NDVI data [60].The high STD values seem to confirm this.We also have some cross calendar year or cloud/snow issues in the south of Chile if the high STD values are a good indication.The main reason that the median of the SOS is displayed in Figure 6 instead of the mean SOS is that the mean of the SOS is not a useful value due to the fact that, for certain pixels, the start of the growing season is for some years in November or December (Day of year ~300-365) and for other years in January or February (Day of year ~1-60).Hence the mean SOS of some vegetation types will not result in the correct mean of growing seasons that cross or skip years.The patterns and the highly variable amplitude of the STD for the 30 years of SOS values demonstrate both areas where the mean SOS value is based on end and start of calendar years as well as where the SOS is variable.Further work needs to be done to accurately visualize the interannual variation in SOS data.The spatial distribution of the timing of the SOS for the study area (Figure 6) is consistent with our expert knowledge of the region.The SOS, in concert with the STD, will give an indication where the SOS is likely to cross calendar years (e.g., for SOS values in December and January and possibly for other SOS values representing months further removed from these two months), and where the SOS data are very variable (e.g., tropical forests-cloud frequency; and deserts-rainfall events).The mean of 30 years of LOS values and the COV are shown in Figure 7.The LOS values, just as the SOS, also respond variably to temperature, moisture, and radiation limitations [2,27] or a combination thereof.There is a strong latitudinal and elevation gradient in the LOS along Chile and the Andes Mountains.It should be noted that the LOS values for the tropical regions are likely less reliable due to the strong interference by aerosols and cloud cover that impacts the NDVI data [60].Extraction of the LOS from the seasonal NDVI curves seems to underestimate the LOS, particularly in the tropical forested areas.We also have some cross calendar year or cloud/snow issues in the south of Chile if the high COV values are a good indication.The patterns and the highly variable amplitude of the COV values for the 30 years of LOS values are demonstrating both areas where the mean LOS value is variable and where clouds could cause some less reliable retrieval from the seasonal curves.

Interannual 30-Year Trend of Annual Productivity and Phenology
Significant (p ≤ 0.05) positive and negative trends in the annual mean NDVI data over a 30-year time span are shown in Figure 8.The trends are in percent NDVI units per year.Green colors relate to positive trends, while yellow to dark red correspond to negative trends.The strongest negative trends (>5%) are observed in southeast Bolivia and northwest Paraguay and Central Argentina.Intermediate negative trends (−2% to −5%) are observed for many parts of the shrublands and some of the agricultural areas in Argentina.The Atacama Desert also had an unexpected, but significant, weak negative trend.Since this desert is one of the driest deserts of the world, we suspect that this negative trend is not due to less vegetation but rather due to multi-sensor spectral response and possibly atmospheric effects.Positive trends in the annual NDVI are seen in central and southern Chile dominated by a mix of agriculture (some tree plantations) and temperate forests.Various areas with significant positive trends can be observed within Brazil, northeast Bolivia and parts of Peru, Paraguay and Uruguay.The coastal desert in Peru has a long term NDVI mean of close to zero and the positive trend is likely an anomaly related to the NDVI time series data.Some of the negative trends in Bolivia, Paraguay and Argentina have been shown with less detail in other studies that were using shorter time series data [6,18,26,[61][62][63][64].One study used annual GIMMS NDVI data (1982-1999) and found the positive trends in south Chile but did not find as many negative trends [18].A detailed summary of the positive and negative trends in the annual mean NDVI during the last 30 years are provided per land cover/use type and country in Supplementary Table S3.
We also applied a Kendall and Spearman test (results are not included in this manuscript) and found the same regions and patterns to have significant trends in productivity and phenology as the linear regression results (Figures 8 and 9).It should be noted that most regions have SOS values that differ 30-60 days or even more among the 30 years.Some areas have a tendency to cross calendar years.We did not account for this interannual cross year variability in SOS, which will likely cause an underestimation of the number of pixels with significant trends.This is a drawback of this methodology that requires more attention in future research approaches.LOS is not impacted by cross calendar year issues as it is calculated separately from SOS within the TIMESAT software.

The Impact of MEI and AAO on Productivity and Phenology
The 30-year NDVI time series of a closed shrubland are shown together with MEI and AAO to demonstrate the nature of these data (Figure 10). Figure 10 also shows some general indication of lower seasonal NDVI responses during lower MEI values or La Niña phases and higher during positive MEI values or El Niño phases.The AAO and MEI seem to have a tendency to have opposite signs; while the MEI is positive the AAO index is often negative.
Figure 10.The AAO data tends to be more variable than the MEI data.We therefor plotted a 3-month moving average of the AAO data with the monthly MEI and the NDVI data.The latitude and longitude for this closed shrubland sample are indicated in the legend.

Vegetation Response (Mean NDVI) as a Function of MEI and AAO
After running the correlations between antecedent and coincident seasons of MEI and AAO, the winter MEI and summer AAO index (p ≤ 0.05), coincident with the July to June growing season, were determined to have the strongest correlations and coherent spatial patterns.Although not presented in this manuscript, the only other significant season that showed a good relationship with productivity was the antecedent summer (the summer of the year before) AAO index.Figure 11 shows where the annual mean NDVI values were significantly (p ≤ 0.05) correlated with the MEI winter and the AAO summer index.The MEI winter is mostly positively correlated with the annual mean NDVI for mixed shrubland and agricultural areas in central Chile, open shrubland areas in southeastern Argentina and the pampas with grasslands and agriculture in Uruguay.Most of these regions seem to respond strongly to the increased precipitation amounts associated with higher MEI winter values or El Niño phases.We observe both positive and negative significant correlations (p ≤ 0.05) between the annual mean NDVI and AAO summer index.The positive correlations are found in the Patagonia region of southern Chile (Figure 11).Since the AAO tends to be associated with significantly higher temperatures for the southern cone [34], this area has higher productivity for positive AAO summer phases.This part of Chile is more temperature limited [2], and thus will increase its productivity with higher temperatures.Many of the areas with negative correlations are associated with some of the same areas where positive correlation with MEI winter was measured.Some of these areas seem to become more water limited with increases in temperature or decreases precipitation related to the AAO summer phases in these regions [34] and will reduce productivity.The open shrublands of southern Argentina do not seem to be impacted by warmer temperatures during higher AAO index values, which could indicate that this vegetation is more resilient to temperature shifts.

Phenology Response (SOS and LOS) as a Function of MEI and AAO
The areas with significant phenological responses (p ≤ 0.1; annual SOS and LOS) to the winter MEI and summer AAO index are shown in Figure 12.The MEI winter is both positively and negatively correlated with the SOS in the study area, indicating a later and earlier SOS respectively.Some of the stronger spatial patterns related to a later SOS (positive correlation with MEI winter ) can be observed in the pampas of Argentina and Uruguay.Some other notable spatial patterns are related to an earlier SOS (negative correlation with MEI winter ), which can be observed in some of the shrublands of Argentina, central Chile and Bolivia.For the aforementioned regions, the earlier SOS often corresponds to longer LOS and vice versa.The patterns in the SOS and LOS phenological timing variables show that both the MEI winter values or El Niño phases and the AAO summer have an impact.We observe both positive and negative significant correlations (p ≤ 0.1) between the annual SOS and LOS and the AAO summer index.These SOS-AAO summer related patterns are often opposite those of the SOS-MEI winter related patterns.Later SOS patterns in northern Argentina and Paraguay are more unique to the SOS-AAO summer .It is clear that the timing of the start and length of the growing season show many complex relationships with MEI winter and AAO summer values.These are likely related to the timing of precipitation and temperature fluctuations that are impacted by climate circulation patterns.The next phase of this research could further explore these relationships using spatially and temporally distributed precipitation and temperature data, where available, and possibly cloud cover frequency and daytime length.However these analyses are not within the scope of this manuscript.

Land Use Change Impacts on Trends in Productivity and Phenology
The 1981-2011 3G NDVI time series data provides many new insights about the productivity and phenological characteristics and dynamics for a diverse range of land cover types across South America.
New are the trends of increasing productivity for some of the areas in Chile with temperate forests mixed in with agriculture.The reasons for this increase in productivity are not always clear.However, some evidence in this research showed that warmer temperatures in the sub-antarctic and temperate vegetation cover types in Chile are associated with the AAO, possibly causing more productivity.Land use change due to an increase in tree plantations (Pinus radiata and Eucaliptus globulus) in this region in Chile [65] could also contribute to increased productivity.Higher rainfall events during El Niño phases impacted central Chile but did not result in a significant productivity increase during the last 30 years.This same AAO might cause a more drought related decrease in productivity in some of the open and closed shrublands of southern Argentina.
Large regions with larger > 2% decreasing NDVI trends are found east of the Andes.Some of the negative NDVI trends in the dry Chaco regions of central and northern Argentina, northwestern Paraguay and southeastern Bolivia are mostly related to changes in land use [66][67][68].
Many of these natural areas have been degraded over time due to cattle ranching and wood extraction, while other areas have seen an increase in soybean area/production and other agricultural crops over the last decade [69,70], replacing some of the open deciduous forest areas.Soybean, pasture and natural vegetation mosaics in Brazil [71] just east of Argentina and Paraguay, on the other hand show an increase in the NDVI trend.If we examine some of the smaller change locations in Brazil, north of Bolivia, a pattern of tropical deforestation coupled with expansion of agricultural becomes evident, leading to reduced productivity while replacing more productive tropical and subtropical forests.Some of the new areas where decreasing interannual trends in productivity are reported raise more questions than answers.What causes the slight decrease in NDVI in the Atacama Desert?Is this decrease related to land cover which is predominantly bare soil, atmospheric aerosols, or multi-sensor continuity?
The main areas with a later start of growing season are in the Patagonia steppe the Chaco region and Agricultural regions in Brazil.Most of these areas also show significant productivity trends.Land management and agricultural practices likely play a role in these areas; however it is not clear what the impact of climate is on these trends since this was not explored in this research.Shorter lengths of the growing seasons over the last 30 years are clearly associated with agricultural land use in central Argentina.Many of the significant shifts in SOS and LOS range from 30 to about 90 days over this 30-year period.More research will be needed to explore the causes of these phenological shifts further.

Impacts of MEI and AAO Climate Patterns on Productivity and Phenology
MEI and AAO explain interannual variability in productivity and phenology for some regions in the study area (Figures 11 and 12).MEI related precipitation variability in several areas in Chile [72] and South America [34] are clearly associated with the interannual variability in annual mean NDVI, predominantly shrubland areas.These locations also coincide with areas where the COV of the annual mean NDVI (Figure 5) tends to be higher.The productivity of rain fed agriculture in the Pampas are also sensitive to the variability in MEI [73], but do not show a high COV value for the annual mean NDVI values.The AAO reduces productivity for areas that are associated with higher temperatures and limited precipitation.Diagonally across the Andes, higher temperatures in central Chile, the Pampas, and shrublands in southern Patagonia will increase evapotranspiration and reduce productivity for these more water limited systems.This is exacerbated by reduced rainfall in the Pampas and Central Chile during positive AAO phases that often seem to coincide with La Niña phases.
Many of the significant patterns associated with the SOS and LOS as a function of the MEI and the AAO index coincide with large scale patterns of climate variability.However, more somewhat random patterns with some significance are associated with the correlation between phenology and MEI or AAO data are shown as well (Figure 12).It is not always clear what all these more random patterns mean.A more detailed or local scale analysis with finer resolution data seem to be required to explore this further.

Conclusions
The third generation of the 8 km resolution Normalized Difference Vegetation Index (NDVI3g) time series data (1981-2011) provides unprecedented opportunities to examine changes and variability in Earth system response and processes over a 30 year time period.This research provided new baselines and characteristics of South America's of productivity and phenological trends and the impact of climate indicators.Both positive and negative trends were observed for Chaco and Patagonia.Land use and climate variability and change were related to the observed productivity and phenology trends and pattern.Importantly, this research showed how the Multivariate El Niño Southern Oscillation (ENSO) Index (MEI) and Antarctic Oscillation (AAO) index explained the variability in annual productivity and phenology for some distinct areas in South America that are impacted by precipitation and temperature variability associated with these climate circulation patterns.
We highlighted significant land cover trends and variability in annual productivity (represented by the mean annual NDVI) and land surface phenological (represented by the start and length of the annual growing seasons) response.Negative trends in productivity are observed in arid and semi-arid Some of the seasonal and annual NDVI baseline data and phenological variables are variably impacted by clouds and aerosols.Some caution should be observed to not over interpret the results for the more tropical areas.Similarly for the desert areas where slight decreases in NDVI trends are likely more related to the multi-sensor continuity characteristics of the NDVI3g data set than actual changes in the landscape.
This research highlighted firsthand the 30-year long-term trends and spatial variability in productivity and phenology and their response to the MEI and AAO indicators, providing critical insights into South America's ecosystem dynamics and functions.This study provides ample motivation and a clear context to further explore and attribute causes of land surface change and variability, possibly using observation at finer spatial scales and a focus on certain time periods.A logical future step in this effort is to explore the impact of precipitation and temperature on the trends in NDVI and phenology time series data, and further explore the coupled land use and climate impacts on these trends and patterns of change and interannual variability.

Figure 1 .
Figure 1.Land cover types based on data from Eva et al. [33].The black circles are the sample areas representative of the major land cover types that are used to show the phenology and some examples of the seasonal and interannual NDVI time series data within the indicated study area extent.

Figure 2 .
Figure 2. Seasonal NDVI response and variability for representative vegetation cover types-sample areas, 9 pixel averages, based on the mean NDVI value composite period over a 30-year period.The standard deviations (vertical bars) are indicative of the 30-year interannual variability.

Figure 3 .
Figure 3. Interannual variability in phenological trajectories for some of the sample points.The closed deciduous forest point represents a decrease in NDVI of more than 5% while the barren point represents a decrease of less than 2%, significant at p ≤ 0.05.

Figure 5 .
Figure 5. Spatial distributions of the annual mean NDVI and its coefficient of variation (COV) representing the interannual variability.

Figure. 6 .
Figure. 6. Median of the SOS based on thirty years of NDVI data.The STD is provided for the annual SOS data.It should be noted that the STD is very variable because of cross calendar year SOS values or high interannual variability.

Figure 7 .
Figure 7. Thirty-year mean LOS and coefficient of variation based on the phenology of thirty years of NDVI data.

Figure 8 .
Figure 8. Significant positive and negative trends (p ≤ 0.05) in annual mean NDVI based on 30 years of continuous NDVI time series data (1982-2011) are shown in the right side figure.Land cover is provided as a reference.

Figure 9 .
Figure 9. Thirty-year significant trends (p ≤ 0.05) in the SOS and LOS indicating regions where SOS is earlier or later and where the LOS is longer or shorter.Many of these shifts span 30 to 90 days over 30 years, but are often very variable from year to year.

Figure 12 .
Figure 12.The start and length of the growing seasons (SOS and LOS) as a function of MEI and AAO displays great variability for the period of 1982-2011(p ≤ 0.1).
and sub humid vegetation types across Patagonia-Argentina, the Atacama Desert in northern Chile, and the Chaco region in northwest Uruguay and southeast Bolivia.Positive trends in productivity are observed for parts of the temperate forest and agricultural areas in Chile and sub humid and humid areas in Brazil, Bolivia and Peru.The start of the growing season (SOS) results highlights the large phenological variability that exists from 1982 to 2011.Years with later SOS values often coincide with reduced productivity.Growing seasons that have lengthened over the 30 year study period are mostly found for areas in the south of Chile (sub-Antarctic forest) and Argentina (Patagonia steppe).The central Pampa in Argentina with mixed grasslands and agriculture tends to have shorter growing seasons.The seasonal Multivariate ENSO Indicator (MEI) values for the winter months and the Antarctic Oscillation (AAO) index values for the summer months have a significant impact on vegetation productivity and phenology of the following regions: mixed shrubland in central Chile and temperate and sub-Antarctic forest in southern Chile, Patagonia in southeastern Argentina and the Pampa in northeastern Argentina, and the Chaco region in Paraguay.