Recent Changes in Terrestrial Gross Primary Productivity in Asia from 1982 to 2011

Past changes in gross primary productivity (GPP) were assessed using historical satellite observations based on the Normalized Difference Vegetation Index (NDVI) from the Advanced Very High Resolution Radiometer (AVHRR) onboard the National Oceanic and Atmospheric Administration (NOAA) satellite series and four terrestrial biosphere models to identify the trends and driving mechanisms related to GPP and NDVI in Asia. A satellite-based time-series data analysis showed that approximately 40% of the area has experienced a significant increase in the NDVI, while only a few areas have experienced a OPEN ACCESS Remote Sens. 2013, 5 6044 significant decreasing trend over the last 30 years. The increases in the NDVI are dominant in the sub-continental regions of Siberia, East Asia, and India. Simulations using the terrestrial biosphere models also showed significant increases in GPP, similar to the results for the NDVI, in boreal and temperate regions. A modeled sensitivity analysis showed that the increases in GPP are explained by increased temperature and precipitation in Siberia. Precipitation, solar radiation and CO2 fertilization are important factors in the tropical regions. However, the relative contributions of each factor to GPP changes are different among the models.


Introduction
Asia, which is characterized by a rapidly growing economy, human-induced disturbances, and climate changes, is an important region for understanding global terrestrial carbon cycles [1][2][3].Siberian regions have experienced a 1.0-2.0°C increase in temperature from 1961 to 2004 [4].China and India have recently experienced rapid economic growth and a large increase in CO 2 emissions.For example, CO 2 emissions in China and India accounted for 28% and 7%, respectively, of the global CO 2 emissions in 2011 [5].In addition, many regions, for example, China, Southeast Asia, and South Asia, have experienced land-use changes, including changes in the use of cropland and forest areas [6,7].
Terrestrial gross primary productivity (GPP) is the largest atmosphere-land carbon flux, driving terrestrial carbon cycles as a primary input for vegetation carbon.An accurate estimate of the magnitude of temporal GPP changes is required to quantify and predict terrestrial carbon budgets.Numerous approaches have been used to analyze GPP at both global and regional scales, which include analyzing satellite-based remote sensing data as an indicator of GPP [8][9][10][11], using empirical upscaling with machine learning techniques [12][13][14][15], and assembling diagnostic [16][17][18] and prognostic [19,20] terrestrial biosphere model simulations.
These approaches have also been used to analyze spatio-temporal GPP variations in Asia.For example, numerous studies have analyzed satellite-based time-series data, including data that from the Advanced Very High Resolution Radiometer (AVHRR) onboard the National Oceanic and Atmospheric Administration (NOAA) satellite series, to infer GPP variations in relation to environmental changes.Analyses of NOAA AVHRR data have revealed that vegetation productivity is increasing and the growing season has lengthened in boreal and Arctic regions [8,9,21].In East Asia, including China, several studies have revealed that the vegetation greenness has increased over the last 20-30 years [22][23][24][25][26].In Southeast Asia, interannual variations in vegetation indices were found to be largely determined by the El Niño Southern Oscillation [27].Furthermore, vegetation greenness was found to have increased in India because irrigated areas and fertilization use increased [28].However, this increase in vegetation greenness has stagnated in recent years [29].Terrestrial biosphere models and empirical models have also been used to analyze GPP variations in many regions in Asia.Empirically upscaled estimates of GPP in East Asia revealed a significant reduction in the summer of 2003 in response to meteorological anomalies [30].Increased GPP has also been identified using diagnostic terrestrial biosphere models in China [31] and India [28].
Further analyses could utilize a combined application of satellite-based observations and terrestrial biosphere models.Moreover, satellite-based observation data can be used for model evaluation.Once the model output successfully reproduces the satellite-based observations, mechanisms from the observed results can be analyzed with process-based models using sensitivity experiments.Recent vegetation greening trends in northern high-latitude regions have been analyzed using the satellite-based leaf area index (LAI) and a terrestrial biosphere model; increases in temperature were related to the vegetation greening trend [32].The combined use of remote sensing data and terrestrial biosphere models have helped explain the observed changes in the Normalized Difference Vegetation Index (NDVI) in the Amazon from 1984 to 1999.In this region, increased solar radiation has been found to enhance the terrestrial net primary productivity [33].Vegetation growth in the Northern Hemisphere (>25°N) has been explained by a process-based terrestrial biosphere model and has been found to have resulted from changes in atmospheric CO 2 , temperature, and precipitation [34].An application of the above approaches to Asian regions could elucidate the changes in vegetation greening trends in Asia in a process-based manner.
In this study, we used both satellite-based time-series observations and four process-based terrestrial biosphere models to identify and understand changes in terrestrial GPP in Asia (the regions are defined in Figure 1).Using the updated satellite-based vegetation index products, that is, NOAA AVHRR NDVI3g [37,38], we first quantified changes in the NDVI from 1982 to 2011 in Asia and identified the regions with increasing or decreasing trends in vegetation activity.Then, we applied four different terrestrial biosphere models to simulate recent changes in terrestrial GPP; we verified these results with satellite-based NDVI variations as a proxy of terrestrial GPP at sub-continental scales.Finally, using the model simulation results, we determined the causes of the simulated trends in terrestrial GPP.

