Evaluating and Quantifying the Climate-Driven Interannual Variability in Global Inventory Modeling and Mapping Studies (GIMMS) Normalized Difference Vegetation Index (NDVI3g) at Global Scales

: Satellite observations of surface reflected solar radiation contain information about variability in the absorption of solar radiation by vegetation. Understanding the causes of variability is important for models that use these data to drive land surface fluxes or for benchmarking prognostic vegetation models. Here we evaluated the interannual variability in the new 30.5-year long global satellite-derived surface reflectance index data, Global Inventory Modeling and Mapping Studies normalized difference vegetation index (GIMMS NDVI3g). Pearson’s correlation and multiple linear stepwise regression analyses were applied to quantify the NDVI interannual variability driven by climate anomalies, and to evaluate the effects of potential interference (snow, aerosols and clouds) on the NDVI signal. We


Introduction
Over the past 30 years atmospheric CO2 levels have been rising as a result of fossil fuel emissions. Superimposed on the atmospheric CO2 growing trend are large interannual excursions (nearly 100% of the trend), and evidence strongly suggests that this interannual variability in the atmospheric CO2 growth rate is driven largely by the terrestrial biosphere [1,2]. Global carbon cycle models are challenged to reproduce this CO2 growth rate variability, and explanations largely involve independent responses of terrestrial net primary production (NPP), heterotrophic respiration (RH), and fire emissions to climate anomalies such as those associated with ENSO events [3][4][5][6]. Understanding the potential for the terrestrial biosphere to mitigate or perhaps aggravate CO2 accumulation in the atmosphere may be as important for predicting future climate change as projections of future fossil fuel emissions. Improved forecasts of future trends in atmospheric CO2 will be possible, through the development of models that can realistically represent the processes controlling current and past terrestrial carbon fluxes [7].
Precise modeling of NPP is key in terrestrial carbon cycle models, because NPP not only represents the net uptake of carbon from the atmosphere by vegetation but also affects other carbon fluxes by providing the substrate for RH and fuel for combustion by fire. Terrestrial NPP has been showed to play a central role in determining the local and global CO2 content of the atmosphere at temporal scales spanning hours [8] to epochs [9].
Long-term satellite observations of vegetation "greenness" provide a record of variability that greatly improves modeling NPP in terrestrial carbon cycle models. Satellite derived greenness refers to various surface reflectance indices that can be measured with satellite instruments. These indices express the amount of green vegetation absorbing solar radiation, a state that is often nearly proportional to photosynthetic productivity [10][11][12] and net ecosystem carbon exchange [13]. Global carbon cycle models often rely on satellite observations of vegetation greenness to either prescribe phenological variability [14] or to evaluate vegetation dynamics in prognostic models [15].
Polar orbiting satellites have provided global measurements of the greenness index Normalized Difference Vegetation Index (NDVI) at sub-monthly time steps and/or moderate spatial resolutions (30 m) for decades. NDVI is commonly used to derive canopy solar radiation absorptance for estimating NPP [3,16,17]. It has also been used in global climate models to prescribe vegetation influences on the hydrological and surface energy cycle [18][19][20][21]. NDVI has been shown to indicate variability in vegetation productivity in response to climate variability [22] and human management [23]. NDVI may not capture all the climate driven variability in NPP, however, and includes noise and biases that interfere with diagnosing true vegetation responses to climate [24]. Even different processing of the same NDVI data set can result in large differences in modeled global carbon fluxes [25].
The purpose of this study is to quantify how much of the interannual variability in a new NDVI data set reflects real responses of vegetation to interannual climate variability, and to evaluate whether the responses of NDVI to climate variability are meaningful enough to be of value for modeling vegetation activity. This NDVI data set is the latest version of Global Inventory Modeling and Mapping Studies (GIMMS) NDVI3g. It is derived from NOAA's Advanced Very High Resolution Radiometer (AVHRR) satellite record, and is the longest global sub-monthly time series of greenness index available (July 1981-present), more than twice as long as those available from newer sensors such as NASA's Moderate Resolution Imaging Spectroradiometer (MODIS, February 2000-present). Since this NDVI data set is used as input to or for benchmarking global carbon, water and energy cycle models, it is important that we understand how much variability in vegetation growth is captured by this NDVI data set and if the causes of variability are linked to climate (precipitation, temperature, and solar radiation), disturbance (fires and large area outbreaks of pests), human management (e.g., irrigation and fertilization), or residual errors. Other studies have been published that reveal correlations between climate and NDVI [26][27][28][29][30], but none have combined analyses on monthly anomalies, lead time dynamics, cumulative climate effects and non-climate signal interference at global scales [31,32], and none have used global NDVI time series longer than 20 years.
We explored the sources of variability in the GIMMS NDVI3g, in particular the ability of the GIMMS NDVI3g to capture vegetation responses to climate variability by: (1) evaluating the performance of GIMMS NDVI3g anomalies in comparison to the NDVI anomalies from the more advanced MODIS instrument; (2) performing correlation analysis between the anomalies of GIMMS NDVI3g and climate variables (precipitation and temperature) and comparing our results to previous studies; and (3) applying multiple linear regression to estimate the amount of GIMMS NDVI3g variance that is driven by climate versus that caused by interference in surface reflectance signal (e.g., snow, aerosols, clouds, and residual effects). Our results corroborate and expand in detail previously reported responses of NDVI to climate variability. We argue that a significant part of the variability in GIMMS NDVI3g reflects a true response of vegetation to climate variability at monthly to interannual time scales. Our analyses identified unique regional and seasonally specific dynamics in the response of vegetation to antecedent precipitation. Identification and quantification of the climate responses of vegetation represented in the GIMMS NDVI3g data will lead to improved understanding and prediction of interannual variability of, among other things, terrestrial carbon fluxes and atmospheric CO2 growth rate.

