Satellite observation of tropical forest seasonality: spatial patterns of carbon exchange in Amazonia

Determining the seasonality of terrestrial carbon exchange with the atmosphere remains a challenge in tropical forests because of the heterogeneity of ecosystem and climate. The magnitude and spatial variability of this flux are unknown, particularly in Amazonia where empirical upscaling approaches from spatially sparse in situ measurements and simulations from process-based models have been challenged in recent scientific literature. Here, we use satellite proxy observations of canopy structure, skin temperature, water content, and optical properties over a period of 10 years (2000–2009) to constrain and quantify the spatial pattern and seasonality of carbon exchange of Amazonian forests. We identify nine regions through an optimized cluster approach with distinct leaf phenology synchronized with either water or light availability and corresponding seasonal cycles of gross primary production (GPP), covering more than 600 million ha of remaining old growth forests of Amazonia. We find South and Southwestern regions show strong seasonality of GPP with a peak in the wet season; while from Central Western to Northeastern Amazonia cover three regions with rising GPP in the dry season. The remaining four regions have significant but weak seasonality. These patterns agree with satellite florescence observations, a better proxy for photosynthetic activity. Our results suggest that only one-third of the patterns can be explained by the spatial autocorrelation caused by intra-annual variability of climate over Amazonia. The remaining two-thirds of variations are due to biogeography of the Amazon basin driven by forest composition, structure, and nutrients. These patterns, for the first time, provide a complex picture of seasonal changes of tropical forests related to photosynthesis and influenced by water, light, and stomatal responses of trees that can improve modeling of regional carbon cycle and future prediction of impacts of climate change.