Study Region
We analyzed a region that covers a wide range of vegetation and climate conditions in Asia (Figure 1a), including a temperature gradient from tropical forest regions in Southeast Asia (annual average temperature > 25 °C) to tundra regions in Siberia (annual average temperature < −10 °C) as well as a precipitation gradient from coastal and tropical regions (>2,000 mm•year −1 ) to semiarid or arid regions in Central Asia (<200 mm•year −1 ).The vegetation cover also exhibits large differences with evergreen tropical forests in Southeast Asia; temperate forests in East Asia; boreal forests in Siberia; croplands in China, India and Russia; grasslands in central Asia; and open shrublands (tundra) in northern Siberia.
For a sub-continental scale analysis, we divided the region into nine sub-regions based on the definition used in the atmospheric CO 2 inversion analysis (Figure 1b).We adopted the sub-continental division of the Greenhouse Gases Observing Satellite (GOSAT) Level 4A products (i.e., CO 2 flux for each of the 64 global regions) [36].The areal definition is based on the continental-scale vegetation type.Therefore, this sub-region definition is suitable for assessing the large-scale terrestrial carbon cycle changes in specific ecosystems at the sub-continental scale.

Satellite-Based Time-Series Data: NDVI, LAI, and FPAR
The latest version of the NDVI dataset generated from the AVHRR sensor onboard a series of NOAA satellites (Global Inventory Modeling and Mapping Studies (GIMMS) NDVI3g) from January 1982 to December 2011 was used in this study [37,38].The NDVI data were used as a proxy for terrestrial GPP [39] to compare with the model outputs.The datasets were provided with a spatial resolution of 8 km × 8 km and a 15-day temporal resolution over the period from July 1981 to December 2011.To create spatially and temporally consistent data with the terrestrial biosphere model simulations, we averaged the GIMMS NDVI3g data to create quarter-degree spatial resolution and monthly temporal resolution data.
The third generation data of the LAI and Fraction of Photosynthetically Active Radiation absorbed by a canopy (FPAR) datasets (LAI3g and FPAR3g) [40] were used as inputs for the BEAMS model [18,41].These datasets are integrated products of GIMMS NDVI3g data [37,38] and Moderate Resolution Imaging Spectroradiometer (MODIS) LAI/FPAR products [42] that include updated algorithms for data screening [43].In generating the LAI3g and FPAR3g data, a set of neural networks were first trained on best-quality and post-processed MODIS LAI and FPAR products as well as GIMMS NDVI3g data for the overlapping period from 2000 to 2009.The trained neural networks were then used to produce the LAI3g and FPAR3g datasets using GIMMS NDVI3g data for the remaining period.The original data were provided at 15-day temporal intervals and 8-km spatial resolutions.We converted these data into monthly temporal resolution and quarter-degree spatial resolution by averaging that data using the same procedure used for the NDVI dataset.