Data Sets
The GIMMS NDVI3g data set (here after GIMMS NDVI or NDVI unless specifically defined otherwise) was processed in a way consistent with and quantitatively comparable to NDVI derived from improved sensors such as MODIS and SPOT-4 Vegetation, and was corrected for dropped scan lines, navigation errors, data drop outs, edge-of-orbit composite discontinuities, and other artifacts [33]. The data processing algorithm also minimized solar zenith angle effects introduced by orbital drift and effects of changes in the sun-target-sensor geometry, corrected for bias using Sea-viewing Wide Field-of-view Sensor (SeaWiFS) data, and corrected for volcanic stratospheric aerosol effects from the El Chichon (1982)(1983)(1984) and Mt Pinatubo (1991-1993) volcanic eruptions [33,34]. Calibration and correction for orbital drift artifacts were accomplished by empirical mode decomposition (EMD); cloud contamination was filtered when the NDVI values drop beyond 95% confidence interval; and maximum value compositing (MVC) was used to minimize atmospheric (e.g., aerosols) and radiative geometry effects [33,34]. For MODIS NDVI, in addition to MVC a more sophisticated approach was used for atmospheric corrections: NDVI was derived from surface reflectances corrected for atmospheric effects [35]. The atmospheric correction algorithm of MODIS surface reflectances used an atmospheric radiative transport model and accounted for the effects of gaseous (O3, O2, CO2, and water vapor) and aerosol scattering and absorption, adjacency effects caused by variation of land cover, surface/atmosphere Bidirectional Reflectance Distribution Function (BRDF), atmosphere coupling effects, and contamination by thin cirrus [36,37].
We acquired GIMMS NDVI at the native resolution of 0.083°, bi-monthly. Since long term observational global climate data sets are at much coarser spatial scales, we focused our analysis on spatial scales of 0.25° to 1°. We aggregated the data by averaging to 0.25°, 0.5°, and 1° spatial resolutions, and by averaging bi-monthly to produce monthly time steps, depending on the native resolutions of the comparison data sets. The temporal and spatial resolutions considered here reflect those typical to global carbon, water, and energy cycle models especially those coupled to atmospheric transport. Monthly MODIS Terra (MOD13C2 V5, 2001-2010) and Aqua (MYD13C2 V5, 2003-2010) NDVI were also acquired (http://reverb.echo.nasa.gov) and regridded by averaging from 0.05° to 0.25° for comparisons.
We used the newly released monthly Global Precipitation Climatology Project (GPCP) V2.2 (http://precip.gsfc.nasa.gov). This data set has complete global coverage from January 1979 to December 2010, and a native spatial resolution of 2.5° [38]. We disaggregated the 2.5° data to 0.5° with no interpolation, and then averaged 0.5° to obtain the 1° data. We assume that the scales of significant regional climate anomalies are on the order of or larger than 2.5° × 2.5°. For comparison we also obtained Tropical Rainfall Measuring Mission (TRMM, 1998-2010) 0.25° daily precipitation product (3B42 V6, http://mirador.gsfc.nasa.gov) with nominal spatial domain of 50°N-50°S, and aggregated it from daily to bi-monthly and monthly temporal resolutions.
We constructed the temperature data set from the 0.5° Climate Research Unit (CRU CL 1.0, http://www.ipcc-data.org/obs/get_30yr_means.html) mean monthly climatology [39] and 2° Goddard Institute for Space Studies (GISS, 1200 km smoothing radius, http://www.giss.nasa.gov) monthly surface anomalies which extend for more than 130 years from 1880 to present [40]. The base period of the CRU climatology and GISS anomalies was 1961-1990 [39]. The derived temperature data set has a spatial resolution of 0.5° and a monthly temporal resolution.
Monthly MODIS snow, aerosol and cloud data for 2000-2010 were obtained to evaluate how these variables interfere with estimates of vegetation variability derived from NDVI. The snow product is MOD10CM V5 (http://reverb.echo.nasa.gov) at 0.5° spatial resolution aggregated from 0.05° resolution [41]. Aerosol and cloud data were regridded to 0.5° from the 1° MOD08_M3 V5.1 product (ftp://ladsweb.nascom.nasa.gov). All regridding from fine to coarser scales was accomplished by averaging finer data within the coarser grid cell. Regridding from coarse to finer scale was achieved by assigning the value of the coarser grid cell to the finer grid cells within. Interpolation was not used in regridding.

Methods
where y and m are the year and month of the monthly value, respectively, and N is the number of years of the data set. This procedure removes the seasonal cycles of the data sets and allows us to study the relationship between the interannual variabilities of NDVI and climate variables. The trends in the anomaly time series were not eliminated because they are considered part of the interannual variabilities, and they are strongly linked to climate change in some regions such as northern high latitudes [26]. The trends are generally very small (within ± 0.01 yr í1 , [42][43][44][45][46]) compared to the interannual variabilities in NDVI, accounting for <10% of the anomalies. For temperature, the GISS anomalies were not directly used because their base period  is not consistent with that of the NDVI data (1982-2010). All the statistical analyses (correlations and regressions) were conducted on monthly anomalies unless specified (e.g., bi-monthly for GIMMS NDVI and TRMM precipitation anomaly correlations). NDVI and precipitation data used in the analyses were GIMMS NDVI3g and GPCP precipitation, respectively, unless specified.
We conducted Pearson's correlation analyses at the grid cell level between GIMMS NDVI and MODIS Aqua NDVI anomalies at 0.25° spatial resolution for 2003-2010, between GIMMS NDVI and GPCP precipitation anomalies at 1° spatial resolution for 1982-2010, and between GIMMS NDVI and GISS temperature anomalies at 0.5° spatial resolution for 1982-2010. Selected spatial resolutions reflect the compromise between fine spatial scale NDVI signal and the coarser scale climate data. Correlation between GIMMS and MODIS Aqua NDVI was conducted on anomalies of all months (N = 96 months of 8 years), and the results were compared with the correlation between MODIS Aqua and Terra NDVI anomalies (N = 96 months of 8 years). Aqua NDVI, instead of Terra NDVI, was used to correlate with GIMMS NDVI because of the sensor degradation issue detected on the Terra platform [47]. For GIMMS NDVI-GPCP precipitation and GIMMS NDVI-GISS temperature, correlation analyses were performed on the anomalies of both all months (N = 348 months of 29 years) and individual months (N = 29 years), each with varying climate lead times to examine how fast vegetation reacts to climate anomalies. Only the correlation coefficients that were significant at 95% confidence level (p < 0.05, two-tailed) are reported.
We recognized that precipitation can only approximate soil moisture limitations on vegetation. Actual soil moisture depends not only on precipitation but also on interception, evaporation, lateral and vertical runoff, and soil characteristic such as depth and water holding capacity. To account for an energy constraint on soil moisture we analyzed correlation patterns for individual months between the anomalies of GIMMS NDVI and a water availability index defined as GPCP precipitation minus potential evaporation derived from net radiation data (2000-2009, from CERES instruments, Wielicki et al. [48], http://ceres.larc.nasa.gov/) and GISS temperature following the Priestley-Taylor method [49]. The NDVI-water availability anomaly correlation patterns (not shown) were generally the same as those for GPCP precipitation alone though the NDVI-water availability anomaly correlations were degraded by the shorter time length of the data.
We tested the assumption of Gaussian distribution of variance and possible bias caused by outliers using a non-parametric approach (Spearman's rank correlation) on the NDVI-precipitation and NDVI-temperature correlations. We found good agreement between the results obtained from Pearson's correlation and Spearman's rank correlation (see Figure S1). All other results presented are based on the Pearson's correlation.
Stepwise regression analysis was applied per grid cell to estimate the amount of NDVI variance that could be attributed to climate (antecedent precipitation and current monthly temperature), known interference for which data were available (current monthly precipitation, snow, aerosols and clouds), and a residual term that combined other unidentified sources of NDVI variability that might include disturbance, human management and underestimated errors. This analysis was done for each month to reveal the seasonality of the relationship of the variables. Only the coefficients of determination that were significant at 95% confidence level (p < 0.05) were reported.

Correlation Analyses on the NDVI-Climate Anomalies of All Months
The NDVI and climate anomalies were aggregated to a time series each (N = 348 months of 29 years), and correlation was carried out for climate lead times up to 6 months for precipitation and up to 4 months for temperature. For the lead 0 case, correlation analysis was conducted on the full time series of NDVI and precipitation or temperature anomalies. For lead k (precipitation or temperature leading NDVI by k months), we dropped the first k values from the NDVI anomaly array and the last k values from the precipitation or temperature anomaly array and conducted the correlation analysis.
We evaluated whether finer temporal and spatial resolutions influenced the strength or spatial distribution of correlations between NDVI and precipitation anomalies by comparing results from the 1°, monthly NDVI and GPCP precipitation data which we used for the remainder of our analyses with results from 0.25°, bi-monthly NDVI and TRMM precipitation. The spatial distributions of correlations were largely unaffected by the temporal and spatial scale changes though NDVI-TRMM precipitation correlations were weaker due to the shorter time span of the TRMM data.
Since we observed significant positive NDVI-precipitation correlations for multiple lead times, we further examined the correlations between NDVI anomaly and precipitation anomaly cumulated for varying lead time durations. We applied correlation analysis to NDVI anomaly of the current month and precipitation anomaly summed over a duration ranging from 1 to 10 months (cumulative precipitation anomaly hereafter, N = 348 -duration + 1). Details about correlation analysis for varying durations can be found in Wang et al. [50]. Briefly, for duration of k months, the cumulative precipitation anomaly was the sum of precipitation anomalies over a k-month period corresponding to lead = 0 to k-1 months. The results of these analyses were presented in detail later but they showed negative correlations between NDVI anomalies and precipitation anomalies at 1-month duration (i.e., lead 0) which were not ecologically meaningful and were likely the result of cloud interference and were thus not used in further analyses. The highest correlations were found between anomalies in prior 6 months cumulative precipitation and NDVI. From these results we identified 6-month duration from lead 1 to 6 months (lead 1 + 2 + , … , + 6) to be the optimum duration to represent the total effects of precipitation on vegetation growth, and this 6-month duration was applied to correlation analyses on NDVI-cumulative precipitation correlations for individual months in Section 2.2.2.
When dealing with correlations between two time series (e.g., the NDVI and precipitation anomalies of all months) we have to account for the effect of temporal autocorrelation (also called serial correlation), otherwise the significance of correlation would appear to be higher than it actually is (i.e., giving a smaller p value). Ordinary correlation assumes that the observations of a variable are independent of each other. However, this assumption is generally invalid for time series. In a time series, very often the observation at time i is closely associated with the observations at time i í 1, i í 2, … , meaning that temporal autocorrelation exists in the time series. For example, NDVI anomaly in June generally does not deviate much from NDVI anomaly in May. For a time series with N observations Xi (i = 1, 2, … , N), if there is a significant correlation between the subsets Xi (i = 1, 2, … , Ník) and Xi+k (i + k = 1+k, 2+k, … , N), this correlation is called the kth-order temporal autocorrelation and the associated correlation coefficient is the kth-order temporal autocorrelation coefficient of this time series [51]. Few previous studies have corrected their correlation results for temporal autocorrelation. Without such correction, the strength of the correlations would be overestimated and the interpretation of the results would be misleading [51].
For our correlation analyses between NDVI and climate variable anomalies of all months, we accounted for the effect of the 1st-order temporal autocorrelation following the procedure by Dawdy and Matalas [51]. Ideally higher orders of temporal autocorrelation should also be accounted for until the temporal autocorrelation of at least one of the two time series becomes insignificant. However, this process is very laborious given the large number of land grid cells. Practically only the 1st-order temporal autocorrelation is considered [52]. This is reasonable because for NDVI and climate variables the correlation between values of current month and last month is likely the strongest [53].
Briefly, 3 steps [51] were applied to each land grid cell to estimate the correlation coefficient (r value) and significance of the correlation coefficient (p* value) that is corrected for the 1st-order temporal autocorrelation. Take the GIMMS NDVI-GPCP precipitation anomaly correlation at lead 0 as an example, so here the two time series are 1982-2010 monthly NDVI and precipitation anomalies (N = 348). First, the correlation coefficient (rNDVI-precipitation) and the associated p value were calculated. Second, the 1st-order temporal autocorrelation coefficients rNDVI (1) and rprecipitation(1) of the NDVI anomaly time series and the precipitation anomaly time series, respectively, were computed. Finally, if both rNDVI(1) and rprecipitation(1) are significant at 95% confidence level, we derived the effective number of degrees of freedom (n*) and the adjusted p value (p*) from the sample size N and the correlation coefficients above using Equations (2)-(4): here n = N í 2, and f(t) is the probability density function of Student's t-distribution. If either rNDVI (1) or rprecipitation(1) is insignificant at 95% confidence level, no correction is needed and p* is the same as p. The NDVI-precipitation anomaly correlation is considered significant if the p* value is below 0.05. In summary, the correlation coefficient rNDVI-precipitation remains the same but the significance of the correlation coefficient is degraded (i.e., p* > p) during this process.

Correlation Analyses on the NDVI-Climate Anomalies of Individual Months
Correlation analysis was also conducted for each month of the year (January to December) for climate lead times up to 7 months for precipitation and up to 4 months for temperature (N = 29). Correlations generally deteriorated after ~4 months, but in some places NDVI-precipitation correlations remained significant positive for lead times up to 6 months. Therefore, similar to the correlation of all-month data, for each month we accounted for the long-lasting effect of precipitation on NDVI by further analyzing the correlations between NDVI anomalies and precipitation anomalies cumulated over a 6-month duration from lead 1 to 6 months (N = 29), and the correlation results will be used in Section 2.2.3 to estimate the fraction of NDVI variance driven by climate.
The NDVI-temperature correlations of northern high latitudes need to be interpreted with caution for the winter and spring months due to the interference of snow. In warmer winters and springs, there is less snow cover and the land surface produces a higher NDVI signal. However, this does not necessary mean a higher rate of vegetation growth. We detected potential snow interference in the NDVI-temperature correlations when NDVI-snow correlations were significantly negative (p < 0.05, N = 132 for all months and N = 11 for individual months). For both all months and individual months we corrected snow interference by setting the NDVI-temperature correlation of a grid cell to zero if significant negative NDVI-snow correlation was observed for this grid cell.

NDVI Variance Explained by Climate
The fraction of the NDVI variance explained by climate (R 2 ) can be approximately computed from the correlation coefficients (calculated in Section 2.2.2) of NDVI-cumulative precipitation (lead 1-6) and NDVI-temperature at lead 0 corrected for snow interference. NDVI-cumulative precipitation correlation was used to account for the long-lasting effect of precipitation on vegetation, and NDVI-temperature correlation at lead 0 was selected because it's the strongest among all the lead times we tested (from 0 to 4 months).
For each land grid cell in each month, if only one of the two correlation coefficients (cumulative precipitation or temperature) was significant positive, R 2 was the square of the significant positive correlation coefficient; if both were significant positive, R 2 was essentially the coefficient of determination of the multiple linear regression of NDVI on cumulative precipitation and temperature and could be calculated using the following equation [ where r is the correlation coefficient. For all other cases, R 2 was zero. Only significant positive correlations were considered here because they were ecologically meaningful: increased vegetation activity with warmer temperatures in mid-and high-latitudes and with higher rainfall in water limited regions.

NDVI Variance Explained by Accounting for Atmospheric Interference
After estimating the variance in NDVI caused by climate variability we briefly addressed the variance not caused by climate. Snow, aerosols and clouds among others can interfere with detection of vegetation greenness by satellites. The GIMMS NDVI has been adjusted to account for many of these effects, but such effects have not been fully eliminated. The need for snow correction in the interpretation of the temperature correlations is one example of interfering elements that might contribute to NDVI variance. Precipitation of current month was considered as another potential interfering element because it was negatively correlated with NDVI in many regions (details will be presented in Section 3.2) which is not ecologically meaningful.
To examine how much of the NDVI variance was explained by the effects of snow, aerosols, clouds and precipitation of the current month, forward stepwise multiple linear regression analysis [55] was applied to each grid cell and each of the 12 months. NDVI anomaly was the dependent variable, and anomalies of snow, aerosols, clouds and precipitation of the current month were independent variables (N = 11, 2000-2010). All of the independent variables included in the regression models were negatively correlated with NDVI anomalies. More than 85% of the final regression models (p < 0.05) were found to have only one independent variable, so the multicollinearity problem was negligible and no correction was needed.

NDVI Variance Unexplained
Even after taking into account climate and the interference identified above, there was still some of the NDVI variance that was not explained. To determine the extent of this fraction, forward stepwise multiple linear regression analysis was employed to the anomalies of NDVI, cumulative precipitation (lead 1-6), temperature, snow, aerosols, clouds and precipitation of the current month (lead 0) on all land grid cells where data were available, with NDVI anomaly as the dependent variable and the rest as independent variables (N = 11, 2000-2010). The coefficient of determination (R 2 ) of the regression represents the percentage of NDVI variance driven by climate anomalies and atmospheric interference listed above, and the fraction of the NDVI variance unexplained was estimated as 1 í R 2 . More than 79% of the final regression models (p < 0.05) had only one independent variable, so again the multicollinearity problem was negligible.

Comparisons with MODIS
Though the GIMMS AVHRR data span almost three times the length of MODIS, the latter was designed to measure specific narrow reflectance bands, to have explicit atmospheric corrections, to be highly calibrated, and to minimize sun-surface-platform radiative geometry effects. Many design features of MODIS were implemented to overcome problems originally identified in the AVHRR record. The new version of the GIMMS NDVI data set represents the culmination of 3 decades of experience with AVHRR NDVI focused on identifying and minimizing known problems. To gain insight about the robustness of the GIMMS NDVI for analysis of its longer time period we compared GIMMS and MODIS monthly NDVI anomalies for the 2003-2010 Aqua period. Seventy six percent (76%) of land grid cells showed statistically significant (p < 0.05) positive correlations between GIMMS NDVI and MODIS Aqua NDVI anomalies ( Figure S2a). This number was close to the percentage (87%) of land grid cells that had statistically significant (p < 0.05) positive correlations between the anomalies of MODIS Aqua and Terra NDVI ( Figure S2b). Correlations between GIMMS and MODIS Aqua were generally high except over tropical forests, eastern temperate regions and a band across boreal Asia ( Figure S2a). These same regions also showed lower correlations between Terra and Aqua NDVI ( Figure S2b), indicating high uncertainty in the interannual NDVI signal in these regions for all sensors. Gallo et al. [56] and Fensholt and Proud [43] reported similar agreement between older versions of GIMMS NDVI and MODIS NDVI.
These comparisons build confidence that at the scales of our analysis the GIMMS NDVI is comparable to the more advanced but shorter time length MODIS data. GIMMS NDVI should perform nearly as well as MODIS in detecting climate driven variability but in the case of GIMMS NDVI over the longer 30.5-year record.

Precipitation Controls on NDVI Variability
To characterize the response of vegetation measured as NDVI to precipitation we examined what regions and what latencies occur in the correlation between NDVI and precipitation. For each grid cell we aggregated all monthly NDVI anomalies (N = 348, January 1982 to December 2010) and identified significant correlations with precipitation anomalies at various precipitation anomaly lead times. We found that globally positive correlations were the strongest and most spatially extensive at lead 1 (Figure 1a (Figure 1d, [59]). Generally NDVI-precipitation anomaly correlations were strongest for grasslands, shrublands, croplands, and savannas and lowest for forests. The lower correlations between NDVI and precipitation in forests probably reflect in part lower water limitations compared to the regions of high correlations.    NDVI correlations with cumulative precipitation summed over various lead times were examined to account for the long lasting effects of precipitation. Figure 2a-c shows the spatial patterns in correlations with cumulative precipitation of all months (N = 348 í duration + 1) for durations of 1 to 3 months (see Figure S3 for correlations for durations up to 10 months). The number of global land grid cells with significant positive correlations increased markedly from 1-month duration (lead 0) to 2-month duration (lead 0 + 1) and leveled off at 7-month duration (lead 0 + 1 + , … , + 6), and 80% of the land grid cells that have significant positive correlations at 7-month duration were significant at 2-month duration (Figure 2d). The largest increase in the number of pixels with significant positive correlations occurred at lead 1. Cumulative precipitation anomalies at 1-month duration (lead 0) were negatively correlated with NDVI anomalies in some boreal regions and in part of Southeast Asia and South America (Figure 2a). Such negative correlations are not ecologically meaningful and thus precipitation of the current month was treated as an interfering element in the NDVI signal in Section 3.5.
Results from lead correlation analyses for individual months showed that globally the pattern of strongest correlations at lead 1 was valid for all 12 months of the year, and the strongest correlations generally occurred during each region's respective rainy season (Figure 3). For example, in central North America, central Eurasia, India and the Sahel, the strongest correlations with NDVI anomalies in June, July, August, September and October were for precipitation anomalies in May, June, July, August, and September, respectively (Figure 3f-j). Another example is eastern and southern Africa and Australia, where NDVI anomalies in each month of December to April were most strongly correlated with precipitation anomalies 1 month before (Figure 3l and 3a-d). Figure 4 shows the percentage of global land grid cells that have positive correlations (p < 0.05) between NDVI and precipitation anomalies as a function of lead times from 0 to 6 months for each month of the year. Often correlations were negative at lead 0 for some areas (Figure 2a), but strong positive correlations emerged at leads of 1 and 2 months and peaked at lead 1 with 8-12% of total land grid cells having significant positive NDVI-precipitation correlations for all 12 months of the year (Figure 4). Longer leads showed diminishing correlations. These global patterns are consistent with regional studies using either NDVI or enhanced vegetation index (EVI) and either ground-or satellite-based climate data [50,[60][61][62], and with the global analysis of Lotsch et al. [32].
Despite the general global pattern of strongest correlations at 1-month lead, in some temperate and subtropical herbaceous dominated regions NDVI responses to precipitation lead times varied seasonally, reflecting the dynamics of soil moisture availability. Prominent examples are described below.
Correlation patterns over western and central North America exhibited seasonally varying dependence of NDVI on precipitation (Figures 5a and 3e-h): NDVI anomalies in May were most strongly correlated with winter precipitation anomalies; in June and July significant correlations were observed between NDVI and recent precipitation with some residual dependence of NDVI on winter/spring precipitation; by August only most recent precipitation history showed signification correlations with NDVI. These results are consistent with expected soil moisture storage dynamics: spring and early summer plant growth is dependent upon precipitation received during winter when precipitation exceeds evaporation and the excess precipitation is stored in snow and soils, whereas in late summer soil moisture from spring is largely depleted and thus plant growth depends on most recent precipitation. Although not pointed out in Méndez-Barroso et al. [61], similar patterns can also be found in their results from the station-level time series of absolute EVI and precipitation values in northwestern Mexico which is among the regions presented here. An example of a different seasonal dynamic in NDVI responses to precipitation was seen in sub-tropical regions such as Australia. During the rainy season correlations were strongest for recent precipitation (October-May, Austral spring-fall), while in all dry season months (e.g., June-August, Austral winter) NDVI anomalies were correlated most strongly with previous rainy season precipitation (Figure 5b, Figure 3 and Figure S4). In other words, correlations with NDVI anomalies during the dry season were highest for precipitation in the previous February-April wet season. Another example of this pattern was found in Southern Africa ( Figure S4).

Temperature Controls on NDVI Variability
Combining all monthly NDVI and temperature anomalies (N = 348, January 1982 to December 2010) produced positive correlations between NDVI and temperature anomalies in northern mid-and high-latitudes ( Figure 6). Zhou et al. [63] reported similar results in their northern latitude study. Correlations with temperature were highest at lead 0 ( Figure 6). Central Europe shows the strongest temperature effect (Figure 6a) that lasts for 3 to 4 months (Figure 7), while in the other northern mid-and high-latitude regions the temperature effect lasts for only 1 to 2 months (Figure 7). This is likely a result of the much shorter snow cover duration in central Europe compared to other regions in the same or higher latitudes. MODIS Terra snow cover data show that in March snow cover has almost completely disappeared in central Europe while in other regions in the same or higher latitudes snow cover remains above 70%. The shorter snow cover duration in central Europe is also reflected in the earlier start of the growing season in central Europe compared to other regions in the same or higher latitudes (Figure 8).
Seasonal positive correlations between NDVI and current temperature anomalies were evident and supported ecological interpretations. Positive correlations emerged in spring (March) and progressed northward into June (Figure 7), consistent with the field observations in northeastern Siberian tundra by Blok et al. [64] showing that early summer temperature is the most important factor influencing vegetation growth in mid-and high-latitudes. The spring-early-summer northward wave of correlation matches closely with the progress of the MODIS-derived start of the growing season (MOD12Q2, [65,66], Figure 8). A similar but weaker wave of correlation moved south from mid-summer to early fall (July-September, Figure 7).

Climate Driven Variance in NDVI
As expected, we found that precipitation was rarely negatively correlated (p < 0.05) with NDVI for leads of 1 month or more. Temperature was always positively correlated with NDVI anomalies at a lead of 0 month except where temperature and precipitation anomalies were negatively correlated. The latter were water limited regions and low precipitation was associated with higher temperature and drought. Combining the variance explained by positive correlations with cumulative antecedent precipitation (1-6 months lead) and by positive correlations with temperature at lead 0 (excluding negatively correlated NDVI-snow grid cell-months) following the method described in Section 2.2.3 and Equation (5), we estimated the variance in monthly NDVI explained by climate (Figure 9). Percentage of monthly NDVI variance (R 2 ) driven by climate ranged from about 30% to 70% in some boreal regions of Eurasia and North America in spring and summer due to temperature, and in central and western North American, central Eurasia, Australia, Africa, and South America due to precipitation.
Our approach to characterizing climate effects on NDVI was simple in that it considered only monthly precipitation and temperature (though we have looked at bi-monthly precipitation and NDVI and found little difference in results) and it used a single variant approach. Inclusion of the Priestley-Taylor energy constraint did not improve correlations with NDVI compared to precipitation-NDVI correlations. We have also applied a multi-variant approach but the results did not change our overall conclusions and were more difficult to interpret.
Our simple statistical approach was applicable because the regions influenced by temperature versus precipitation tended to be distinct: temperature correlations were strongest above 30°N latitude while precipitation correlations were strongest for herbaceous vegetation, < 1000 mm per year precipitation below 60°N latitude. In addition we evaluated the correlations between temperature and precipitation (not shown) which showed little interactions in terms of effects on NDVI though the analysis helped to explain noise and bias in the NDVI.

Influence of Snow, Aerosols, Clouds and Precipitation of the Current Month on NDVI Variability
Some of the important conditions that interfere with satellite measurements of vegetation activity using NDVI include snow cover, atmospheric aerosols, and clouds. Precipitation of the current month showed negative correlations with NDVI anomalies in some regions ( Figure 2a) and was also treated as an interfering element here. The readily available MODIS snow, aerosol and cloud products and precipitation data motivated exploration into how much variance in the NDVI signal is caused by these interfering elements that linger in the GIMMS NDVI3g. Figure 10 is an example month (September) from the series of monthly maps in Figure S5 that compare total monthly NDVI variance (Figures 10a and S5a), fractions of NDVI variance contributed by climate (Figures 10b and S5b, same as Figure 9g) and by atmospheric interference (Figures 10c and S5c), and fraction of unexplained NDVI variance (Figures 10d and S5d). Atmospheric interference includes snow, aerosols, clouds and precipitation of the current month all of which are negatively correlated with NDVI at lead 0 month (Section 2.2.4). Unexplained NDVI variance represents the residual NDVI variance that is not significantly correlated with cumulative precipitation (1-6 months lead) or current (lead 0) temperature, snow, aerosols, clouds, and precipitation anomalies (Section 2.2.5). Generally NDVI variance was <0.01 (Figures 10a and S5a), which translates to a standard deviation of <0.1. In other words NDVI anomalies were less than 0.1 for most land grid cells most of time. Figure 10c suggests that interference by snow, aerosols, clouds and precipitation remains in the NDVI signal after the GIMMS NDVI3g processing. For some regions and some periods (e.g., North America, Europe and central Asia in winter months), more than 80% of the NDVI variance was associated with these atmospheric and surface interfering elements. Some notable patterns in the monthly maps in Figure S5 include (1) snow in April-May and climate (temperature) in June explained a large fraction of NDVI variance at high northern latitudes; and (2) the variability explained in the Amazon was mostly from clouds in February during the rainy season and from aerosols in September caused by fires (Figure 10c).   Even the combination of climate and the interference listed above cannot fully explain all the NDVI variance. In regions such as Europe, southeastern China, parts of South America, and sometimes parts of boreal Asia, a significant fraction of NDVI variance remains unexplained (Figures 10d and S5d). Figure 11 summarizes regional (regions as defined in Transcom [67]) and seasonal patterns in % of land grid cells in which NDVI anomalies (1) are significantly correlated with climate (cumulative precipitation of lead 1-6 months and temperature at lead 0, all positive), (2) are accounted for variance (from snow, aerosols, clouds and precipitation of the current month combined, all negatively correlated with NDVI anomalies), and (3) show no correlations with the variables tested in this study. Significant climate correlations are generally most prevalent during the growing season (e.g., boreal and temperate regions, Southern Africa). Fraction of land grid cells in tropical regions (Tropical America and Asia) showing significant positive correlations with climate is low with little seasonal variability. Australia had high correlations with climate (precipitation) for all 12 months. Fraction of land grid cells with no significant climate or interference correlations was high in non-growing season in boreal and temperate regions. Tropical regions and temperate South America had high fractions of land grid cells with no significant correlations with the variables tested. These results (Figures 9-11 and S5) show when and where NDVI was more reliable at capturing vegetation responses to climate anomalies.

Uncertainties in the Data and Results
All correlations and regressions presented here were those with probabilities >0.95 that were obtained solely by random chance. We tested the assumption of normality of distributions and the potential for bias caused by outliers by performing non-parametric Spearman's rank correlation and found little influence on the results and conclusions. We tested for potential bias as a result of temporal and spatial aggregation by using a different precipitation data set (TRMM) with finer temporal and spatial resolutions (bi-monthly and 0.25° as opposed to monthly and 2.5° for original GPCP) and found that it did not change the results or conclusions. We tested to a limited extent the assumption that precipitation was an adequate surrogate for water limitation and found the assumption to be supported. We identified regions and months where correlations were most robust as well as where/when uncertainty was large. We evaluated the GIMMS NDVI against two other satellite data (MODIS Aqua and Terra) and identified where all three products were in agreement and where satellite NDVI data in general were not adequate to address vegetation responses to climate.
There was still a significant fraction of NDVI variance unexplained in some regions and some periods. The unexplained NDVI variance may be partly due to the uncertainty in the NDVI and climate data sets. It could, however, be a result of the variability in other drivers of NDVI that we didn't account for in this study (see the next section).

Significant Drivers of NDVI Variance
Other than climate driven variability, we have not taken into account processes that are likely to also influence NDVI variability, such as stimulation of plant growth caused by increasing atmospheric CO2, land use change such as deforestation and biomass burning, irrigation and fertilization, and N deposition. Particularly, land management practices such as irrigation and fertilization contribute to a significant fraction of the NDVI variability in croplands [68].
From our analysis we identify a number of significant causes of NDVI variance and their likely effects on NDVI.
(1) Climate driven variability in NDVI as quantified in this study. Northern latitude NDVI anomalies were positively correlated with temperature. NDVI anomalies in temperate to tropical semi-arid and arid regions were positively correlated with precipitation with regionally varying seasonal dynamics. (2) Climate driven variance missed by our simple representation of climate variability (e.g., soil moisture, solar radiation, vapor pressure deficit (VPD), and evaporation) and by our linear spatial averaging. (3) Atmospheric interference from snow, aerosols and clouds reported here are negatively correlated with NDVI. Though these are addressed in part by the GIMMS NDVI3g processing, we found some residual negative correlations with these variables.
(4) Variance caused by sensor-sun-surface radiation geometry. This uncertainty is also addressed as part of the GIMMS NDVI3g processing prior to our use. (5) Other mechanistic processes not accounted for: irrigation, fertilization, fires and land degradation. (6) Atmospheric and surface interfering factors not accounted for in our analysis (e.g., water vapor and cloud type) and from aerosol, water vapor, and cloud effects not captured in the GIMMS NDVI3g processing.

Future Work
Our results support the use of satellite data for the representation of vegetation responses to climate variability in models. A few studies have confirmed that vegetation anomalies observed from satellites are associated with regional climate anomalies and are large enough to explain the anomalies in measured and modeled carbon and water fluxes [69,70]. However, some global scale studies show little correspondence between satellite-based estimates of FPAR (fraction of photosynthetically active radiation) and anomalies in global carbon fluxes from atmospheric inversions [25,71]. Our next work will be to characterize the impacts of climate driven NDVI anomalies on modeled carbon fluxes and to evaluate these impacts using independent methods such as eddy covariance data and atmospheric CO2 inversions. We expect that variances of up to 0.01 in the NDVI signal in some regions during growing season can represent relative standard deviations of up to 10% in the FPAR which itself is proportional to NPP in some global carbon cycle models [17]. Preliminary evaluation of how this amount of climate driven variability influences carbon fluxes (not shown) suggests that climate driven variability in NDVI identified in this study could produce variability in simulated NPP of 20% or more in some regions and seasons. Examples of this level of variability in gross primary production and leaf area index have been reported at eddy covariance sites [72,73], and further work evaluating model responses to climate driven NDVI variability against other independent flux data is underway.

Conclusions
Our analyses confirm that GIMMS NDVI3g (Global Inventory Modeling and Mapping Studies normalized difference vegetation index, version 3) captures vegetation responses to climate nearly as well as more advanced sensors such as MODIS (Moderate Resolution Imaging Spectroradiometer) Aqua, in agreement with results from other studies [43,56]. The 30.5-year GIMMS AVHRR (Advanced Very High Resolution Radiometer) NDVI3g record thus provides a robust time series to explore causes of variations in land surface vegetation greenness as expressed by NDVI.
In addition to taking advantage of the longer time series we also provided a more thorough month-by-month analysis of drivers of NDVI variability than previous global studies. Temperate and subtropical arid and semi-arid regions and northern high latitude spring and summer showed strong interannual climate driven variability in NDVI. Seasonal variability in the importance of antecedent precipitation represented seasonal variability in available soil moisture as in central North America where dependence of vegetation growth on water stored from the winter gradually diminishes into the summer. These patterns of variability could be used to evaluate soil moisture predictions from more sophisticated, complex soil moisture models. As expected, evergreen systems showed reduced seasonal and interannual variations. This behavior was true for precipitation, but our results showed that 50% or more of the NDVI variance in May can be attributed to temperature variability (Figure 9c) for some regions in the boreal forests. Thirty to sixty percent (30-60%) of the NDVI variance for the Sahel-Sahara border of northern Africa and in Southern Africa can be explained by precipitation variability during their respective rainy seasons. During their respective dry seasons NDVI variance was low in general. Temperate grasslands showed strong climate (precipitation) controls during the growing season. Here again NDVI variance tended to be low outside the growing season so low explained variance by climate was expected for these periods. The climate driven variability in NDVI shown in this study demonstrates the ability of the GIMMS NDVI3g data record to provide realistic representation of interannual variability in solar radiation absorbed by vegetation.