Large-scale heterogeneity of Amazonian phenology revealed from 26-year long AVHRR/NDVI time-series

Depiction of phenological cycles in tropical forests is critical for an understanding of seasonal patterns in carbon and water fluxes as well as the responses of vegetation to climate variations. However, the detection of clear spatially explicit phenological patterns across Amazonia has proven difficult using data from the Moderate Resolution Imaging Spectroradiometer (MODIS). In this work, we propose an alternative approach based on a 26-year time-series of the normalized difference vegetation index (NDVI) from the Advanced Very High Resolution Radiometer (AVHRR) to identify regions with homogeneous phenological cycles in Amazonia. Specifically, we aim to use a pattern recognition technique, based on temporal signal processing concepts, to map Amazonian phenoregions and to compare the identified patterns with field-derived information. Our automated method recognized 26 phenoregions with unique intra-annual seasonality. This result highlights the fact that known vegetation types in Amazonia are not only structurally different but also phenologically distinct. Flushing of new leaves observed in the field is, in most cases, associated to a continuous increase in NDVI. The peak in leaf production is normally observed from the beginning to the middle of the wet season in 66% of the field sites analyzed. The phenoregion map presented in this work gives a new perspective on the dynamics of Amazonian canopies. It is clear that the phenology across Amazonia is more variable than previously detected using remote sensing data. An understanding of the implications of this spatial heterogeneity on the seasonality of Amazonian forest processes is a crucial step towards accurately quantifying the role of tropical forests within global biogeochemical cycles.