Climate Time-Series Data
Daily climate data were generated from 1901 to 2011 with quarter-degree spatial resolution by merging reanalysis-based products and observation-based products based on a similar approach used with the CRUNCEP data [44].For example, daily air temperature data from 1901 to 2011 were generated based on the Climate Research Unit Time Series 3.2 (CRU TS3.2) data (available period: 1901-2011) [45] and National Centers for Environmental Prediction (NCEP)/the National Center for Atmospheric Research (NCAR) reanalysis (NCEP/NCAR reanalysis) data (available period: 1948 to present) [46].The daily variations in air temperature were based on NCEP/NCAR reanalysis data; the monthly averages were adjusted to fit the CRU TS3.2 data by adding an offset.Because NCEP/NCAR reanalysis data are not available for the 1901-1947 period, the reanalysis data for 1948 were used instead.Similarly, precipitation data were generated by merging daily NCEP/NCAR reanalysis data and CRU TS3.2 monthly precipitation data [45] by multiplying by an adjustment factor.In addition, we found that the frequency of rainy days in the NCEP/NCAR reanalysis data was too high compared with the CRU TS3.2 wet-day frequency data.Therefore, the frequency of wet days was also adjusted using the CRU TS3.2 wet-day frequency data by redistributing the rain.Incoming surface shortwave radiation data were generated by merging daily shortwave radiation created by the Mountain Climate Simulator (MTCLIM) [47] (available period: 1901-2011) and monthly averages of Japan Aerospace Exploration Agency (JAXA) NASA/MODIS global 5-km irradiance (JAXA radiation) data [30,48] (available period: 2000-2010) by multiplying by an adjustment factor.The MTCLIM-based shortwave radiation relies on the assumption that variations in diurnal temperature ranges are related to incoming solar radiation.Because JAXA radiation data are not available for 1901-1999 and 2011, the average of the data for 2000-2010 was used instead.Daily vapor pressure deficit data were generated using MTCLIM.For longwave radiation and wind speed, NCEP/NCAR reanalysis data were directly used; the data from 1948 were used for the 1901-1947 period.

Atmospheric CO 2 Time-Series Data
Atmospheric CO 2 data were obtained from the observed monthly mean at the Mauna Loa station from 1959 to 2011 [49] and from the observed annual mean from the ice core at Law Dome, Antarctica, from 1901 to 1958 [50].The monthly mean atmospheric CO 2 concentration data from 1901 to 1958 were created by subtracting differences in the annual mean between 1959 and the corresponding year from the monthly data in 1959.The annual mean atmospheric CO 2 concentrations were used as inputs to the Biome-BGC, LPJ and TRIFFID models.The monthly mean atmospheric CO 2 concentrations were used as model inputs for BEAMS.

Static Data
Static data used for the analysis included land cover, elevation, soil texture, and soil depth.Land-cover data were derived from MODIS products (MCD12Q1 collection 5.0) [35] by taking the most dominant land cover in each quarter-degree grid cell.Elevation data were obtained from GTOPO30 [51].Soil texture data and soil depth data were created using the Harmonized World Soil Database version 1.2 and ISLSCP II Ecosystem Rooting Depths (mean 95% ecosystem rooting depth in meters [52]), respectively.All data were re-gridded to quarter-degree resolution using the nearest neighbor method.

