Distinguishing the drivers of trends in land carbon fluxes and plant volatile emissions over the past 3 decades

The terrestrial biosphere has experienced dramatic changes in recent decades. Estimates of historical trends in land carbon fluxes remain uncertain because long-term observations are limited on the global scale. Here, we use the Yale Interactive terrestrial Biosphere (YIBs) model to estimate decadal trends in land carbon fluxes and emissions of biogenic volatile organic compounds (BVOCs) and to identify the key drivers for these changes during 1982–2011. Driven by hourly meteorology from WFDEI (WATCH forcing data methodology applied to ERA-Interim data), the model simulates an increasing trend of 297 Tg C a in gross primary productivity (GPP) and 185 Tg C a in the net primary productivity (NPP). CO2 fertilization is the main driver for the flux changes in forest ecosystems, while meteorology dominates the changes in grasslands and shrublands. Warming boosts summer GPP and NPP at high latitudes, while drought dampens carbon uptake in tropical regions. North of 30 N, increasing temperatures induce a substantial extension of 0.22 day a for the growing season; however, this phenological change alone does not promote regional carbon uptake and BVOC emissions. Nevertheless, increases of leaf area index at peak season accounts for ∼ 25 % of the trends in GPP and isoprene emissions at the northern lands. The net land sink shows statistically insignificant increases of only 3 Tg C a globally because of simultaneous increases in soil respiration. Global BVOC emissions are calculated using two schemes. With the photosynthesis-dependent scheme, the model predicts increases of 0.4 Tg C a in isoprene emissions, which are mainly attributed to warming trends because CO2 fertilization and inhibition effects offset each other. Using the MEGAN (Model of Emissions of Gases and Aerosols from Nature) scheme, the YIBs model simulates global reductions of 1.1 Tg C a in isoprene and 0.04 Tg C a in monoterpene emissions in response to the CO2 inhibition effects. Land use change shows limited impacts on global carbon fluxes and BVOC emissions, but there are regional contrasting impacts over Europe (afforestation) and China (deforestation).


Introduction
The terrestrial biosphere interacts with the atmosphere through photosynthesis and biogenic volatile organic compound (BVOC) emissions.Annually, terrestrial ecosystems assimilate ∼ 120 petagrams of carbon (Pg C) from the atmosphere (Beer et al., 2010), most of which reenters atmosphere through respiration and decomposition, resulting in a net global land carbon sink of 2.6 ± 0.7 Pg C a −1 (Le Quèrè et al., 2009;Sitch et al., 2015).Global BVOC emissions are estimated to be about 1 Pg C per year (Carslaw et al., 2010).These emissions are important precursors of atmospheric oxidants and aerosols, both of which affect surface air quality and exert additional regional and global chemical climate forcings (Scott et al., 2014;Unger, 2014).Observations and simulations have shown significant changes in terrestrial carbon assimilation and BVOC emissions in the past 2-3 decades (Lathiere et al., 2006;Sarmiento et al., 2010;Sindelarova et al., 2014;Sitch et al., 2015).Understanding drivers of these trends is important for the projections of future carbon fluxes, water cycle, air quality, and climatic responses.
Trends in land carbon assimilation and BVOC emissions are related to the changes in atmospheric CO 2 , meteorology, and human land use and land cover change perturbations.Elevated CO 2 promotes plant photosynthesis (Ainsworth and Long, 2005) but can directly inhibit isoprene productions (Arneth et al., 2007).Warming accelerates both carbon uptake and BVOC emissions when temperature is not above the thermal optimum (25-30 • C for photosynthesis and 35-40 • C for isoprene emission) for ecosystems that are not water stressed (Farquhar et al., 1980;Guenther et al., 1993;Piao et al., 2013).Additional warming above thermal optimum may decrease photosynthesis but still promote respiration, reducing net carbon uptake by plants (Liang et al., 2013).Increased temperatures also indirectly influence carbon exchange and BVOC emissions through the extension of growing season (Piao et al., 2007).Drought decreases gross primary productivity (GPP) and net primary productivity (NPP) (Zhao and Running, 2010) but may temporally enhance isoprene emissions (Monson et al., 2007).Land use change affects the regional carbon budget and BVOC emissions through either additional emissions or land cover changes due to deforestation, forest management, and agricultural activities (Lathiere et al., 2006;Houghton, 2010).
Estimates of recent decadal global trends in the land carbon budget and BVOC emissions are limited and uncertain due to the lack of observations.The earliest site-level measurements of land carbon fluxes were set up in the 1990s (Wofsy et al., 1993).The flux-tower data sets provide longterm records of regional carbon exchange with high precision but low spatial representation.In contrast, satellite products, such as GPP and NPP retrievals from the Moderate Resolution Imaging Spectroradiometer (MODIS) (Zhao et al., 2005) and isoprene emissions based on tropospheric formaldehyde columns from the Global Ozone Monitoring Experiment (Palmer et al., 2006), improve the spatial coverage but usually are available for only a relatively short time period (months to several years) and suffer from systematic biases when compared with ground measurements (e.g., Heinsch et al., 2006;Marais et al., 2012).Terrestrial biosphere models, evaluated with both site-level and satellitebased observations, are useful tools to estimate trends and attribute drivers of changes in land carbon fluxes and BVOC emissions (e.g., Mao et al., 2013;Stavrakou et al., 2014;Sitch et al., 2015).
In this study, we use the Yale Interactive terrestrial Biosphere (YIBs; Yue and Unger, 2015) model driven with long-term reanalysis meteorology to study the global trends of land carbon fluxes and BVOC emissions over the past 3 decades.The YIBs model is a process-based vegetation model that simulates complete land carbon cycle, including photosynthesis, plant/soil respiration, carbon allocation, and tree growth.Simulated carbon fluxes has been fully validated with carbon fluxes from 145 flux-tower sites and multiple satellite products (Yue and Unger, 2015).The YIBs model provides a unique tool to identify drivers of decadal trends in carbon fluxes and BVOC emissions because of the distinct treatments of plant phenology and BVOC emissions.
Many state-of-the-art vegetation models suffer from poor representations of phenology, which may lead to large biases in the simulated carbon fluxes (Richardson et al., 2012).The optimized phenology in YIBs is based on assessment of 13 existing models (nine for spring and four for autumn) and has been validated against both ground-based records and multiple satellite retrievals (Yue et al., 2015).In addition, YIBs incorporates two independent isoprene emission schemes within the exact same host model framework (Unger et al., 2013;Zheng et al., 2015): (i) a photosynthesisdependent isoprene emission scheme (Niinemets et al., 1999) and (ii) the Model of Emissions of Gases and Aerosols from Nature (MEGAN) isoprene scheme (Guenther et al., 2012) that is widely used in chemistry-transport modeling.Therefore, the YIBs model allows us to investigate modeling uncertainties due to the differences in the BVOC emission algorithms themselves.The major goals of this study are to identify: (1) the dominant drivers of the 30-year trends in carbon fluxes and BVOC emissions from elevated CO 2 , changes in meteorology (temperature, radiation, and soil moisture), and human land use change; (2) the feedback of biosphere, including changes in phenology and leaf area index (LAI), to the trends of land carbon uptakes and BVOC emissions; and (3) the discrepancies in BVOC trends due to application of different isoprene emission schemes.