There is an error in figure 7; figure 7(a) is the same as figure 6(a). The corrected figure and caption are given below.
Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Forest phenology is a critical process regulating fluxes of carbon and water between the biosphere and the atmosphere.
Biophysical and structural changes in the forest canopy driven by phenological cycles are closely tied to vegetation-climate feedbacks (Richardson et al 2013, Wu et al 2012, Gonsamo et al 2012. Therefore, depiction of the spatial configuration of phenological cycles in multiple biomes is essential for improvement of the representation of ecophysiological processes in Earth system models (Bradley et al 2011) and for evaluation of the impacts of climatic variation on these processes (Richardson et al 2013).
Undisturbed Amazonian vegetation is characterized by a heterogeneous mosaic of forest types (Silman 2007, Gond et al 2011 originating from complex interactions between three main factors: (1) climate, which generates gradients in the timing of the onset and ending of the dry season (Liebmann andMarengo 2001, Sombroek 2001), the length of the dry season (Marengo et al 2001) and the seasonality of solar radiation (Schafer et al 2002); (2) environmental factors, which present a large heterogeneity in soil physical and chemical characteristics (Quesada et al 2009), drainage and topography; and (3) biological factors, which are characterized by changes in species composition across the basin (Sombroek 2001, Sternberg 2001, ter Steege et al 2006, by structural differences in canopy sizes (Barbier et al 2009), by gradients of net primary productivity (Aragão et al 2009) and by canopy litter production (Chave et al 2010).
Previous studies using the normalized difference vegetation index (NDVI) derived from Advanced Very High Resolution Radiometer (AVHRR) data have shown little variation in the phenology of Amazonian forests (Batista et al 1997, Maignan et al 2008. Distinct changes in the amplitude of the NDVI, however, were observed for eastern (high amplitude) and north-western (low amplitude) Amazonia (Asner et al 2000). Negative anomalies in the NDVI during the dry season of eastern Amazonia were associated with decreases in canopy greenness related to dry season water stress. Moreover, the reduction in the amplitude of the NDVI from El Niño to La Niña years suggested that Amazon forest phenology is responsive to rainfall variations. An important breakthrough using MODIS data was the identification of seasonal changes in forest greening during the dry season as a response to increased solar radiation for large areas of primary forests in Amazonia (Huete et al 2006, Myneni et al 2007, Xiao et al 2006, Samanta et al 2012. However, generalization of the dry season green up over Amazonian forests must be interpreted cautiously. Bradley et al (2011), using similar data, showed that 1.27 × 10 6 km 2 of the forested area analyzed was in phase with radiation, while 1.58 × 10 6 km 2 was in phase with precipitation. Most of the forested area had a weak seasonality (Bradley et al 2011). Moreover, inter-annual variability in changes in the VI could not be explained by climatic variation in Amazonia (Brando et al 2010). This variability indicates that other environmental controls may play an important role in determining the phenology of Amazonian forests. Our knowledge about Amazonian phenology is still limited to basin-wide generalizations, mainly because of unsuccessful attempts to identify spatially variable coherent regional patterns of phenology. This is in part due to limitations linked to the length of the time-series available for previous analysis (<10 yr) and field information confined to a single area in eastern Amazonia (Tapajós National Forest;Doughty and Goulden 2008, Malhado et al 2009, Brando et al 2010 that may not be representative of other areas of Amazonia, and therefore limiting our capacity to interpret the patterns observed from satellite data. Based on the principle that phenological variation of tropical trees is shaped through adaptations to biotic and abiotic factors and, hence, phenological traits are an evolutionary reflection of the influence of these external factors (van Schaik et al 1993, Borchert 1994, Borchert et al 2002, we expect that the phenological cycles within Amazonian vegetation types will vary greatly, reproducing the environmental diversity that can emerge from the spatial arrangement of potential drivers of phenology identified above for this biome. Despite the difficulty in determining the ultimate causes of forest phenological cycles, analysis of pixel-based long-term seasonality of VI is likely to provide an indication of spatially explicit homogenous phenological responses of Amazonian vegetation to different evolutionary pressures. The longest satellite-derived VI time-series currently available is the NDVI from AVHRR (1981-2011) (Zhu et al 2013). A recent comparison between the NDVI and the enhanced vegetation index (EVI) from MODIS revealed that bidirectional reflectance effects are more exacerbated in the EVI than the NDVI, potentially influencing the interpretation of phenological cycles in forest ecosystems (Galvão et al 2011). To identify distinct ecological responses of Amazonian forests to localized biotic and abiotic controls, in this study we propose to analyze 26-year long time-series of NDVI data from AVHRR using a novel pixel-level pattern recognition technique based on temporal signal processing concepts. Specifically, we aim to (1) characterize the spatial variability of forest phenology in Amazonia by identifying regions where the forest phenological cycle is homogeneous (phenoregions); (2) quantify the relationships between the seasonality of the NDVI signal of each distinct phenoregion and climate by analyzing the covariance between NDVI, rainfall and solar radiation; and finally (3) investigate how the NDVI changes in relation to field-derived data on leaf flushing and litterfall based on a new compilation of published information for Amazonia.

Datasets
We used two long-term NDVI time-series  with 8 km spatial resolution from the Advanced Very High Resolution Radiometer (AVHRR) available from the Global Inventory Monitoring and Modeling Studies (GIMMS) group. The first corresponds to biweekly data and the second to monthly data. These two datasets have been constantly improved to minimize noise resulting from residual atmospheric effects, orbital drift, inter-sensor variations, and stratospheric aerosol effects . The biweekly data were used for the pattern recognition analysis due to their higher temporal resolution. The monthly data were used to perform the correlation analysis between monthly rainfall and incoming radiation and NDVI datasets, as well as to investigate the seasonality of field phenology, NDVI and rainfall data.
Monthly precipitation and incoming radiation time-series covering the period from 2000 to 2008 were created for each phenoregion based on the Global Land Data Assimilation System (GLDAS) precipitation grids (0.25 • spatial resolution) (Rodell et al 2004). Additionally, we compiled published information on field observations of peak of new leaf and litterfall production for 12 sites in Amazonia to investigate how the NDVIs for selected phenoregions respond to the seasonality of these two phenological processes: (1) leaf flushing and (2) abscission (table S1 available at stacks.iop. org/ERL/8/024011/mmedia).
To reduce the effect of land use and land cover change on the phenological signal, we masked all areas corresponding to agriculture, water body, bare soil and urban areas based on the GLOBCOVER product developed by the European Space Agency (ESA) (Arino et al 2007).

Data processing
Our rationale for detecting forest phenology was based on three assumptions: (1) NDVI time-series can provide cyclical spectra that are directly related to the distinct phenological patterns of each phenoregion; (2) the periodicity of the cycles can be mathematically separated, allowing the identification of phenologically homogenous regions; and (3) the multi-temporal nature of the NDVI time-series is analogous to the characteristics of hyperspectral data, permitting the use of similar image processing techniques.
We applied a signal processing approach to the NDVI time-series aiming to automatically identify NDVI standard curves, which represent phenological forest patterns, in order to map these regions. This signal processing approach consisted of 'noise-adjusted' principal component analysis (minimum noise fraction), automated extraction of end members (pixel purity index) and an automated classification method based on spectral similarity (spectral angle mapper).
The minimum noise fraction (MNF) is a linear transform that provides an optimal ordering of images in terms of image quality and increases the signal to noise ratio, generating principal component images that are unaffected by noise (Green et al 1988).
After MNF transformation, the pixel purity index (PPI) was then applied to the principal component images. The PPI is an algorithm based on convex geometry concepts used to find the most 'spectrally pure' or extreme pixels (end members), which examines multidimensional envisioned data in n-dimensions. Each point (pixel) within this data space can be examined as a linear combination of an unknown number of pure components (Boardman 1995). The PPI is computed by repeatedly projecting n-dimensional scatter plots onto a random unit vector. The extreme pixels in each projection are recorded and the total number of times each pixel is marked as extreme is noted. A threshold of 5000 pixels was used to define how many pixels are marked as extreme at the end of the projected vector.
The selected end members (phenological signatures) were then used as input for the spectral angle mapper (SAM) analysis. This pattern recognition technique is based on the measurement of spectral similarity between two or more spectra obtained considering each spectrum as a vector in n-dimensional space (Kruse et al 1993). Once similar spectra are identified they are aggregated into a single category. Lastly, the resulting map was postprocessed and pixels with less than nine neighbors and phenoregions with less than 100 pixels were excluded from the final map.
We computed for each phenoregion the mean monthly NDVI. These time-series were smoothed using the Savitzky-Golay filter (Chen et al 2004) and the coefficient of variation was calculated as a homogeneity index. For each NDVI time-series extracted from each phenoregion, phenophases corresponding to the start (SS), end (ES) and peak (PS) of the growing season were calculated according to the methodology proposed by Jönsson and Eklundh (2002).
Finally, correlation analyses were performed to quantify the relationships between NDVI and rainfall and radiation for each phenoregion.

Results and discussion
The pattern recognition technique enabled the characterization of 26 phenoregions with distinct long-term NDVI temporal cycles (figure 1). These phenoregions tended, in general, to reflect the spatial distribution of main vegetation types described for Amazonia (table 1).
Previous studies used approaches based on automated cluster analysis or training techniques using multisensor data (Saatchi et   forest types have structural and climate differences that have been previously described (Sombroek 2001, Gond et al 2011 but not yet clearly separated using remote sensing data. Another important distinction, probably influenced by climatic gradients, was observed among dense forest, open forest and dry forest in the north-south direction.

Reconciling NDVI variation with field-based phenology
For all phenoregions analyzed, the peak of new leaf production reported from field observations was followed by a steady increase in NDVI or occurred during the maximum NDVI (figure 2). The peak of litterfall production corresponded to time periods in which NDVI had a descending trend (figure 3).
In contrast to previous studies using MODIS data that reported widespread leaf flushing during the dry season in Amazonia (Huete et al 2006, Myneni et al 2007, Xiao et al 2006, Samanta et al 2012, field-based information from our literature compilation revealed that the peak of new leaf production occurs during the wet season in 66% of the sites. The NDVI from the AVHRR follows this same dynamics with values increasing from the beginning or middle of the rainy season to the end.
Specifically, field-based reports showed that the pluvial forest in high lands (Phenoregion 18) and forest transition between dense, open forest and savanna (Phenoregion 24) presented just one annual peak of new leaf flushing in the dry season (Boubli 2005, Pinto andSetz 2004). The closed evergreen forest (Phenoregion 9) presented two peaks of flushing in the wet season (Peres 1994). The savanna forest (Phenoregion 6), mosaic of evergreen lowland forest, swamp forest, and savanna (Phenoregion 25) and high forest with regular canopy (Phenoregion 16) presented just one peak of flushing in the wet season (Norconk and Conklin-Brittain 2004, Basset et al 2001, Bonal et al 2000. The NDVI of each phenoregion followed these events in both the dry and the wet season. The coefficient of variation of the NDVI values did not exceed 40% in any of the regions analyzed. The automated PPI method captured NDVI temporal signatures consistently with the main phenological features observed in the field.

Timing of phenological phases
We identified a strong rainfall and radiation gradient, with maximum rainfall occurring in parallel with minimum radiation. Both variables followed a southeast-northwest direction. The phenophases clearly followed this climatic gradient (figure S2 available at stacks.iop.org/ERL/8/024011/ mmedia). In the northwest of the basin the onset of the growing season occurred mainly between May and July ( figure 5(a)), during the middle to the end of the dry season when solar radiation is at its maximum ( figure 4(c)). This region is characterized by a non-existent or short dry season. On the other hand, in the central-north Amazon the onset of the growing season occurred from March through May, during the middle to the end of the wet season, following the increased availability of solar radiation. In the extreme south of the Amazon basin, Phenoregions 17 and 26, characterized by submontane forest and savanna, respectively, had the onset of the growing season between October and December, two months before the maximum rainfall (beginning of the wet season). In the southwest-northeast direction the peak of the growing season ( figure 5(b)) occurred following the middle-to-end of wet season gradient, from April (southwest) to September (northeast). The end of the growing season (figure 5(c)) also followed the end of dry season gradient in the southwest-northeast direction (September-February) with a lag of about three months.
Correlation analysis revealed that without lags only savanna forests presented positive correlation with rainfall. However, 76.9% of all phenoregions had a significant positive correlation between NDVI and rainfall with a three month time lag (figure 6). The correlation between NDVI and incoming radiation was positive in a large part of the basin without lags (figure 7). These results indicate that large areas of Amazonian forests strategically flush new leaves from the middle to the end of the wet season when incoming radiation is increasing and soil water availability is still high (Juárez et al 2007). A delay effect of the NDVI in response to rainfall is expected and is likely to explain the occurrence of negative correlation coefficients when lags are not considered. An exception was savanna forest and submontante forest with palms and lianas that presented positive correlation between NDVI and rainfall without time lags.
The NDVI amplitude between maximum and minimum values reached 16.67% in the mosaic of evergreen lowland forest, swamp forest, and savanna (Phenoregion 22) and 39.07% in the mixed high and open forest with a tropical savanna climate (Phenoregion 25). In the alluvial dense evergreen forest domain in the low lands (Phenoregion 3), where Tapajos National Forest and Caxiuana National Forest are located, the NDVI amplitude reached 21.72% (table S1 available at stacks.iop.org/ERL/8/024011/mmedia). These results were similar to previous studies that used EVI (Samanta et al 2012, Myneni et al 2007, Huete et al 2006. Our results showed a basin-wide phenological gradient delay between the southeast and northwest portions of four months for the start of the growing season and six months for the peak of growing season, when the forest reached the highest NDVI values at the middle-to-end of the wet season. This reveals that the phenological responses to climatic variation are heterogeneous across Amazonia. Based on our NDVI and field-based analysis we argue that the flush of new leaves does not occur simultaneously during the dry season across all Amazonian forest as previously stated (Samanta et al 2012, Myneni et al 2007, Saleska et al 2007, Huete et al 2006, Asner et al 2004. By combining the long-term NDVI time-series (26 years) with ground-based studies produced in the last 26 years, we highlighted that leaf flushing in the majority of Amazonia, as approximated by the NDVI, starts after the soil is fully recharged with water and incoming radiation starts to increase after the peak of the wet season. This tendency is characterized in our data by a peak in the NDVI by the end of the rainy season. We demonstrated that the onset and end of the growing seasons in Amazonia cannot be considered to occur simultaneously throughout the basin from July to September. We showed a southeast-northwest direction gradient of rainy and dry seasons (figures 4 and 5). It was also possible to note a difference of five months in the peaks of the dry and wet seasons in the southeast-northwest direction and four months in the start of the growing season (September to December) in the west-east direction of the Amazon basin, which is in agreement with Marengo et al (2001).
It is interesting to note that in a large fraction of southern Amazonia, the growing season occurs almost exclusively during the wet season for a period of three months. In contrast, Figure 7. The Pearson correlation coefficient between NDVI and incoming radiation (a) without lag, (b) with one month lag, (c) with two month lag and (d) with three month lag. The black color corresponds to unclassified areas in the phenoregion map (see figure 1). in the northern, and wetter, fraction of the basin the growing season occurs during the transition between wet and dry seasons for a period of six months. The shorter time of the growing season in the south could be associated with a 'sprint' of greenness as an evolutionary strategy to reduce losses through herbivory during this susceptible phase of foliar development (Patino et al 2012, Moles andWestoby 2000). This 'sprint' of greenness could be either associated to the flushing of new leaves (Sinimbu et al 2012) or to the increased rate of leaf expansion (Brando et al 2010).
Previous studies of Amazonia greening tend to assume a dry season between July and October, concluding that Amazonia greens up during the dry season (Samanta et al 2012, Huete et al 2006, Xiao et al 2006. Therefore, due to the variability observed here in the timing of phenological phases and climatic gradient, more caution should be taken when linking climatic conditions with canopy phenology in Amazonia.

Conclusions
This work provided a critical understanding of the phenological mosaics in Amazonia and how they respond to climatic variability. The fundamental aspect of our methodology was the implementation of a pixel by pixel pattern recognition technique that could successfully identify the differences in the natural cycles of forest phenology. Our phenoregion map clearly distinguishes regions that were not detected in previous studies. This is probably because spatially pixels are spectrally similar but temporally they produce distinct responses to biotic and abiotic factors.
Differently from previous studies, using MODIS data, which generated contradictory results, creating an unresolvable debate (Saleska et al 2007, Asner and Alencar 2010, Samanta et al 2012, Atkinson et al 2011, our results added a new and independent perspective on the phenology of Amazonian forests. The coherent phenological patterns found in our analysis are a result of the combination of the data characteristics and the methodological approach used. In spite of inherent limitations of optical remote sensing, in terms of saturation of NDVI signal over tropical canopies (Aragão et al 2005) and interference of clouds (Asner and Alencar 2010), it appears that, in particular, the length of the time-series allowed the detection of cyclical behavior in the AVHRR/NDVI data. This indicates that saturation of NDVI was not a problem for the analysis. Moreover, the identification of peak NDVI during the dry season in some areas and the wet season in other areas indicates that cloud cover is not a major factor influencing the seasonality of this dataset. Because our methodology explicitly distinguishes curves based on their shape and not their amplitude, the influence of NDVI saturation and potential cloud or aerosol interference is minimized.
Forest phenology is a complex process involving several environmental factors which enables similar vegetation types to develop distinct phenological cycles or, alternatively, different vegetation types to develop similar phenological cycles. For this reason we opted to compare the phenoregions with the main vegetation types on a macroscale.
Our results show that Amazon forest phenology is more heterogeneous than previously anticipated. To better understand the drivers of this heterogeneity it is critical to perform a formal evaluation of the contribution of other variables across the basin, including soil, geomorphology, vegetation and climatic information.