Terrestrial Biosphere Models
We used four terrestrial biosphere models to analyze changes in GPP from 1982 to 2011.The models included BEAMS Version 1.3 [18,41], Biome-BGC Version 4.2 [53], LPJ Version 3 [20,54], and MOSES/TRIFFID [55,56].These models differ in terms of their model structure, parameters, and required input data.The BEAMS model is categorized as a diagnostic terrestrial biosphere model, which simulates spatial and temporal patterns in the carbon cycle more realistically using satellite-based LAI and FPAR inputs than prognostic terrestrial biosphere models.The model calculates GPP using a function of FPAR, shortwave radiation and an environmental regulation factor.The environmental regulation factor is calculated by combined Farquhar's photosynthesis model [57,58] and the Ball-Berry-type stomatal conductance model [59,60] using CO 2 , temperature, soil moisture, and relative humidity.The Biome-BGC model is a prognostic terrestrial biosphere model with fixed vegetation cover.The model runs with only climate data as a time-variant input.The model calculates GPP using a function of the CO 2 , shortwave radiation, temperature, soil moisture, and vapor pressure deficit using the Farquhar biochemical photosynthesis model [57] with a Jarvis-type stomatal conductance model [61].LPJ and TRIFFID are dynamic global vegetation models that explicitly simulate vegetation types and run with only climate data as a time-variant input.The LPJ model calculates GPP using a function of CO 2 , shortwave radiation, temperature, and soil moisture using the Farquhar biochemical photosynthesis model [62] with stomatal conductance estimated from a Ball-type conductance model [59].The TRIFFID model calculates GPP as a function of CO 2 , shortwave radiation, temperature, soil moisture, and relative humidity using a leaf-level photosynthesis model [58,63] with a coupled photosynthesis and stomatal conductance model [64].
The models were used under a consistent protocol.The 1901-1930 period was used for the spin-up period dependent on each model's spin-up setting.The models were run from 1901 to 2011 with time-variable inputs of climate data (for all models) and satellite-based LAI and FPAR data (for BEAMS only).BEAMS was spun-up for 3000 years using 1901 climate data.Biome-BGC was spun-up using the 1901-1930 climate data repeatedly until the soil carbon reached equilibrium, with 6,000 years as the maximum spin-up length.LPJ was spun-up for 1,000 years using the 1901-1930 climate data repeatedly with the fire module disabled.Finally, TRIFFID was spun-up for 1000 years using the 1901-1930 climate data repeatedly.
In addition to the simulations with time-variable inputs for all climate parameters, we conducted sensitivity experiments based on an existing study to identify the causes of temporal GPP variations [65].
For each simulation, the models were run using only one climate time-series with an interannual variability in temperature, precipitation, solar radiation and atmospheric CO 2 concentration; the other climate components were held at their 1901 values.We tested four climate parameters (i.e., temperature, precipitation, solar radiation and atmospheric CO 2 concentration) that are commonly incorporated in terrestrial biosphere models.Using the NDVI data from 1982 to 2011, trends in the NDVI were estimated in each grid and each subcontinental region in a manner similar to an existing study [33].The Mann-Kendall trend test [66] was used for statistical analysis because this technique does not require an assumption of normality in variance and is more robust against anomalous outliers compared with linear regression analysis.The slopes of the trends were calculated based on the Kendall robust line-fit method [67], which is a method for robust linear regression that chooses the median slope among all lines through pairs of two-dimensional sample points.The method is a non-parametric method and is relatively insensitive to outliers.
To avoid including changes in non-vegetation effects (e.g., snow) in the calculation of temporal changes, we focused on the annual average during the vegetation growing season in each grid cell.The vegetation growing season is defined in each grid as those months in which the 1982-2011 monthly average air temperatures were greater than 5 °C and the monthly average NDVI was greater than 0.1.
Climate trends (temperature and precipitation) were also estimated from 1982 to 2011 to interpret the NDVI trends using the same approach as for the NDVI analysis, that is, the Kendall robust line-fit method [62].We used the growing season average for temperature and the annual total for precipitation in the trend estimation.

Evaluation of Trends in NDVI and Modeled GPP from 1982 to 2011
Recent trends in GPP from 1982 to 2011 were analyzed using model simulation results with time-variable inputs.Changing trends in the annual total GPP, which were estimated using the same methods as described in Section 2.4.1, were compared with satellite-based NDVI trends.Furthermore, we also evaluated the interannual variations in the NDVI and modeled GPP based on the correlation analysis between the NDVI and modeled GPP using Kendall's rank correlation coefficients (Kendall's τ) [66] to determine the ability of each model to simulate interannual variability consistently with the satellite-based NDVI.

Attribution of the Detected Changes in NDVI and GPP
Factors determining the modeled GPP trends were analyzed using the results of sensitivity experiments for each terrestrial biosphere model.A series of model simulations was conducted by changing a single time-series climate variable-temperature, precipitation, solar radiation, or atmospheric CO 2 concentration-while keeping the other climate variables fixed.To ensure that the comparison was consistent among the models, we used outputs from three terrestrial biosphere models that can only be run with climate data (i.e., Biome-BGC, LPJ and TRIFFID).The BEAMS model was not used in this analysis because the model uses time-variable satellite-based LAI and FPAR, which already contains climate variation effects, and isolating the unique effect of a specific climate parameter is difficult.
We evaluated the relative contribution of each individual factor-that is, temperature, precipitation, solar radiation, and atmospheric CO 2 concentration-to the simulated GPP trends under actual climate variations.Using the model results simulated by changing a single time-series climate variable, we evaluated the relative contribution of each individual factor by plotting them in stacked bar graphs in each region and for each model.