Observations and benchmark products
We use long-term global measurements of LAI, GPP, and NPP to validate the simulated trends.The LAI data set for 1982-2011 is retrieved based on the Normalized Difference Vegetation Index (NDVI) from Global Inventory Modeling and Mapping Studies (GIMMS) with 1/12 • resolution and a 15-day interval (Zhu et al., 2013).We also use LAI data for 2000-2011 from the MODIS (http://modis.gsfc.nasa.gov/).GPP benchmark products of 1982-2011 are upscaled from the FLUXNET eddy covariance measurements using an ensemble of regression trees (Jung et al., 2009).As a comparison, we also use the GPP and NPP data sets for 2000-2011 from the MODIS, which have been developed based on remote sensing of biome parameters and assimilated meteorology (Zhao et al., 2005).All the data sets are interpolated to the monthly interval at the 1 • × 1 • offline YIBs model resolution.

Model
The YIBs model is a process-based terrestrial vegetation model that simulates the land carbon budget and dynamic tree growth (Yue and Unger, 2015).The model adapts routines from the mature TRIFFID (Cox, 2001) and CASA (Schaefer et al., 2008) models with special updates in the parameterizations of ozone vegetation dam-age (Yue and Unger, 2014), plant phenology (Yue et al., 2015), and the photosynthesis-dependent isoprene emission (Unger et al., 2013).The model simulates carbon uptake for nine plant functional types (PFTs) including tundra, C3/C4 grass, shrubland, deciduous broadleaf forest (DBF), evergreen needleleaf forest (ENF), evergreen broadleaf forest, and C3/C4 cropland.The vegetation biophysics calculates leaf-level photosynthesis using the well-established Farquhar scheme (Farquhar et al., 1980;von Caemmerer and Farquhar, 1981) and the stomatal conductance model of Ball and Berry (Collatz et al., 1991).The canopy radiative transfer scheme computes direct and diffuse photosynthetically active radiation (PAR) for sunlit and shaded regions for an adaptive number of layers.The leaf photosynthesis is then integrated over all canopy layers to generate the GPP.
Part of the assimilated carbon is used for maintenance and growth respiration, and the rest is allocated among different pools for plant development.The model calculates phenology for deciduous forests using cumulative temperature summation with additional constraints from chilling and photoperiod (Yue et al., 2015).The phenology of shrubland and grassland is jointly determined by the temperature-and drought-dependent metrics.The LAI is then updated daily based on phenology and the net carbon assimilation.The soil respiration scheme considers carbon flows among 12 biogeochemical pools, including three live pools and nine dead pools.The land carbon source or sink is calculated as the difference between the net carbon assimilation and soil respiration.
The YIBs model incorporates two independent leaflevel isoprene emission schemes embedded within the exact same host model framework (Zheng et al., 2015).The photosynthesis-based (PS_BVOC) isoprene scheme calculates emissions based on the electron transport-limited photosynthesis rate, canopy temperature, and intercellular CO 2 concentrations (Niinemets et al., 1999;Arneth et al., 2007;Unger et al., 2013).The MEGAN scheme applies commonly used leaf-level empirical functions of light and canopy temperature (Guenther et al., 1993).Both schemes implement CO 2 inhibition effects on BVOC emissions parameterized as a reciprocal empirical function of intercellular [CO 2 ] following the observations from Possell et al. (2005).For monoterpene emissions, the YIBs model applies the same temperature-dependent scheme as Lathiere et al. (2006) but with CO 2 -inhibition effects.The leaf-level BVOC emissions are integrated over the multiple canopy layers following the same approach as GPP to obtain the total canopy-level emissions.
YIBs can be used in three different configurations with increasing complexity: (1) offline local site level, which is driven with hourly measurements of CO 2 concentrations and meteorology at flux-tower sites; (2) offline global forced with spatially uniform but annually updated CO 2 concentrations and hourly gridded reanalysis meteorology; (3) online coupled to the NASA ModelE2 driven with simulated meteorology by the GCM every half-hour.At the site level, YIBs simulates reasonable seasonality (correlation coefficient R > 0.8) of GPP at 121 out of 145 flux-tower sites with biases in magnitude ranging from −19 to 7 % depending on PFTs.On the global scale, the offline model simulates an annual GPP of 125 ± 3 Pg C and net ecosystem exchange of −2.5 ± 0.7 Pg C for 1982-2011, with seasonality and spatial distribution consistent with both satellite observations and benchmark synthesis products (Yue and Unger, 2015).However, the model does not include a fully coupled carbon-nitrogen cycle, which may overestimate CO 2 fertilization effects.In addition, phenology of evergreen trees is set to constant value of 1, leading to underestimation of phenological feedbacks to flux trends.In this study, we use the (2) offline global version of the model, which is driven with global meteorology reanalysis data and observed CO 2 concentrations.

Simulations
We apply observed historical atmospheric CO 2 concentrations from the fifth assessment report (AR5) of the Intergovernmental Panel on Climate Change (IPCC) (Meinshausen et al., 2011).We apply an annually varying historical transient land cover data set (Oleson et al., 2013), which is developed based on a combination of remote sensing data from both MODIS (Hansen et al., 2003) and the Advanced Very High Resolution Radiometer (AVHRR) (Defries et al., 2000) and from land use change from Hurtt et al. (2011).We use hourly meteorological variables for 1980-2011 from the WATCH forcing data (WFD) methodology applied to ERA-Interim data (WFDEI; Weedon et al., 2014).The WFDEI reanalysis is an update of the WFD, which is developed based on the European Centre for Medium-range Weather Forecasts (ECMWF) ERA-40 reanalysis (Uppala et al., 2005).Meteorological variables applied include surface air temperature, specific humidity, wind speed, surface pressure, total PAR, and soil temperature and wetness.All of the forcing data are interpolated to the 1 • × 1 • model resolution at the hourly interval.
We perform 10 sensitivity simulations to distinguish driving factors for the changes in land carbon fluxes and BVOC emissions in the past 3 decades (Table 1).The control simulation (CO2_MET_LUC) uses interannually varying meteorology, [CO 2 ], and land cover for 1980-2011.The CO2_MET run is the same as the control simulation but prescribes land cover at the year 1980.Three single-factor runs prescribe most boundary conditions at the year 1980 but allow the interannual variations of [CO 2 ] (CO2_ONLY), land cover (LUC_ONLY), and meteorology (MET_ONLY).Results from these runs are compared with that of control simulation to determine the dominant drivers of simulated trends.To understand the impact of individual meteorological variables, three additional runs are performed with fixed (or recycled) [CO 2 ], land cover, and all meteorol- ogy for the year 1980 but one field varying for 1980-2011 each time, including temperature (TEMP_ONLY), PAR (PAR_ONLY), and soil wetness (SOILW_ONLY).Finally, two runs are performed to examine feedback of biospheric changes.LAI_ONLY prescribes all boundary conditions at the starting year 1980 but implements the year-to-year LAI simulated by the control run.PHEN_ONLY also prescribes all forcings at the starting year except for the year-to-year phenology from control simulation.All simulations are initialized following the same spin up process (Yue and Unger, 2015) and are integrated for 1980-2011.

Drivers of trends in LAI
Observations show an increasing trend of LAI on most of vegetated continents, especially in Europe, northern and eastern Asia, central Africa, and southeastern USA in the past 3 decades (Fig. 1a).The simulation with year-to-year [CO 2 ], land cover, and meteorology reproduces the magnitude of trend in Europe and the sign of trend in northern Asia, eastern USA, central Asia, and Australia (Fig. 1b).The model predicts negative changes in central Africa, western USA, eastern Asia, and eastern South America, which are inconsistent with satellite observations.These negative trends are mainly contributed by the changes in meteorology (Fig. 1e), except for that in East Asia where land cover changes due to human activities result in the decline of LAI (Fig. 1f).Without the land use perturbation, the negative LAI trend in East Asia is weakened and the prediction is closer to observations (Fig. 1c).For the individual drivers, CO 2 fertilization leads to widespread increases in LAI (Fig. 1d), meteorology causes dipole changes on most continents (Fig. 1e), and land use change generally results in negative trends (Fig. 1f).Regionally, simulation CO2_MET_LUC shows a positive trend of 0.0035 m 2 m −2 a −1 in Europe (Table 2), close to the observed value of 0.0049 m 2 m −2 a −1 (Fig. 1a).In other areas, simulated LAI trends are either underestimated (by 87 % in Amazon, 78 % in North America, and 48 % in Central Africa) or opposite in sign (East Asia and Indonesia) compared to observations.Such inconsistencies indicate the limit of model simulations but may also in part result from the uncertainties in the satellite measurements (see Sect. 4.1).

Drivers of trends in land carbon fluxes
Predicted GPP and NPP trends show similar spatial pattern as that of LAI (Fig. 2a and c).However, regional trends are all positive in the main continents and on the global scale (Tables 2 and 3).Tropical areas are experiencing maximum changes, especially in Central Africa (GPP by 83.3 Tg  2. C a −2 and NPP by 51.7 Tg C a −2 ) and the Amazon (52.7 and 27.1 Tg C a −2 ).In the Northern Hemisphere (NH), changes are significant in Europe (53.4 and 33.2 Tg C a −2 ), East Asia (42.4 and 27.2 Tg C a −2 ), and North America (13.6 and 9.7 Tg C a −2 ).Thirty-year historical observations of GPP and NPP are not available.Therefore, we compare YIBs predictions with MODIS land carbon fluxes over the more recent period of 2000-2011 (Fig. 3).Different from the 30year trend, land carbon fluxes over the recent decade show negative trends in southeastern USA, southern Africa, eastern Australia, and central and northern Asia (Fig. 3a and c).
Most of these changes are consistent with the MODIS observations (except for the USA, Fig. 3b and d) and are attributed to the drought tendency in the past decade (Zhao and Running, 2010).
For the 30-year trend, both CO 2 and meteorology are playing important roles (Fig. 2b and d).CO 2 fertilization dominates the GPP and NPP trends of tropical forests in the Amazon, central Africa, and Indonesia, and ENF and DBF in boreal North America, eastern Europe, and central and northern Asia.Land use change plays a limited role in land carbon cycle flux trends over the past 3 decades, except for some areas in northern Africa.Meteorological forcing drives changes in land carbon fluxes for tundra in subarctic regions, C3 grasslands in the central USA and southern Africa, C4 grasslands in central Africa and eastern South America, and shrublands in Australia and southern Asia.Soil wetness plays the dom- inant role in the tropical and subtropical areas (Fig. 4b).
The drought tendency in the western USA, central Africa, and the east of South America (Fig. S1d) results in the regional decline of land carbon fluxes (Fig. 4a).In contrast, the increasing wetness in the northern Amazon and southern Africa leads to the enhancement of regional GPP.Warming is the main cause for the GPP trends over the subarctic areas (Fig. 4b).Contribution of PAR is limited except for some areas in eastern Europe.
The simulated net ecosystem productivity (NEP) shows weaker trends compared with GPP and NPP (Fig. 2e), because NEP is offset by the significant trends in heterotrophic respiration (Rh) (Table 2).Regionally, the YIBs model predicts enhanced net land carbon uptake in boreal North America, northern Asia, and southern Africa but reduced NEP in the central USA, the Amazon, central Africa, eastern Europe, and East Asia.The simulated global NEP trends (Fig. 5d) are in broad agreement with the comprehensive bottom-up estimates by Pan et al. (2011), who found slightly decreasing net carbon uptake by global established forests (without human perturbations in the tropics but with afforestation in subtropical areas) in 2000-2007 relative to that in 1990-1999.Attribution analysis shows that the NEP trends are mainly driven by the changes in meteorological forcings (Fig. 2f),  because CO 2 fertilization enhances both NPP and Rh with similar magnitude (Fig. 5).
On the global scale, GPP, NPP, and Rh increase, respectively, by 298, 185, and 181 Tg C a −2 in the past 3 decades (Table 3).The long-term trends of carbon fluxes are mainly driven by CO 2 fertilization, while the interannual variability is related to meteorological forcings (Fig. 5).Warming alone decreases GPP especially in tropical forests (not shown) but increases autotrophic respiration (Ra), leading to global reductions of 56 in NPP and 10 Tg C a −2 in NEP (Table 3).Drought alone strongly decreases GPP, especially for tropical grassland and shrubland (Fig.   carbon fluxes, though it induces relatively large reductions in NEP (Table 3).

Drivers of trends in BVOC emissions
Simulated isoprene emission trends are sensitive to the choice of modeling scheme.With the PS_BVOC scheme, global isoprene emissions increase by 0.4 Tg C a −2 during 1982-2011.Large enhancements are predicted in central Africa (0.25 Tg C a −2 ) and Europe (0.16 Tg C a −2 ), while moderate reductions are found in the western USA, eastern South America, and East Asia (Fig. 6a).Drought accounts for the decline of isoprene emissions in the USA and South America, but land use change is the main driver for the reductions in East Asia (Fig. 6b).Increasing [CO 2 ] promotes photosynthesis but meanwhile inhibits BVOC emissions, leading to offsetting CO 2 effects on isoprene.Consequently, the global isoprene emission is mainly driven by meteorological changes (Fig. 6b).In contrast, using MEGAN scheme, the YIBs model simulates a global reduction of 1.1 Tg C a −2 for isoprene emissions (Fig. 6c).Strong declines are found in the tropical rainforest, for example in the Amazon (−0.43 Tg C a −2 ), central Africa (−0.14 Tg C a −2 ), and Indonesia (−0.16 Tg C a −2 ) (Fig. 6c).The MEGAN scheme is sensitive to both light and temperature (Guenther et al., 1995).The strong positive brightening trends in PAR in Europe (Fig. S1b in Supplement) promote isoprene emissions there.The positive impacts of NH warming (Fig. S1a) are compensated by CO 2 inhibition, leading to small changes in isoprene emissions (Fig. 6c).In the tropical areas, where trends of temperature and PAR are limited, CO 2 inhibition results in strong reductions of BVOC emissions.Monoterpene emissions show a global reduction of 0.04 Tg C a −2 over the past 3 decades (Fig. 6e).

Feedback of biospheric changes to the trends
Due to the changing climate and CO 2 fertilization, the biosphere is experiencing significant changes in the past 3 decades.The most evident alterations include LAI changes in peak season and phenological changes in growing and falling seasons.In this section, we explore the feedback of these biospheric changes to the carbon uptake and BVOC emissions.

Impacts of LAI changes
Sensitivity run LAI_ONLY retains the trends in LAI but prescribes other forcings.In this simulation, trends in GPP (Fig. S2a) and NPP (Fig. S2c) generally follow that in LAI (Fig. 1b), but with smaller magnitude relative to those in control simulations (Fig. 2a and c).LAI in the north of 30 • N shows widespread increases in both observations and simulations (Fig. 1a and b).Over these northern lands, the unit change in leaf area leads to enhancement of regional GPP by 32 Pg C a −1 , much lower than the response of 116 Pg C a −1 LAI −1 for the simulation including CO 2 fertilization and climate forcings (Fig. 7a).In the tropical areas, both positive and negative LAI trends are predicted due to the competition between CO 2 fertilization and drought effects (Fig. 1).As a result, LAI-induced GPP and NPP changes show patchy  distributions at tropics (Fig. S2a and c), leading to moderate changes in the global carbon assimilations (Table 3).Trends in isoprene emission (calculated with the PS_BVOC scheme) also follow that of LAI, except that leaf expansion results in decreased emissions at high latitudes (∼ 60 • N, Fig. S2e).The cause for such inconsistency is un-clear but might be because the denser leaves reduce radiation penetrating to lower canopy layers.Such an impact would only affect BVOC emissions at high latitudes because PAR is usually limiting near subarctic areas.In most of the subtropical areas, increased LAI leads to enhanced isoprene emissions.On average, unit change in LAI at north of 30 • N leads to enhanced isoprene emissions by 43 Tg C a −2 , only 25 % of the magnitude in simulation CO2_MET (Fig. 7b).A similar ratio of 23 % is achieved for MEGAN isoprene emissions.These results are consistent with that for GPP (Fig. 7a), suggesting that CO 2 fertilization and meteorological changes are the main drivers for the changes in carbon uptake and BVOC emissions, even over the northern lands where the most evident changes in LAI are observed.

Impacts of phenological changes
Plant phenology, which is the timing of budburst and leaf fall, is closely related to temperature, moisture, and photoperiod and thus is experiencing significant changes in the past decades following climate change (Jeong et al., 2011;Keenan et al., 2014;Buitenwerf et al., 2015;Yue et al., 2015).Extension of the growing season has the potential to promote carbon uptake of forests (e.g., Piao et al., 2007;Richardson et al., 2009).However, such inference requires careful interpretation because the phenological changes are usually accompanied with warming and elevated [CO 2 ], both of which are also contributing to the enhancement of carbon fluxes.Phenological changes are also expected to affect BVOC emissions; however, such investigations are still missing (Richardson et al., 2013).With the YIBs model, we evaluate the impacts of the growing season extension on both carbon uptake and BVOC emissions by isolating long-term phenological trends from changes in temperature and [CO 2 ].
The YIBs model simulates advanced spring and delayed autumn over most areas in NH (Fig. S3).Budburst dates advance on average by 0.16 days a −1 in Europe and 0.15 days a −1 in East Asia (Table 2) but with moderate changes or even delays in northwestern Asia and eastern Siberia (Fig. S3a).Spring is earlier by 0.14 days a −1 in eastern USA while delayed by 0.15 days a −1 in northwestern USA and southeastern Canada, leading to a minor advance of 0.01 days a −1 over North America.Dormancy onset dates are largely delayed in eastern Europe and northwestern Asia (∼ 0.3 day a −1 ), western USA (∼ 0.1 day a −1 ), boreal Canada (∼ 0.1 day a −1 ), and northeastern China (∼ 0.1 day a −1 ) (Fig. S3b).Advanced autumn (∼ 0.1 day a −1 ) is predicted in northern Asia.Most of these changes are consistent with observations from remote sensing data (Jeong et al., 2011) except for some discrepancies in the magnitude.The predicted phenological trends mainly follow the long-term changes of surface air temperature, especially that in April (for spring) and September (for autumn) (Fig. S4).Sensitivity tests without chilling requirement and photoperiod limit show similar changes (Yue et al., 2015), suggesting that temperature changes dominantly drive the trends of forest phenology in the past 3 decades.
On average, the YIBs model simulates advanced budburst by 0.12 day a −1 and delayed dormancy onset by 0.09 day a −1 at north of 30 • N in the past 3 decades (Fig. 8a and b).Observations based on remote sensing greenness show trends of −0.11 day a −1 for onset and 0.25 day a −1 for offset during 1990-2009 (Zhu et al., 2013).An ensemble prediction based on nine terrestrial models yields an advance of 0.08 ± 0.13 day a −1 for onset and a delay of 0.22 ± 0.1 day a −1 for offset (Sitch et al., 2015).Our predictions are in broad agreement with these estimates though the autumn delay is less, likely because the positive trend of offset is weaker for the recent decade (Jeong et al., 2011).
We plot the annual total GPP and isoprene emissions at north of 30 • N against the length of growing season for 1982-2011 (Fig. 8c and d).In the CO2_MET run, the 1-day extension is correspondent to increases of 0.17 Pg C a −1 in GPP and 0.34 Tg C a −1 in isoprene emissions.If only temperature is allowed to vary, the phenological trend remains the same while the increases of GPP and isoprene emissions are largely weakened.In the TEMP_ONLY run, the 1-day extension in growing season is accompanied by increases of 0.05 Pg C a −1 in GPP and 0.25 Tg C a −1 in isoprene emissions.The changes in BVOC emissions are not as dramatic as those of GPP because CO 2 has both enhancing and suppressing impacts on the former.If we further exclude temperature effects (PHEN_ONLY run), GPP increases only by 0.01 Pg C a −1 while isoprene emissions decrease by 0.1 Tg C a −1 , both of which are not statistically significant, suggesting that the phenological change alone does not promote either GPP or isoprene emissions.There are two reasons for this apparent contradiction.First, the extension of the growing season occurs in shoulder months, usually in May and September, when both GPP and BVOC emissions and their changes are much smaller compared to that in peak months (Fig. S5).Second, phenological changes are not uniform in space.As Fig. S3 shows, both positive and negative changes are predicted for budburst and dormancy onset dates.Such spatial inhomogeneity, in combination with the discrepancies in regional vegetation types and meteorological conditions, result in varied responses in GPP (Fig. S2b) and isoprene emissions (Fig. S2f).
Plant phenology at lower latitudes (30 • S-30 • N) is also experiencing dramatic changes, though such changes are diverse in phase, magnitude, or both (Buitenwerf et al., 2015).In the model, tropical phenology is mainly driven by soil wetness and as a result exhibits large changes in the past 3 decades (not shown).These changes lead to a reduction of 42 Tg C a −1 in GPP at the tropics (Fig. S2b), which accounts for 14 % of global GPP trend but with the opposite sign (Table 3), suggesting additional inhibition of drought on carbon cycle.A similar conclusion applies for BVOC emissions (Fig. S2f), though experiments suggest that isoprene production has some tolerance to mild drought conditions (e.g., Pegoraro et al., 2006).However, changes in droughtdependent phenology are very uncertain and observations are not available for evaluation.We assume that phenological changes may have larger impacts on both carbon assimilation and BVOC emissions at tropical areas than that at higher latitudes.

Uncertainties in observations
Terrestrial biosphere modeling is a useful tool to identify drivers of long-term changes in land carbon fluxes.The reliability of simulations is dependent on the availability of observations for model validation.In this study, we use 30year LAI observations from the LAI3g product (Zhu et al., 2013) and 12-year GPP from MODIS (Zhao et al., 2005), both of which are remote sensing retrievals, to validate the simulated trends (Figs. 1 and 3).We found the offline global model biases against both fields, especially for LAI (Fig. 1).Such discrepancies may in part result from the uncertainties in measurements themselves.As a check, we compare the derived LAI trends from LAI3g with retrievals from MODIS for the overlap period of 2000-2011 (Fig. S6a and  b).Global LAI significantly increases in LAI3g but show widespread reductions in MODIS, especially over subtropical areas.Simulated trends (CO2_LUC_MET) are closer to the estimates with MODIS, especially for the changes in the NH (not shown).Meanwhile, we compare the derived GPP trends from MODIS with that upscaled from FLUXNET data using an ensemble of regression trees (Jung et al., 2009) for 2000-2011 (Fig. S6c and d).The two products show similar trends over most areas except for some discrepancies in tropical areas and the eastern USA.Simulated GPP trends match results from Jung et al. (2009) better than those from MODIS (Fig. 3a).However, we do not use Jung et al. (2009) to validate simulations for 1982-2011 because the earliest fluxtower observations began only in middle 1990s.The large discrepancies in the observed trends among different data sets not only indicate the importance of model evaluations with multiple products but also put forward the necessity of data inter-comparisons and algorithm improvements to alleviate uncertainties in observations.

Comparisons with other modeling studies
We assess the YIBs simulations within the context of other models and/or previous multi-model studies to evaluate the robustness of the predicted trends in land carbon fluxes and BVOC emissions.The YIBs model predicts NPP trends of 67.4 Tg C a −2 in northern land (25-90 • N) and 98.1 Tg C a −2 in tropical land (15 • S-25 • N), similar to the ensemble estimates of 63 ± 22 and 102 ± 34 Tg C a −2 for 1990-2009 based on nine terrestrial biosphere models (Sitch et al., 2015).However, the simulated NPP trend is only 19.8 Tg C a −2 in southern land (15-90 • S), much lower than the ensemble mean value of 53 ± 31 Tg C a −2 in Sitch et al. (2015).As for the NEP, the YIBs predicts trends of 2.0 Tg C a −2 in northern land, 1.0 Tg C a −2 in tropical land, and −0.3 Tg C a −2 in southern land, much smaller in magnitude compared with the −2.0 ± 12, 36.0 ± 13, and 21 ± 17 Tg C a −2 estimated by Sitch et al. (2015).However, their predictions are insignificant (p > 0.05) for nine, five, and seven out of nine models in the northern, tropical, and southern land, respectively, suggesting that the strengthening uptake by terrestrial ecosystem is not robust.
For the BVOC, Stavrakou et al. (2014) investigated isoprene emissions over Asia during 1979-2012 using the MEGAN scheme and taking into account both climate and land use changes.Their results showed widespread increases in the emissions over China but moderate decreases in Indonesia.In contrast, the YIBs model with the MEGAN scheme simulates widespread reductions in the same areas for 1980-2011 (Fig. 6c).The discrepancies between studies are accounted for by differences in the drivers including land cover change, meteorology, and CO 2 inhibition effects.The YIBs model is driven with land cover data from Hurtt et al. (2011), which estimates an increase of crop (non-isoprene emitter) fraction in eastern China by 0.32 % per year in the last 3 decades, at the cost of the coverage loss by 0.12 % a −1 for DBF and 0.14 % a −1 for ENF (strong BVOC emitters).However, the data from Ramankutty and Foley (1999), used by Stavrakou et al. (2014) with updates to 2007, show a reduction of the crop fraction over eastern China for a similar period.In addition, the ERA-Interim PAR used in Stavrakou et al. (2014) shows an increasing trend in southeast China (cf.their Fig.5c).On the contrary, the WFDEI PAR for YIBs exhibits a declining trend in the same region (Fig. S1b), leading to a reduction in isoprene emissions.The WFDEI surface solar radiation is based on the ERA-Interim radiation but is adjusted using the average cloud cover from the Climatic Research Unit (CRU) and takes into account the effects of interannual changes in atmospheric aerosols (Weedon et al., 2011).Finally, the YIBs simulations include CO 2 inhibition effects on BVOC emissions, which were neglected in Stavrakou et al. (2014).Naik et al. (2004) predicted a global trend of 1.3 Tg C a −2 for isoprene emissions during 1971-1990 using the Integrated Biospheric Simulator (IBIS) driven with monthly mean CRU meteorology.Lathiere et al. (2006) Müller et al. (2008) reported a global increase of 4.5 Tg C a −2 for 1995-2006 using a canopy environmental model and the NCEP meteorological data.In contrast to these previous studies, YIBs with the MEGAN scheme simulates a decreasing trend of ∼ 1 Tg C a −2 in the past 3 decades.The main cause of the discrepancy in the sign of change is the missing CO 2 inhibition effects in the previous studies.In addition, differences in vegetation models, meteorological forcings, and time frames of investigation also likely contribute.The YIBs result is consistent with a recent study by Sindelarova et al. (2014), who reported a decreasing trend of

Impacts of CO 2 effects
Similar to the multi-model ensemble predictions (Sitch et al., 2015), we found that global trends in carbon fluxes are dominantly driven by CO 2 fertilization (Figs. 2 and  5).In the YIBs, the global responses to elevated [CO 2 ] is 0.2 % ppm −1 for GPP and 0.27 % ppm −1 for NPP, with relatively uniform spatial distribution (Figs.S7a and b).The GPP response falls within the range of 0.05-0.21% ppm −1 predicted by 10 terrestrial models (Piao et al., 2013) and that of 0.01-0.32% ppm −1 observed from multiple free-air CO 2 enrichment (FACE) sites (Ainsworth and Long, 2005).The NPP response is higher than the model ensemble of 0.16 % ppm −1 (Piao et al., 2013) and the observed median value of 0.13 % ppm −1 (Norby et al., 2005), suggesting that CO 2 fertilization to NPP may be overestimated in the YIBs.One possible cause is the omission of N limitation in the model, which could reduce CO 2 responses by half (Piao et al., 2013).Elevated [CO 2 ] leads to increases of 0.023 Pg C a −1 ppm −1 in NEP, within the multi-model range of 0.003-0.06Pg C a −1 ppm −1 (Piao et al., 2013).
Responses of BVOC emissions to elevated [CO 2 ] are different between PS_BVOC and MEGAN schemes (Fig. S7c  and d).PS_BVOC includes both CO 2 fertilization (on photosynthesis) and inhibition (on isoprene) effects, leading to moderate but generally positive changes in isoprene emissions.In contrast, emissions from the MEGAN scheme are not dependent on foliar photosynthesis and as a result only CO 2 inhibition is enforced.Chamber experiments show contrary tendencies for photosynthesis and isoprene in response to elevated [CO 2 ] (Possell et al., 2005), supporting the simulations with MEGAN.In addition, the magnitude of CO 2 inhibition implemented in MEGAN (−0.25 % ppm −1 ) is close to observations (−0.26 % ppm −1 ) in Possell et al. (2005).However, most of these experiments are conducted for a short-term period and cannot detect LAI changes due to the long-term CO 2 fertilization.In addition, the impacts of CO 2 are dependent on species and environmental conditions (ambient temperature and light availability).For example, Buckley (2001) found almost no responses in isoprene emissions to the elevated [CO 2 ] for oak trees.Furthermore, experiments with high temperature and/or light density show increasing isoprene at elevated [CO 2 ] (Sun et al., 2013).These studies suggest that the real responses of isoprene emissions to CO 2 under long-term climate change may not be so linear as predicted in MEGAN scheme.More sensitivity experiments and long-term samplings are required to identify CO 2 -isoprene relationships on broad range of biomes and locations.

Impacts of meteorology
Predicted long-term trends show large deviations against observations at tropical areas (Fig. 3), where meteorology plays important and complex roles.Responses of carbon fluxes to temperature are more diverse than to CO 2 (Fig. S8a and b).In the YIBs, negative responses of GPP and NPP are predicted in tropical areas, where soil moisture availability limits plant functions (e.g., stomatal conductance) to the increased temperature.Furthermore, for tropical rainforests where ambient temperature is higher than optimal photosynthetic temperature (25-30 • C), additional warming decreases carbon assimilation, especially for NPP because of simultaneous increases in plant respiration (Liang et al., 2013).On the contrary, warming leads to enhanced GPP and NPP at wetter and cooler areas in the NH subtropics.Such spatial pattern is consistent with multi-model ensemble predictions (Piao et al., 2013).On the global scale, warming results in changes of −0.7 % • C −1 for GPP in YIBs, falling within the range of −1.6-1.4 % • C −1 estimated by 10 models (Piao et al., 2013).Predicted NPP responses of −15-6 % • C −1 (Fig. S8b) is not so positive as the measurements of −8-40 % • C −1 , probably because most current warming experiments are located in the subtropics of NH (Wu et al., 2011) Piao et al., 2013).Simulated isoprene emissions with PS_BVOC show similar warming responses as that of carbon fluxes (Fig. S8c), except for tropical rainforests where the former is positive while the latter is negative.Such decoupling is attributed to the differences in optimal temperatures between isoprene (35-40 • C) and photosynthesis (25-30 • C).Simulations with MEGAN scheme show very strong temperature dependence of 6-15 % • C −1 (Fig. S8d), consistent with measurements of 5-20 % • C −1 for aspen (Niinemets and Sun, 2015) and 9-12 % • C −1 for oak (Li et al., 2011).However, experiments with some other species (e.g., spruce in Kivimaenpaa et al. (2013)) show no responses or moderate ones, suggesting that warming sensitivity of isoprene emissions might be dependent on species and ambient conditions.
Responses to PAR are mostly positive and distributed evenly, with global sensitivity of 0.3 % W −1 m 2 for GPP and 0.5 % W −1 m 2 for NPP (Fig. S9a and b).Isoprene emissions from both PS_BVOC and MEGAN schemes show similar responses to PAR, with larger sensitivity in the subtropics than that in the tropics (Figs.S9c and S9d), likely because the ambient PAR is higher at lower latitude, leading to slower responses of isoprene emissions to the unit changes of light (Guenther et al., 1993).YIBs simulations show that PAR is not the driver of long-term trends in carbon fluxes and BVOC emissions (Fig. 4), likely because changes in solar radiation is limited in the past 3 decades (Fig. S1b).
Soil moisture dominates climate-driven flux changes in tropical areas (Fig. 4).In YIBs model, changes in soil water of leaf stomatal conductance and plant phenology (especially for shrublands and grasslands in arid regions).Both GPP and NPP show strong responses to soil wetness variations, especially over tropics where > 10 % changes are found for every increase of 0.01 in soil wetness at 1.5 m (Fig. S10a and  b).On the global scale, GPP changes by 4.7 % 0.01 −1 and NPP by 5.5 % 0.01 −1 in response to soil wetness.Although experiments also show rapid reductions in carbon assimilation due to drought stress (e.g., Ruehr et al., 2012;Xia et al., 2014), the magnitude of such influence is difficult to evaluate because different metrics and depths of soil water are used in measurements.Isoprene emissions from PS_BVOC show similar soil-wetness responses to that of GPP (Fig. S10c), indicating that drought reduces BVOC emissions.However, observations show insignificant changes of isoprene with mild drought stress (e.g., Pegoraro et al., 2006), though such drought tolerance is strongly weakened at severe drought and/or warm conditions (Centritto et al., 2011).Consistent with these experiments, the MEGAN scheme does not include drought inhibition on isoprene emissions.Simulations with YIBs show large responses of BVOC to soil wetness in tropical areas (Fig. S10d) mainly because of the changes in drought-dependent phenology.

Impacts of land use change
Changes of land use show moderate impacts on global carbon budget (Fig. 2) and BVOC emissions (Fig. 6) in the past 3 decades, though regional perturbations are found in China and Europe.The afforestation in Europe helps promote regional carbon uptake, resulting in more reasonable trends in LAI compared with remote sensing data (Fig. 1).However, the expansion of crop in China leads to a reduction in LAI, which is not supported by the satellite data.One possible cause is the uncertainty in crop fraction, because data from Hurtt et al. (2011), used by YIBs, show crop expansion while data from Ramankutty and Foley (1999) suggest reductions of the crop fraction over eastern China over the similar period.The role of land use change in our simulation might be conservative because we consider only land cover changes.Perturbed emissions from land use management, such as forest lodging, cropping practice, use of fertilizer, fire management and so on (Houghton, 2010) may alter regional carbon budget by changing carbon sinks to sources.Studies including gross emissions of land use perturbation estimated a global net land source to atmosphere, which shows decreasing trend in the last 3 decades (Ciais et al., 2013).Such change may help strengthen net land carbon sink but is missing in our study.

Impacts of biospheric changes
The land biosphere has experienced significant changes in the past 3 decades.At north of 30 • N, changes in LAI account for 25 % of the trends in regional carbon fluxes and isoprene emissions.However, the extension of growing season alone makes insignificant contributions to the increased carbon assimilation.This conclusion is inconsistent with site-level observations that show evident increases in carbon assimilation at early spring and/or late autumn in recent decades (Dragoni et al., 2011;Keenan et al., 2014).There are two causes for such discrepancies.First, phenology at specific location may exhibit much more intense changes than that at larger scale.For example, Dragoni et al. (2011) estimated extensions of growing season by 2.3-3.3 day a −1 in Morgan-Monroe State Forest in south-central Indiana, USA, for 1998USA, for -2008.The magnitude of this change is ∼ 10 times larger than the observed value of 0.36 day a −1 from satellite and simulated value of 0.22 day a −1 with YIBs for the northern lands.Second, enhanced temperature also contributes to the stronger uptake at early spring and late autumn.One difficulty for the observation-based estimate of phenological impacts is that extension of growing season is accompanied by warmer climate, which may stimulate both carbon assimilation and BVOC production.In a recent study, Barlow et al. (2015) found invariant length of land carbon uptake period at high northern latitudes based on the first order of time differential of atmospheric CO 2 concentrations, suggesting that increased greenness is not necessarily equal to enhanced carbon uptake in shoulder seasons.Furthermore, Barlow et al. (2015) showed that enhanced peak uptake is the main driver for the strengthened carbon sink at high northern latitudes over the past 4 decades.These conclusions are supportive of our simulations for the monthly trends at subtropical regions (North America, Europe, and East Asia) (Fig. S5).

Conclusions
With YIBs model, we estimated global increases of carbon assimilation especially at tropical areas for 1982-2011.This trend is mainly attributed to the widespread CO 2 fertilization effect and jointly affected by changes in meteorology and land cover.Increase of temperature promotes carbon uptake of forest ecosystems at high latitudes (> 30 • N) while drought tendency dampens GPP and NPP of grasslands and shrublands at low latitudes (30 • S-30 • N).The widespread increases of LAI at northern lands account for ∼ 25 % of the regional GPP trends.Significant changes in phenology are found at north of 30 • N; however, this temperature-driven phenological change alone is not promoting regional carbon assimilation.Changes in land use show limited influences on global carbon fluxes, except for some regional impacts over Europe (afforestation) and China (deforestation).Due to the simultaneous enhancement in soil respiration, land carbon sink has remained almost stable in the past 3 decades.The YIBs model does not yet include a fully coupled carbonnitrogen cycle; thus the model may overestimate CO 2 fertilization effects.On the contrary, implementation of droughtdependent phenology may amplify drought inhibition effects on photosynthesis and result in an underestimation of carbon uptake.
We estimated global trends of BVOC emissions with two schemes.Simulations with PS_BVOC scheme show increasing isoprene emissions, mainly attributed to the increases of temperature.For this scheme, CO 2 effects are neutralized due to both fertilization (on photosynthesis) and inhibition (on isoprene).Simulations with MEGAN scheme show decreasing emissions of isoprene and monoterpene because of CO 2 inhibition, especially in the tropics.In subtropical areas, both schemes predict regional increases of BVOC emissions in Europe following the warming trend and afforestation but reductions in the USA and China due to cropland expansion.
The Supplement related to this article is available online at doi:10.5194/acp-15-11931-2015-supplement.

Figure 1 .
Figure 1.Comparison of trends in (b-f) simulated leaf area index (LAI) with (a) observations for 1982-2011.Observations are derived from GIMMS NDVI.Simulations are performed with either (d, e, f) single forcings or (b, c) the combinations of these forcings.Forcings considered include meteorology from WFDEI reanalysis (MET), CO 2 fertilization (CO 2 ), and land use change (LUC).For every forcing included in the simulation, the year-to-year fields are utilized.Otherwise, the forcing is prescribed at the year 1980.Only significant trends (p < 0.05) are presented.The six box regions in (a) indicate areas for statistical analyses in Table2.

Figure 2 .
Figure 2. Simulated trends in (a) gross primary productivity (GPP), (c) net primary productivity (NPP), and (e) net ecosystem productivity (NEP), and (b, d, f) the dominant drivers for these changes during 1982-2011.Simulations are performed with WFDEI reanalysis.Three factors, meteorological forcing, CO 2 fertilization, and land use change, are considered as the potential drivers of flux trends.For each grid in figures (b, d, f), the factor generating the largest (either maximum or minimum) trend with the same sign as the net change (a, c, e) is selected as the driving factor.Only significant trends (p < 0.05) are presented.

Figure 3 .
Figure 3. Comparisons of trends in (a, b) GPP and (c, d) NPP for 2000-2011 between (a, c) simulations and (b, d) observations.Observed fluxes are retrieved from the Moderate Resolution Imaging Spectroradiometer (MODIS).
4), leading to reductions of 51 in NPP and 13 Tg C a −2 in NEP.Trends in PAR do not affect GPP and NPP but may decrease NEP by 23 Tg C a −2 because soil respiration is slowly increasing to reach the equilibrium.Land use change has very limited impacts on the trends of www.atmos-chem-phys.net/15/11931/2015/