Introduction
The terrestrial gross primary production (GPP) is considered the largest CO 2 flux (123 ± 8 petagrams) and responsible for driving several ecosystem functions globally (Beer et al 2010). Estimates of magnitude and regional variations of this flux remain uncertain in humid tropical forests, particularly in Amazonia (Huntingford et al 2013, Schimel et al 2014) where limited ground data have caused gross assumptions about seasonality and heterogeneity of these forests (Clark 2004). Forests of Amazonia, although considered evergreen, appear to have seasonal cycles, following the rhythms of rainfall and radiation (Myneni et al 2007, Fu et al 2013 that significantly impact their carbon assimilation and GPP (Clark 2004, Lee et al 2013, Fischer et al 2014. Experimental studies and ecosystem modeling suggest that availability of water and light regulate carbon assimilation of the Amazon forests (Saleska et al 2003, Baker et al 2008, Kim et al 2012. Eddy covariance measurements in Central Amazonia indicate that there is an increased net ecosystem productivity and evapotranspiration during dry season that may suggest the greening of canopy or emergence of new leaves (Hutyra et al 2007, Restrepo-Coupe et al 2013). However, ground observations in Southern and Southwest Amazon support more of a distinct dry season decline in canopy photosynthesis due to water stress (Vourlitis et al 2005, Saleska et  This diversity in results point to a significant heterogeneity in vegetation interaction with climate, creating spatially variable seasonality of GPP that cannot be resolved at landscape to regional scales with limited ground and tower flux data (Restrepo-Coupe et al 2013). In general, GPP is not directly observed at the plot or tower flux sites but approximately estimated from either observations of net primary production (NPP) and autotrophic respiration (R a ) as GPP = NPP + R a from detailed plot level measurements  or derived from tower eddy covariance measurements of the net ecosystem exchange. In either case, the uncertainty in GPP estimation is high and is strongly impacted by the accuracy of ground measurements (Malhi et al 2009) and various sources of errors in tower observations (Restrepo-Coupe et al 2013). Given the lack of reliable direct in situ measurements of GPP and the potential heterogeneity of GPP at the landscape scale, making the few plot or tower based estimates sparse and inadequate (Restrepo-Coupe et al 2013, Malhi et al 2014), consequently, spatially refined estimates of GPP are often derived from remote sensing products or a combination of ecosystem models, climate, and remote sensing observations. The use of various satellite measurements directly related to vegetation biophysical parameters, such as the leaf area index, fraction of absorbed photosynthetically active radiation (fPAR), or vegetation greenness indices like the normalized difference vegetation index (NDVI) (Hashimoto et al 2012) are key variables to examine the variations and seasonality in GPP related to vegetation productivity and photosynthetic activities.
Studies using optical satellite data and vegetation indices have shown that the signals related to vegetation photosynthesis increase during the light-rich dry season in Amazonia (Myneni et al 2007, Samanta et al 2012, suggesting the seasonality of photosynthesis is driven by the light-controlled leaf phenology through the root water dynamics or under no water stress (Kim et al 2012). However, the observed seasonal 'green-up ' (Huete et al 2006, Myneni et al 2007 has been challenged in recent studies using estimates from satellite fluorescence sensors (Lee et al 2013), and also modeling effort showing sun-sensor geometry effects (Morton et al 2014). The multi-year analyses of these data have also shown that optical satellite observations in the tropics suffer from atmospheric contamination (Samanta et al 2010, Maeda et al 2014. These conflicting results point to the need for a comprehensive interpretation of satellite and ground observations for capturing regional variations of vegetation phenology and productivity. Current efforts on examining and mapping the spatial patterns of forest seasonality have been focused on either the non-tropical regions (White et al 2005), or the signal processing of multi-temporal observations from a single satellitederived variable, such as NDVI (Silva et al 2013, Hilker et al 2014. But the use of green index alone is often not sufficient to capture the GPP changes in view of the modeling approaches such as the light-use-efficiency (LUE) model (Monteith 1972, Running et al 2004 or the water-use-efficiency model (Beer et al 2007(Beer et al , 2010. Here, we quantify the spatial patterns of seasonality of the Amazonian forests by combining 3 longterm satellite observations over a decade (from 2000 to 2009). First, we use land surface temperature (LST) (Wan 2008) aboard Moderate Resolution Imaging Spectroradiometer (MODIS) as an approximate measure of the canopy skin temperature and a key control on photosynthetic activity and GPP (Sims et al 2008). LST is strongly correlated with PAR and vapor pressure deficit (VPD) (Sims et al 2008), both of which are essential for quantifying GPP spatial variations and seasonality. Second, we use microwave radar backscatter aboard QuikSCAT (QSCAT) as a proxy for canopy water content and VPD-a variable controlling the photosynthetic activity, which can be monitored all time irrespective of the presence of clouds, aerosols, or seasonality of incoming radiation (Frolking et al 2011, Saatchi et al 2013. And third, we included the Nadir BRDF (bidirectional reflectance distribution function) adjusted reflectance (NBAR) of near-infrared (NIR) band (Schaaf et al 2002(Schaaf et al , 2011 from MODIS satellite, as the proxy for illuminated canopy structure without the influence of sensor geometry (Knyazikhin et al 2013; SI discussion), a variable related to the FPAR seasonality and LUE used in GPP estimation (Monteith 1972). After careful quality checks, we expect that the average seasonality observed from the three different satellite proxy measurements over 10 years can reasonably delineate the spatial patterns of forest seasonality, and separate the forests of Amazonia into dominant phenological regions. The detected patterns are important for understanding the carbon exchange and the gross primary productivity (GPP) of Amazonia and may inspire future regional studies and in situ observations.

Methods
We used multiple satellite-derived products in this study, including greenness data-NDVI and enhanced vegetation index, QSCAT radar measurements, LST and NBAR products from MODIS, suninduced chlorophyll fluorescence (SIF) data from Greenhouse gases Observing SATellite (GOSAT), rainfall estimation from the Tropical Rainfall Measuring Mission (TRMM) product, and downward surface shortwave radiation from the Clouds and the Earth's Radiant Energy System (CERES) (SI materials and methods). The latest MODIS land cover product was used to define the research area and only pixels identified as evergreen broadleaf forests were selected to perform further analysis of seasonality in Amazonia. To evaluate the impact of observed satellite signals on the GPP estimation, we used the upscaled GPP from Max-Planck-Institute for Biogeochemistry (MPI-BGC) (Jung et al 2011) as a reference data set. The upscaled GPP has been developed by using the global network of tower flux measurements to reduce the uncertainty predicting regional variations along climate gradients, even in regions with limited tower measurements as in Amazonia. However, we consider the MPI-BGC GPP data may still have an uncertainty associated with individual grid cell, but using the data to examine seasonal variations at large regional scales could reduce the uncertainty by spatial filtering or averaging.
We made use of TIMESAT algorithm (Jönsson and Eklundh 2002, Eklundh and Jonsson 2011 to fill data gaps, smooth the signals and find seasonality parameters. The monthly climatology averaged from the decade of 2000s (2000-2009) was used as the input to obtain seasonality parameters. The K-means clustering algorithm was performed to create cluster maps from LST, NIR, and QSCAT. Sensitivity analyses using Random Forest (RF) machine learning algorithm (Breiman 2001) and multiple linear regression were performed to identify the important remote sensing variables for GPP prediction. For the regional GPP prediction, we used the RF model to develop a relationship between MPI-BGC GPP and LST, NIR and QSCAT data (figure 2, figure S7, SI materials and methods).

Results
We found nine phenological regions (figure 1(a)) from the average monthly climatology of LST, NIR and QSCAT using K-means clustering (Seber 2004, Jain 2010)-an unsupervised measure to find similar features from multiple bands. The number of phenoregions was internally validated (figure S1) using both stability-based and variance-based statistical measures (Caliński and Harabasz 1974, Davies and Bouldin 1979, Tibshirani and Walther 2005 to ensure the clusters are well-separated (SI materials and methods). The clusters represent regions with (a) similar seasonal variations or phenology related to the amplitude and phase of observations (i.e. intra-annual variability), (c) visualization of the cluster result in PC3 and PC4 domain (SI materials and methods). Colored region is classified as the Evergreen Broadleaf Forests (EBF) according to the latest MODIS land cover map. The boxes on the left panel were selected for regional studies (e.g. figure 2). Numbers in red color of panel (a) are selected regions for further discussion in the main text. Red cross signs denote the field sites available to validate our results (table S2; SI discussion). The direction and length of the vector in panels (b) and (c) indicate how each variable contributes to the two PC components. Each variable (LST, QSCAT and NIR) has 12 such vectors representing monthly climatology from January to December. and (b) similar landscape structure and biogeography related to annual mean values of observations (Steege et al 2013). We used principal component (PC) analysis (Hastie et al 2009) to visualize these two effects and found that the 1st and 2nd PC domain, containing 65% information of all observations, explain mainly variations in annual mean value ( figure 1(b)). This was also confirmed using the clustering map of annualmean-only data (figure S2) that explains 67% of the regional variations (table S1) in our clustered map ( figure 1(a)). The seasonal effect (amplitude and phase) starts dominating the 3rd and 4th PC domain and explaining about 25% of all observations ( figure 1(c)). Although the latter contribution is minor, it provides distinct spatial features (figure S2) highlighting the contribution of intra-annual variability in each pheno-region. The importance of this contribution can be readily demonstrated by statistically optimizing the number of clusters (figure S2).
Spatially, LST, NIR, and QSCAT all have different phases with rainfall or radiation seasonality (figure S3), indicating that climate alone cannot fully explain the geographical differences in vegetation seasonality. The seasonal variations of remote sensing observations appear to be significant even in the extremely wet regions of Amazonia ( figure S4). The inclusion of these observations makes the detected seasonality patterns more spatially correlated and distinct with more than 70% of the cluster results predictable from nearby pixels (table S1). As our method never assumes any spatial autocorrelation, these results indicate the importance of spatial variations in amplitude, phase, and annual mean of remote sensing measurements in defining patterns of seasonality ( figure S5).
Out of the nine pheno-regions, we delineated 7 large contiguous regions for detailed study of seasonality (boxes in figure 1(a)), and among them we selected three boxes with distinct seasonality for comparison: (1) the Southwestern Amazon (region 2) with ecologically dry seasons (EDS; monthly rainfall ⩽ 100 mm) extending 3-5 months, (2) the Central Amazon (region 3) with only 1 to 2 EDS months, and (3) the Eastern Guiana shield (region 6) with similar EDS as in region 2 but with a strong seasonal swing of radiation (figure S6). All three regions show decreasing canopy water content from QSCAT and increasing leaf temperature from LST during the dry season (figures 2(a)-(c)), but have obvious differences in NIR, with peaking in the middle of wet season in region 2 ( figure 2(a)), about the end of dry season in region 3 ( figure 2(b)), and to the middle of dry season in region 6 (figure 2(c)). The range of LST and the amplitude of QSCAT also vary from region to region. In Southwestern Amazonia, QSCAT and LST are strongly out of phase, with peak loss of canopy water content (drop in QSCAT) matching the maximum leaf temperature (increase in LST) during the middle to the end of dry season. In contrast, NIR peaks in wet season but start increasing from the middle of dry season as a response to deeper penetration of NIR signal and strong sensitivity to exposed understory vegetation. From Southwest to Northeast Amazonia, NIR's peak gradually shifts to the middle of dry season as radiation plays increasingly a dominant role in vegetation seasonality and promoting photosynthetic activity in dry season when light is abundant and water is still sufficient (low QSCAT amplitude). The region in Central Amazonia falls in between the two extremes with complex seasonality controlled by a combined water and light influence represented by the lagged behavior of QSCAT, LST with NIR.
These variations in observed vegetation phenology have an aggregate effect on the seasonality of GPP. Using upscaled GPP data as the reference (Beer et al 2010, Jung et al 2011), we found regional GPP seasonality being in phase with NIR during the dry season, remaining high in the wet season in agreement with QSCAT and/or LST, and strongly following the patterns of SIF data from the GOSAT (Frankenberg et al 2011, Lee et al 2013)-an independent measure of GPP directly linking to the photosynthetic activity (figures 2(d)-(f)). GPP peaks during the wet season in Southwest and dry season in Northeast regions of Amazonia respectively responding directly to water and radiation abundance.
These distinct patterns further confirm the heterogeneity of seasonality of forests of Amazonia related to regional variations of water and light constraints, directly addressing the question of whether Amazonian forests are water-or light-limited. We found two regions in Southern Amazonia (regions 1 and 2) with GPP peaks during the wet season (figure 2(d), figure S7); two regions in Central and Western Amazon (regions 3 and 4) showing rising GPP during the dry season (figure 2(e), figure S7); the region in Eastern Guiana shield (region 6) has GPP rising and peaking during the dry season (figure 2(f)); and the remaining two regions (regions 5 and 7) are less seasonal but showing two weak seasons within an annual cycle caused by lagged intensity of water and radiation (figures S7 and S8). The seasonal patterns of GPP are present in North-South and East-West gradient in all vegetation characteristics measured by the remote sensing data (figure S9). The overall agreement of upscaled GPP data with satellite observations to capture the large scale regional patterns of seasonality suggest that the satellite observations may be used to predict regional patterns of the carbon exchange seasonality in Amazonia. These regional patterns are extensive in size and internally heterogeneous in forest cover and landscape features, making the direct ground verification of GPP or vegetation seasonality difficult, unless through systematic and widespread measurements. However, similar patterns have been observed with spatially limited but comprehensive ground measurements of forest productivity, leaf fall and flushing, and seasonal carbon fluxes

Discussion and conclusions
The spatial heterogeneity in GPP seasonality cannot be explained by only a single vegetation variable. The canopy temperature and structure from LST has a nonlinear relationship with GPP ( figure 3(a)), showing a steady increase until reaching a threshold and then followed by a gradual decline. This behavior agrees with observations showing an optimum range of temperature for forest photosynthesis (Berry andBjorkman 1980, Wood et al 2012). GPP shows a strong sensitivity to water deficit as measured by QSCAT backscatter, but quickly decouples from water when forest canopy accumulates enough water and has no limitations for photosynthesis ( figure 3(c)). Over some wet regions of Amazonia, GPP may even drop with increasing water content (regions 3, 4 and 6), pointing to low solar radiation or water-logged vegetation, both impairing photosynthetic activity (figure S10; SI discussion). Though regionally it appears to have tight negative correlation between canopy water and GPP, it may not truly reflect a causal relationship between them because of the small range of variations in the saturation zone of QSCAT (figure 3(c)), and the covariance with other variables such as NIR and LST. We found MODIS NIR being in phase and monotonically increasing with GPP over the Amazon basin ( figure 3(b)), but also have large variations due to the heterogeneity of vegetation composition and seasonality, and potential presence of atmospheric effects during the wet season (figures 2(d)-(f)).
From the multi-sensor data analysis, we found a strong positive relation between NIR and GPP and SIF (figure S10) except in region 5 where there is very low seasonality and all vegetation variables including the relationship between GPP and SIF remain weak (figure S10). Vegetation seasonality captured by three remote sensing data (LST, QSCAT, NIR) together explain 60% of the GPP over the entire Amazon Basin with much larger explanatory significance in more seasonal regions. The variations around the mean in the relationship between the magnitude of GPP and remote sensing data are largely due to the heterogeneity of forest composition, canopy structure, and nutrient availability impacting forest carbon cycling  figure 1(a); right panels show the seasonal cycles of NIR, MPI-BGC GPP, GOSAT sun-induced chlorophyll fluorescence (SIF), and the predicted GPP curve (SI materials and methods) in (d) region 2, (e) region 3, and (f) region 6. The colored circles in the right panels denote the fraction of invalid data in each NIR observation period. Spatial averages in all regions take only pixels from the EBF class into calculation. Precipitation values are from the TRMM product (SI materials and methods). The annual cycle (January-December) is replicated 3 times for better demonstration. and exchange. We also expect the uncertainty in the benchmark GPP due to limited tower data in tropical regions may have contributed to the lack of correlation between GPP and satellite observations of forest canopy. Our results also show that the NDVI vegetation greenness index is not suitable for tropical seasonality analysis in terms of GPP variation (SI discussion).
In analyzing passive optical data (LST and NIR) over each region, we made sure that all regional patterns are extracted from a majority of the pixels in each region that are valid after careful data quality check (SI discussion). However, the atmospheric contamination from clouds and aerosols may still affect the surface retrievals due to small cumulus clouds within a pixel that cannot be readily detected by MODIS cloud detection tools (Koren et al 2008). With the presence of residual atmospheric effect in the pixels labeled as cloud-free, the spectral changes in the optical domain can be complicated. One could expect that NIR may either be elevated due to undetected clouds or reduced due to shadows. And LST, on the other hand, could be suppressed under persistent clouds. Even if we used the MODIS climatologically averaged data over 10 years with both geometric corrections and high quality flag filtering, and working with large area averages for clustering, the use and interpretation of these optical data in the tropical regions should still be treated with extreme caution. The large scale spectral variability of tropical forest due to canopy phenology makes it difficult to investigate whether the extent of changes in optical observations are due to physiological response of leaf fall and flushing, or the residual effect of small scale cloud contamination. A good example is the poor quality of optical data in the wet season, where we see the divergence between NIR and GPP explicitly presented in region 3 (figures 2(e)) and 4 ( figure S7(b)). For seasonal changes of GPP in this period, it is necessary and essential to include the microwave data into the data analysis to compensate for the loss of optical data quality (SI discussion). The evaluation and uncertainty assessment of the optical data use in the tropical regions must be performed by using techniques that include airborne and in situ data.
A handful of flux tower measurements of GPP seasonality are the only ground-based resource existing in the Amazonian forests that are suitable for intercomparison with satellite observations. We selected 5 tropical forests sites in Amazonia and digitized the GPP curves from the research article of Restrepo-Coupe et al (2013) (figure 4). Results show that the in situ GPP curves are in phase with the NIR curves, particularly in the dry season and following the same regional patterns (figure 2). Though the magnitude of All observations were resampled to be the same resolution as MPI-BGC GPP (0.5°× 0.5°) using spatial average. The linear regressions were calculated from all pixels within the selected ranges (i.e. 2 selected ranges for LST: <298K and >300K; 1 selected range for QSCAT: <0.18; and the full range for NIR and SIF), and marked in red for the first selected range and green for the second selected range (LST only). GPP from in situ and gridded products are not comparable due to the normalized units used in in situ GPP, the MPI-BGC GPP are well in phase with the tower data, partially because they are essentially integrated in the MPI-BGC GPP benchmark data. However, the most interesting finding of the comparisons with the in situ data is the overall agreement with the regional pattern in each pheno-region that we defined purely from remote sensing data. Flux sites in Manaus, Tapajós, and Caxiuanã all show a rising GPP during the dry season with a peak at the end of dry season, consistent with their regional average represented by region 3 (figure 2(e)), while the other two tower sites in Ji Paraná and Rio Javaes-Bananal Island, belonging to the pheno-region 1, both exhibit GPP peaking during the middle of wet season and about 2 to 3 months later than the peak shown in region 3. These results, though limited in quantity, prove the validity of our pheno-region clusters and the general agreements between GPP and remote sensing observations of canopy characteristics. With more validation resources being available in the Amazonia (Galbraith et al 2014), the prediction strength of GPP seasonality is expected to be improved significantly.
The carbon exchange of tropical vegetation, whether light-controlled or water-limited, can vary not only in space, but also in time. The illustration using the average GPP response to different incoming radiation and precipitation (figure 5) conceptually shows that the flux tower sites located in different regions can respond to large variations in climatedriven seasonality. The São Gabriel da Cachoeira tower (Saleska et al 2009) in Northwestern Brazilian Amazon (monthly rainfall > 200 mm) has high GPP and less seasonal amplitude due to small changes in radiation and rainfall throughout the year (see region 5, figures S7 and S8). In contrast, the Ji Paraná tower site (Saleska et al 2009, Restrepo-Coupe et al 2013 in Southern Amazonia has enough radiation throughout the year (>200 Wm −2 ) causing GPP varying strongly along the precipitation axis. In the case of the Central Amazon tower site located near the city of Manaus (Restrepo-Coupe et al 2013), GPP has a moderate seasonal variation driven by a hybrid effect of rainfall and radiation causing the GPP increase from October to May rainy season, a short drop in May and June with the decline of rainfall, and a sharp increase from late June to late September from increased radiation.
Using more accurate estimation of chlorophyll fluorescence from GOSAT, we show that while Amazon forests as a whole retain high productivity during the year, there are distinct regional patterns of seasonality. These patterns are potentially driven by variations in forest species, structure, nutrient availability, and climate variations over Amazonia, creating complex mechanisms to regulate photosynthetic carbon assimilation through water and light availability and dynamic stomatal responses (Lee et al 2013). Contrary to recent findings, our analysis of multiple satellite observations along with ground measurements of leaf dynamics (Galbraith et al 2014) indicate that Amazonia does not maintain a consistent canopy throughout the year, and variations of carbon exchange in Amazonia are reflected by potential changes of top canopy structure regionally during dry and wet seasons. Regional differences in vegetation seasonality in the Amazon basin also demand further research and measurements to understand and quantify the vulnerability of each pheno-region to future changes of climate and seasonality. A phase shift off the stable pattern of seasonality these forests could mean increased risks of ecosystem stress and disturbance, water loss, and carbon release, which can have significant impacts at the global scale given the huge amount of carbon stored in the Amazon forests. Figure 5. Relationship between estimated average GPP and climate variation in Amazonia. The GPP value is the average estimation of GPP predicted from the LST, NIR and QSCAT, given the condition of radiation and precipitation in each 1°box. Monthly total precipitation is from the TRMM product, and monthly average downward surface shortwave radiation is from the CERES product (SI materials and methods). The selected seasonal cycles were also plotted at locations where insitu measurements can be obtained from flux towers.