Observed Climate and NDVI Trend for 1982-2011
For the 1982-2011 period, the temperature increased over most regions; precipitation had a geographically variable pattern in the sign of the trends (Figure 2).Distinct increases in temperature were distributed over eastern Siberia (60°N-75°N, 110°E-180°E) and semiarid regions around inner-Mongolia (45°N-50°N, 100°E-130°E).Precipitation had distinct increases in eastern Siberia (60°N-75°N, 110°E-180°E) and part of India (10°N-25°N, 70°E-80°E), with distinct decreases in semiarid regions around inner-Mongolia (45°N-50°N, 100°E-130°E).Increasing trends in the NDVI were dominant in Asia for the 1982-2011 period (Figure 3).Moreover, 41.8% of the vegetated regions had increased NDVI values, with distinct increases in eastern Siberia (50°N-80°N, 100°E-180°E), boreal Eurasia (45°N-60°N zones), parts of China and India, and some scattered regions in tropical Asia.However, the regions with significant NDVI decreases were distributed among small regions (2.9% of the vegetated regions); for example, parts of western Arctic Siberia (60°N-80°N, 80°E-100°E) and semiarid regions around inner-Mongolia (45°N-50°N, 100°E-120°E).In the sub-continental scale analysis, we detected a statistically significant NDVI increase in seven of the nine sub-continental regions (Table 1).The two exceptions were the Boreal Eurasia North West region and the Tropical Asia South region.However, these regions had positive but minor increases in the NDVI.In the Boreal Eurasia North West region, interannual variations in the NDVI gradually increased from 1982 to 2000 and decreased after 2001, except for 2011 (Figure 4a).In the Tropical Asia South region, the NDVI had large interannual variations without significant increasing or decreasing trends from 1982 to 2000 (Figure 4i).1b and Table 1).Normalized NDVI anomalies and modeled GPPs (based on the average for the 1982-2011 period) are shown.

Evaluation of the 30-Year Trends in NDVI and Modeled GPP
The terrestrial biosphere models suggest increasing GPP in Asia from 1982 to 2011, which is mostly consistent with the satellite-based NDVI trends.The Biome-BGC and LPJ models exhibited significant increasing GPP trends in all nine regions.Moreover, the BEAMS and TRIFFID models showed significant increasing GPP trends in eight regions (Table 1).The consistency between models and the NDVI is highest for the BEAMS model because the model uses satellite-based LAI and FPAR as inputs; only the Tropical Asia South region did not exhibit a difference in the trend significance (e.g., the BEAMS model predicts statistically significant increases in GPP, whereas the NDVI does not) likely due to uncertainties in the solar radiation data used for BEAMS models.In both the Boreal Eurasia North West region and the Tropical Asia South region, the prognostic terrestrial biosphere models (Biome-BGC, LPJ, and TRIFFID) suggest significant increases in GPP.However, a similar trend is not detected for the NDVI.
The interannual variations in the NDVI and modeled GPP were highly consistent in temperate and boreal regions and weakly consistent in tropical regions (Table 2 and Figure 4).We found that the BEAMS-based GPP, in which GPP was calculated using satellite-based LAI and FPAR data, was highly consistent with the satellite-based NDVI (e.g., τ = 0.67 -0.83 for boreal and temperate regions).The results from the prognostic biosphere models also showed good agreement with the satellite-based NDVI data regarding interannual variations in GPP (e.g., averages of τ in all boreal and temperate regions are τ = 0.48, 0.50 and 0.43 for Biome-BGC, LPJ, and TRIFFID, respectively).However, we found large deviations between the models and the satellite-based NDVI in tropical forest regions (e.g., τ < 0.3 for most tropical regions).The consistency of interannual variability and 30-year trends in boreal and temperate regions between the NDVI and modeled GPP indicates that these increases are robust for the different methodologies, that is, terrestrial biosphere models and satellite-based observations.Therefore, we can use the models to perform a process-based analysis on the underlying mechanisms.As exceptions, disagreement between the model and the NDVI data was detected in several regions and specific years.
For example, in the Boreal Eurasia North West region, mainly during the 2007-2010 period (Figure 4a), the satellite-based NDVI had negative anomalies for this period, whereas the GPPs in the prognostic models (Biome-BGC, LPJ, and TRIFFID) had positive anomalies for this period.Another distinct difference was found in the Temperate Asia North West over the period 1996-2001.In this region, there was a consistent discrepancy between the three model simulations run with only climate data (i.e., BIOME-BGC, LPJ and TRIFFID) versus the NDVI data and the BEAMS model simulation that incorporates effects of satellite-based observation.In addition, the BIOME-BGC simulations had anomalously high spikes in 2007 for all temperate regions.The disagreement between these data should be analyzed in the future to identify the potential limitations of each methodology.
Disagreements between the models and data for tropical regions are evident in the current study.Possible explanations may relate to both the satellite-based observations and the terrestrial biosphere models.In principle, monitoring tropical forests using shortwave radiations from satellite-based observations is difficult because of atmospheric effects, for example, severe cloud contamination and heavy aerosol loading in tropical regions [68].Moreover, simulating the terrestrial carbon cycle using terrestrial biosphere models is difficult in tropical forests.For example, a site-level model comparison study showed that the model performance in simulating carbon fluxes was poor in tropical forests among various ecosystems in Asia [69].The difficulty in modeling GPP in tropical forest is partly due to complex processes involving water and solar radiation controls on seasonal GPP variations and the requirement for appropriate rooting depth settings [70].