Figure 4 .
Figure 4. Simulated (a) trends in GPP driven alone with WFDEI reanalysis and the (b) drivers for such changes.Simulation in (a) is performed with year-to-year meteorological forcings but prescribed [CO 2 ] and land use in the year 1980.Simulations in (b) are the same as (a) except that the year-to-year variations are allowed only for a single meteorological variable (temperature, PAR, or soil wetness) each time.For each grid, the meteorological variable generating the largest (either maximum or minimum) trend with the same sign as the net change (a) is selected as the driving factor.Only significant trends (p < 0.05) are presented.

Figure 6 .
Figure 6.Simulated trends of (a, c) isoprene and (e) monoterpene and (b, d, f) the dominant drivers for these changes during 1982-2011.Simulations are performed with WFDEI reanalysis.Isoprene emissions are simulated with (a) PS_BVOC and (c) MEGAN schemes.Three factors, meteorological forcing, CO 2 effects (both fertilization and inhibition), and land use change, are considered as the potential drivers of flux trends.For each grid in figures (b, d, f), the factor generating the largest (either maximum or minimum) trend with the same sign as the net change (a-c) is selected as the driving factor.Only significant trends (p < 0.05) are presented.

Figure 7 .Figure 8 .
Figure 7. Responses of (a) GPP and (b) isoprene emissions to the changes in the annual average LAI at the north of 30 • N for simulations CO2_MET (red) and LAI_ONLY (blue).Both GPP and isoprene emissions are the sum of all PFTs.Isoprene is simulated with the PS_BVOC scheme.Units of trends are (a) Pg C a −1 LAI −1 and (b) Tg C a −1 LAI −1 .The spatial distribution of GPP and isoprene changes is shown in Fig. S2.
estimated an increasing global trend of 0.3 Tg C a −2 for 1983-1995 using the ORCHIDEE (Organizing Carbon and Hydrology in Dynamic EcosystEms) vegetation model driven by sub-daily variables from the NCEP/DOE (National Center for Environmental Predictions/Department of Energy) Reanalysis 2.

∼
1.2 Tg C a −2 for global isoprene emissions during 1980-MEGAN scheme and inclusion of a CO 2 inhibition parameterization fromHeald et al. (2009).

Table 1 .
Summary of model simulations driven with WFDEI reanalysis.
2 ] and land cover are prescribed at the year 1980.Phenology varies for 1980-2011.

Table 2 .
Summary of trends in different domains from the simulation CO2_MET_LUC, which is driven with WFDEI meteorology.Significant trends (p < 0.05) are shown in bold.
* Phenology is set to constant for tropical rainforest in the model.

Table 3 .
Summary of simulated trends of global carbon fluxes (Tg C a −2 ) from different experiments.Simulations are using WFDEI meteorology.Significant trends (p < 0.05) are shown in bold.