Attribution of Detected Changes in NDVI and GPP
In boreal Eurasia regions (the Boreal Eurasia North West, North East, South West and South East regions), the modeled increases in GPP are explained by variations in temperature (Biome-BGC and LPJ models), precipitation (Biome-BGC, LPJ, and TRIFFID models), and atmospheric CO 2 concentration (LPJ and TRIFFID models) (Figure 5a-d).Changes in temperature and precipitation provide similar contributions to the GPP increase.The Biome-BGC model showed that increases in the atmospheric CO 2 concentration had a small effect on the GPP increase.However, the LPJ and TRIFFID models showed that changes in atmospheric CO 2 concentration, temperature, and precipitation had similar effects on GPP.
In the temperate Asia regions (the Temperate Asia North West, North East and South East regions), the causes of the GPP increases were not consistent among the models (Figure 5e-g).The GPP increases were primarily due to changes in temperature in the Biome-BGC model.According to the LPJ and TRIFFID models, atmospheric CO 2 concentration change is the main driver of the GPP increases.
In tropical Asia (the Tropical Asia North and South regions), the GPP increases were mostly explained by changes in atmospheric CO 2 concentration in the LPJ and TRIFFID models.In the Biome-BGC model, changes in precipitation and the atmospheric CO 2 concentration equally contributed to the increases in GPP.Changes in incoming solar radiation were another important factor in the GPP increases in the Tropical Asia South region in the Biome-BGC and LPJ models.However, quantifying the solar radiation contribution requires better estimates of solar radiation (see Section 3.4).
Although the estimated GPP trends were similar for the different models and were generally consistent with the satellite-based NDVI trends, the contribution of CO 2 fertilization to the increase in GPP differed among the models and regions (Figure 5).The CO 2 fertilization effect on the GPP trends was lower in the Biome-BGC model compared with the LPJ and TRIFFID models.Therefore, further analysis using data from field observations [71] is required to elucidate the effects of CO 2 fertilization on GPP.

Limitations and Potential Further Studies
Based on the analysis of the four terrestrial biosphere models and the satellite-based NDVI data, we found that the detected NDVI trends are mostly consistent with model-based GPP trends in Asia.These detected NDVI and GPP trends were analyzed to identify the important processes driving changes in the NDVI and GPP using the terrestrial biosphere models.
To fill the gaps in the results from the satellite-based NDVI and model-based GPP, additional analyses are required.First, the terrestrial biosphere models need to be further improved to reproduce observed carbon fluxes in Asia [69,72] and other regions, for example, North America [73] and Europe [74].Further refinements to these models using ground-based observation networks [75] may improve the consistency between the modeled GPP and the satellite-based NDVI and inter-model consistency.Second, an analysis of human-induced disturbances and land-use changes, including afforestation and agriculture, should be assessed in future studies because afforestation in China [22] and the increased irrigation area and fertilization use in India [28] have been found to contribute to increases in terrestrial GPP and the NDVI.In addition, the terrestrial biosphere models used in this study do not include disturbance effects or land-use changes, including plantation and deforestation in the tropical Asia.Third, to assess each model's sensitivity to climate anomalies, terrestrial GPP anomalies in response to certain anomalous climate events should be analyzed.For example, the Boreal Eurasia South East region experienced cooler and cloudier meteorological conditions in the summer of 2003, corresponding to negative anomalies in GPP.All of the models showed the same negative GPP anomalies in 2003 (Figure 4d).The mechanisms should be carefully assessed using observations and models to determine if each model reproduces GPP anomalies through a similar mechanism.Fourth, a further assessment of model-data mismatches should be conducted to clarify the potential issues for both models and satellite-based observations.In this study, we demonstrated several model and data mismatches in Section 3.2.Identifying the mismatch causes may improve the performances of terrestrial biosphere models.
A quantitative comparison of the models and data are also required in the further studies even though current qualitative assessments showed consistent variations in modeled GPPs and satellite-derived NDVI.To date, satellite-based vegetation indices have often been used for the evaluation of terrestrial biosphere models in continental-scale analysis [33].A possible method to conduct a quantitative comparison is to use GPPs of empirically upscaled products that rely on the dense network of eddy-covariance observations and satellite-based data [12].The eddy-covariance network dataset in Asia is insufficient to cover various ecosystems in Asia, especially in Siberia and Southeast and South Asia [70]; further improvements in the dataset are required.
An improved solar radiation dataset is required to interpret the effects of solar radiation on the modeled GPP.In this study, we used empirically estimated surface radiation based on the MTCLIM model with corrections using satellite-based solar radiation data during the 2000-2010 period.The method relies on the assumption that variations in diurnal temperature changes are related to incoming solar radiation, which may produce large uncertainties in estimating solar radiation.Other satellite-based radiation data (e.g., SRB data; [76]) cover the period from 1983 to 2007; no data are available for the past 30 years.Solar radiation is one of the important variables that determines GPP [17]; improved solar radiation data are thus required to improve the analysis of GPP.

Conclusion
Many regions in Asia have experienced vegetation greening from 1982 to 2011.Trends in vegetation greening have been reproduced using multiple terrestrial biosphere models for most regions in Asia.The analysis in this study clearly showed statistically significant increases in the NDVI for almost 40% of the regions in Asia, whereas decreasing vegetation trends were found for only a small percentage of the studied regions (less than 5%).Among the nine sub-continental regions, seven regions showed statistically significant increases in the NDVI.Simulations using the terrestrial biosphere models showed vegetation greening trends that were consistent with the results for the NDVI in boreal and temperate regions.
Three terrestrial biosphere models were used to determine the causes of the detected increasing trend in the NDVI.The determined contributions of climate and atmospheric CO 2 changes varied considerably among the models.Further assessment is required to determine the effect of CO 2 fertilization on GPP and assess the sensitivity of the models to climate anomalies.

Figure 2 .
Figure 2. Trend in (a) growing season average temperature (°C•30 y −1 ) and (b) annual total precipitation (%•30 y −1 ) from 1982 to 2011.The trend slope was estimated using the Kendall robust line-fit method.Solid black lines in the map show the sub-continental division (only for the Eurasian continent) defined in Figure 1b.

Figure 3 .
Figure 3. Trend in growing season average Normalized Difference Vegetation Index (NDVI) from 1982 to 2011.Only grids with statistically significant trends (p < 0.05) are colored.The Mann-Kendall trend test was applied for the significance test; the slope was estimated using the Kendall robust line-fit method.Solid black lines in the map show the subcontinental division (only for the Eurasian continent) defined in Figure 1b.

Figure 4 .
Figure 4. Interannual variations in the NDVI and annual GPP in each subcontinental region (Figure1band Table1).Normalized NDVI anomalies and modeled GPPs (based on the average for the 1982-2011 period) are shown.

Figure 5 .
Figure 5. Contribution of the modeled GPP trend of each individual factor (i.e., temperature, precipitation, solar radiation, and atmospheric CO 2 concentration) as stacked bar graphs in each region ((a) Boreal Eurasia North West, (b) Boreal Eurasia North East, (c) Boreal Eurasia South West, (d) Boreal Eurasia South East, (e) Temperate Asia North West, (f) Temperate Asia North East, (g) Temperate Asia South East, (h) Tropical Asia North, and (i) Tropical Asia South).Cumulative sums of GPP trends for the simulations with each input data temporally changing are approximately consistent with the GPP trends in the simulation with all input data temporally changing (y = 1.21 x + 0.09, R 2 = 0.92; x is the sum of GPP trends for the simulations with each input data changing in each region and for every model, and y represents the GPP trends for the simulation with all input data changing in each region and for each model; 27 data points (3 models × 9 sub-regions)).

Table 2 .
Correlation coefficients between the NDVI and GPP from each terrestrial biosphere model.Correlation coefficients were calculated based on Kendall's rank correlation coefficients (Kendall's τ).