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

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


17
The terrestrial biosphere has experienced dramatic changes in recent decades. Estimates 18 of historical trends in land carbon fluxes remain uncertain because long-term 19 observations are limited on the global scale. Here, we use the Yale Interactive terrestrial 20 Biosphere (YIBs) model to estimate decadal trends in land carbon fluxes and emissions 21 of biogenic volatile organic compounds (BVOCs) and to identify the key drivers for these 22 changes during 1982-2011. Driven with hourly meteorology from WFDEI (WATCH 23 Forcing Data methodology applied to ERA-Interim data), the model simulates an 24 increasing trend of 297 Tg C a -2 in gross primary productivity (GPP) and 185 Tg C a -2 in 25 the net primary productivity (NPP). CO 2 fertilization is the main driver for the flux 26 changes in forest ecosystems, while meteorology dominates the changes in grasslands 27 and shrublands. Warming boosts summer GPP and NPP at high latitudes, while drought 28 dampens carbon uptake in tropical regions. North of 30°N, increasing temperatures 29 induce a substantial extension of 0.22 day a -1 for the growing season; however, this 30 phenological change alone does not promote regional carbon uptake and BVOC 31 emissions. Nevertheless, increases of LAI at peak season accounts for ~25% of the trends 32 in GPP and isoprene emissions at the northern lands. The net land sink shows statistically 33 insignificant increases of only 3 Tg C a -2 globally because of simultaneous increases in 34 soil respiration. Global BVOC emissions are calculated using two schemes. With the 35

Introduction 45 46
The terrestrial biosphere interacts with the atmosphere through photosynthesis and 47 biogenic volatile organic compound (BVOC) emissions. Annually, terrestrial ecosystems 48 assimilate ~120 petagrams of carbon (Pg C) from the atmosphere (Beer et al., 2010), 49 most of which reenters atmosphere through respiration and decomposition, resulting in a 50 net global land carbon sink of 2.6 ± 0.7 Pg C a -1 (Le Quere et al., 2009;Sitch et al., 51 2015). Global BVOC emissions are estimated to be about 1 Pg C per year (Carslaw et al., 52 2010). These emissions are important precursors of atmospheric oxidants and aerosols, 53 both of which affect surface air quality and exert additional regional and global chemical 54 climate forcings (Scott et al., 2014;Unger, 2014). Observations and simulations have 55 shown significant changes in terrestrial carbon assimilation and BVOC emissions in the 56 past 2-3 decades (Lathiere et al., 2006;Sarmiento et al., 2010;Sindelarova et al., 2014;57 emissions or land cover changes due to deforestation, forest management, and 75 agricultural activities (Lathiere et al., 2006;Houghton, 2010). 76 77 Estimates of recent decadal global trends in the land carbon budget and BVOC emissions 78 are limited and uncertain due to the lack of observations. The earliest site-level 79 measurements of land carbon fluxes were set up in the 1990s (Wofsy et al., 1993). The 80 flux tower data sets provide long-term records of regional carbon exchange with high 81 precision but low spatial representation. In contrast, satellite products, such as GPP and 82 NPP retrievals from the Moderate Resolution Imaging Spectroradiometer (MODIS) 83 (Zhao et al., 2005) and isoprene emissions based on tropospheric formaldehyde columns 84 from the Global Ozone Monitoring Experiment (Palmer et al., 2006), improve the spatial 85 coverage but usually are available for only a relatively short time period (months to 86 several years) and suffer from systematic biases when compared with ground 87 measurements (e.g., Heinsch et al., 2006;Marais et al., 2012). Terrestrial biosphere 88 models, evaluated with both site-level and satellite-based observations, are useful tools to 89 estimate trends and attribute drivers of changes in land carbon fluxes and BVOC 90 emissions (e.g., Mao et al., 2013;Stavrakou et al., 2014;Sitch et al., 2015). 91

92
In this study, we use the Yale Interactive Terrestrial Biosphere Model (YIBs, Yue and 93 Unger, 2015) driven with long-term reanalysis meteorology to study the global trends of 94 land carbon fluxes and BVOC emissions over the past three decades. The YIBs model is 95 a process-based vegetation model including complete land carbon cycle (photosynthesis, 96 plant/soil respiration, carbon allocation, and tree growth), plant phenology (Yue et al., 97 2015), and two independent schemes of BVOC emissions (Zheng et al., 2015). Simulated 98 carbon fluxes has been fully validated with carbon fluxes from 145 flux tower sites and 99 multiple satellite products . The major goals of this study are to 100 identify: (1) the dominant drivers of the 30-year trends in carbon fluxes and BVOC 101 emissions from elevated CO 2 , changes in meteorology (temperature, radiation, and soil 102 moisture), and human land use change; (2) the feedback of biosphere, including changes 103 in phenology and leaf area index (LAI), to the trends of land carbon uptakes and BVOC 104 the well-established Farquhar scheme (Farquhar et al., 1980;von Caemmerer and 136 Farquhar, 1981) and the stomatal conductance model of Ball and Berry (Collatz et al., 137 1991). The canopy radiative transfer scheme computes direct and diffuse 138 photosynthetically active radiation (PAR) for sunlit and shaded regions for an adaptive 139 number of layers. The leaf photosynthesis is then integrated over all canopy layers to 140 generate the GPP. 141

142
Part of the assimilated carbon is used for maintenance and growth respiration, and the 143 rest is allocated among different pools for plant development. The model calculates 144 phenology for deciduous forests using cumulative temperature summation with additional 145 constraints from chilling and photoperiod . The phenology of shrubland 146 and grassland is jointly determined by the temperature-and drought-dependent metrics. 147 The LAI is then updated daily based on phenology and the net carbon assimilation. The 148 soil respiration scheme considers carbon flows among 12 biogeochemical pools, 149 including 3 live pools and 9 dead pools. The land carbon source or sink is calculated as 150 the difference between the net carbon assimilation and soil respiration. 151 152 The YIBs model incorporates two independent leaf-level isoprene emission schemes 153 embedded within the exact same host model framework (Zheng et al., 2015). The 154 photosynthesis-based (PS_BVOC) isoprene scheme calculates emissions based on the 155 electron transport-limited photosynthesis rate, canopy temperature, and intercellular CO 2 156 concentrations (Arneth et al., 2007;Unger et al., 2013) simulates an annual GPP of 125 ± 3 Pg C and net ecosystem exchange (NEE) of -2.5 ± 175 0.7 Pg C for 1982-2011, with seasonality and spatial distribution consistent with both 176 satellite observations and benchmark synthesis products . 177 However, the model does not include a fully coupled carbon-nitrogen cycle, which may 178 overestimate CO 2 fertilization effects. In addition, phenology of evergreen trees is set to 179 constant value of 1, leading to underestimation of phenological feedbacks to flux trends. 180 In this study, we use the (2) off-line global version of the model, which is driven with 181 global meteorology reanalysis data and observed CO 2 concentrations. 182 183

Simulations 184 185
We apply observed historical atmospheric CO 2 concentrations from the fifth assessment 186 report (AR5) of the Intergovernmental Panel on Climate Change (IPCC) (Meinshausen et 187 al., 2011). We apply an annually-varying historical transient land cover dataset (Oleson et 188 al., 2013), which is developed based on a combination of remote sensing data from both 189 MODIS (Hansen et al., 2003) and the Advanced Very High Resolution Radiometer 190 (AVHRR) (Defries et al., 2000), and with land use change from Hurtt et al. (2011). We 191 use hourly meteorological variables for 1980-2011 from the WATCH Forcing Data 192 methodology applied to ERA-Interim data (WFDEI, Weedon et al., 2014). The WFDEI 193 reanalysis is an update of the WATCH Forcing Data (WFD), which is developed based 194 on the European Centre for Medium-range Weather Forecasts (ECMWF)  reanalysis (Uppala et al., 2005). Meteorological variables applied include surface air 196 temperature, specific humidity, wind speed, surface pressure, total PAR, and soil 197 temperature and wetness. All of the forcing data are interpolated to the 1°×1° model 198 resolution at the hourly interval. in Europe, northern and eastern Asia, central Africa, and southeastern U.S. in the past 3 226 decades (Fig. 1a). The simulation with year-to-year [CO 2 ], land cover, and meteorology 227 reproduces the magnitude of trend in Europe and the sign of trend in northern Asia, 228 eastern U.S., central Asia, and Australia (Fig. 1b). The model predicts negative changes 229 in central Africa, western U.S., eastern Asia, and the east of South America, which are 230 inconsistent with satellite observations. These negative trends are mainly contributed by 231 the changes in meteorology (Fig. 1e), except for that in East Asia where land cover 232 changes due to human activities result in the decline of LAI (Fig. 1f). Without the land 233 use perturbation, the negative LAI trend in East Asia is weakened and the prediction is 234 closer to observations (Fig. 1c). For the individual drivers, CO 2 fertilization leads to 235 widespread increases in LAI (Fig. 1d), meteorology causes dipole changes on most 236 continents (Fig. 1e), and land use change generally results in negative trends (Fig. 1f). 237 Regionally, simulation CO2_MET_LUC shows a positive trend of 0.0035 m 2 m -2 a -1 in 238 Europe (Table 2), close to the observed value of 0.0049 m 2 m -2 a -1 (Fig. 1a). In other 239 areas, simulated LAI trends are either underestimated (by 87% in Amazon, 78% in North 240 America, and 48% in Central Africa) or opposite in sign (East Asia and Indonesia) 241 compared to observations. Such inconsistencies indicate the limit of model simulations, 242 but may also in part result from the uncertainties in the satellite measurements (see 243 section 4.1). 244 On the global scale, GPP, NPP, and Rh increase respectively by 298, 185, and 181 Tg C 290 a -2 in the past 3 decades (Table 3). The long-term trends of carbon fluxes are mainly 291 driven by CO 2 fertilization, while the interannual variability is related to meteorological 292 forcings (Fig. 5). Warming alone decreases GPP especially in tropical forests (not shown) 293 but increases autotrophic respiration (Ra), leading to global reductions of 56 Tg C a -2 in 294 NPP and 10 Tg C a -2 in NEP (Table 3). Drought alone strongly decreases GPP, especially 295 for tropical grassland and shrubland (Fig. 4), leading to reductions of 51 Tg C a -2 in NPP 296 and 13 Tg C a -2 in NEP. Trends in PAR do not affect GPP and NPP, but may decrease 297 NEP by 23 Tg C a -2 because soil respiration is slowly increasing to reach the equilibrium. 298 Land use change has very limited impacts on the trends of carbon fluxes, though it 299 induces relatively large reductions in NEP (Table 3). 300 America, and East Asia (Fig. 6a). Drought accounts for the decline of isoprene emissions 308 in the U.S. and South America, but land use change is the main driver for the reductions 309 in East Asia (Fig. 6b). Increasing [CO 2 ] promotes photosynthesis but meanwhile inhibits 310 BVOC emissions, leading to offsetting CO 2 effects on isoprene. Consequently, the global 311 isoprene emission is mainly driven by meteorological changes (Fig. 6b). In contrast, 312 using MEGAN scheme, the YIBs model simulates a global reduction of 1.1 Tg C a -2 for 313 isoprene emissions (Fig. 6c). Strong declines are found in the tropical rainforest, for 314 example in the Amazon (-0.43 Tg C a -2 ), central Africa (-0.14 Tg C a -2 ), and Indonesia (-315 0.16 Tg C a -2 ) (Fig. 6c). The MEGAN scheme is sensitive to both light and temperature 316 (Guenther et al., 1995). The strong positive brightening trends in PAR in Europe ( leads to enhancement of regional GPP by 32 Pg C a -1 , much lower than the response of 339 116 Pg C a -1 LAI -1 for the simulation including CO 2 fertilization and climate forcings 340 ( Fig. 7a). In the tropical areas, both positive and negative LAI trends are predicted due to 341 the competition between CO 2 fertilization and drought effects (Fig. 1). As a result, LAI-342 induced GPP and NPP changes show patchy distributions at tropics ( Fig. S2a and S2c), 343 leading to moderate changes in the global carbon assimilations (Table 3). 344 345 Trends in isoprene emission (calculated with the PS_BVOC scheme) also follow that of 346 LAI, except that leaf expansion results in decreased emissions at high latitudes (~60°N, 347 Fig. S2e). The cause for such inconsistency is unclear, but might be because the denser 348 leaves reduce radiation penetrating to lower canopy layers. Such impact would only 349 affect BVOC emissions at high latitudes because PAR is usually limiting near subarctic 350 areas. In most of the subtropical areas, increased LAI leads to enhanced isoprene 351 emissions. On average, unit change in LAI at north of 30°N leads to enhanced isoprene 352 emissions by 43 Tg C a -2 , only 25% of the magnitude in simulation CO2_MET (Fig. 7b). 353 A similar ratio of 23% is achieved for MEGAN isoprene emissions. These results are 354 consistent with that for GPP (Fig. 7a) Plant phenology, which is the timing of budburst and leaf fall, is closely related to 362 temperature, moisture, and photoperiod and thus is experiencing significant changes in 363 the past decades following climate change (Jeong et al., 2011;Keenan et al., 2014;364 Buitenwerf et al., 2015;. Extension of the growing season has the 365 potential to promote carbon uptake of forests (e.g., Piao et al., 2007;Richardson et al., 366 2009). Yet such inference requires careful interpretation because the phenological 367 changes are usually accompanied with warming and elevated [CO 2 ], both of which are 368 also contributing to the enhancement of carbon fluxes. Phenological changes are also 369 expected to affect BVOC emissions, however, such investigations are still missing 370 (Richardson et al., 2013). With the YIBs model, we evaluate the impacts of the growing 371 season extension on both carbon uptake and BVOC emissions by isolating long-term 372 phenological trends from changes in temperature and [CO 2 ]. 373 374 The YIBs model simulates advanced spring and delayed autumn over most areas in NH 375 (Fig. S3). Budburst dates advance on average by 0.16 days a -1 in Europe and 0.15 days a -1 376 in East Asia (Table 2), but with moderate changes or even delays in northwestern Asia 377 and eastern Siberia (Fig. S3a). Spring is earlier by 0.14 days a -1 in eastern U.S. while 378 delayed by 0.15 days a -1 in northwestern U.S. and southeastern Canada, leading to a 379 minor advance of 0.01 days a -1 over North America. Dormancy onset dates are largely 380 delayed in eastern Europe and northwestern Asia (~0.3 day a -1 ), western U.S. (~0.1 day a -381 1 ), boreal Canada (~0.1 day a -1 ), and northeastern China (~0.1 day a -1 ) (Fig. S3b). 382 Advanced autumn (~0.1 day a -1 ) is predicted in northern Asia. Most of these changes are 383 consistent with observations from remote sensing data (Jeong et al., 2011), except for 384 some discrepancies in the magnitude. The predicted phenological trends mainly follow 385 the long-term changes of surface air temperature, especially that in April (for spring) and 386 September (for autumn) (Fig. S4). Sensitivity tests without chilling requirement and 387 photoperiod limit show similar changes , suggesting that temperature 388 changes dominantly drive the trends of forest phenology in the past 3 decades. 389 390 On average, the YIBs model simulates advanced budburst by 0.12 day a -1 and delayed 391 dormancy onset by 0.09 day a -1 at north of 30°N in the past 3 decades (Figs. 8a and 8b). 392 Observations based on remote sensing greenness show trends of -0.11 day a -1 for onset 393 and 0.25 day a -1 for offset during 1990-2009 (Zhu et al., 2013). An ensemble prediction 394 based on 9 terrestrial models yields an advance of 0.08 ± 0.13 day a -1 for onset and a 395 delay of 0.22 ± 0.1 day a -1 for offset (Sitch et al., 2015). Our predictions are in broad 396 agreement with these estimates though the autumn delay is less, likely because the 397 positive trend of offset is weaker for the recent decade (Jeong et al., 2011). September, when both GPP and BVOC emissions and their changes are much smaller 414 compared to that in peak months (Fig. S5). Second, phenological changes are not uniform 415 in space. As Fig. S3 shows, both positive and negative changes are predicted for budburst 416 and dormancy onset dates. Such spatial inhomogeneity, in combination with the 417 discrepancies in regional vegetation types and meteorological conditions, result in varied 418 responses in GPP (Fig. S2b) and isoprene emissions (Fig. S2f). 419 420 Plant phenology at lower latitudes (30°S-30°N) is also experiencing dramatic changes, 421 though such changes are diverse in phase, magnitude, or both (Buitenwerf et al., 2015). 422 In the model, tropical phenology is mainly driven by soil wetness and as a result exhibits 423 large changes in the past 3 decades (not shown). These changes lead to a reduction of 42 424 Tg C a -1 in GPP at the tropics (Fig. S2b), which accounts for 14% of global GPP trend 425 but with the opposite sign (Table 3), suggesting additional inhibition of drought on 426 carbon cycle. A similar conclusion applies for BVOC emissions (Fig. S2f), though 427 experiments suggest that isoprene production has some tolerance to mild drought 428 conditions (e.g., Pegoraro et al., 2006). However, changes in drought-dependent 429 phenology are very uncertain and observations are not available for evaluation. We 430 assume that phenological changes may have larger impacts on both carbon assimilation 431 and BVOC emissions at tropical areas than that at higher latitudes. and 102 ± 34 Tg C a -2 for 1990-2009 based on 9 terrestrial biosphere models (Sitch et al., 466 2015). However, the simulated NPP trend is only 19.8 Tg C a -2 in southern land (15-467 90°S), much lower than the ensemble mean value of 53 ± 31 Tg C a -2 in Sitch et al. 468 (2015). As for the NEP, the YIBs predicts trends of 2.0 Tg C a -2 in northern land, 1.0 Tg 469 C a -2 in tropical land, and -0.3 Tg C a -2 in southern land, much smaller in magnitude 470 compared with the -2.0 ± 12, 36.0 ± 13, and 21 ± 17 Tg C a -2 estimated by Sitch et al. 471 (2015). However, their predictions are insignificant (p > 0.05) for 9, 5, and 7 out of 9 472 models in the northern, tropical, and southern land respectively, suggesting that the 473 strengthening uptake by terrestrial ecosystem is not robust. Their results showed widespread increases in the emissions over China but moderate 478 decreases in Indonesia. In contrast, the YIBs model with the MEGAN scheme simulates 479 widespread reductions in the same areas for 1980-2011 (Fig. 6c). The discrepancies 480 between studies are accounted for by differences in the drivers including land cover 481 change, meteorology, and CO 2 inhibition effects. The YIBs model is driven with land 482 cover data from Hurtt et al. (2011), which estimates an increase of crop (non-isoprene 483 emitter) fraction in East China by 0.32% per year in the last 3 decades, at the cost of the 484 coverage loss by 0.12% a -1 for DBF and 0.14% a -1 for ENF (strong BVOC emitters). 485 However, the data from Ramankutty and Foley (1999), used by Stavrakou et al. (2014) 486 with updates to 2007, show a reduction of the crop fraction over East China for the 487 similar period. In addition, the ERA-Interim PAR used in Stavrakou et al. (2014) shows 488 an increasing trend in southeast China (c.f. their Fig. 5c). On the contrary, the WFDEI 489 PAR for YIBs exhibits a declining trend in the same region (Fig. S1b), leading to a 490 reduction in isoprene emissions. The WFDEI surface solar radiation is based on the ERA-491 Interim radiation but is adjusted using the average cloud cover from the Climatic 492 Research Unit (CRU) and taking into account the effects of interannual changes in 493 atmospheric aerosols (Weedon et al., 2011). Finally, the YIBs simulations include CO 2 494 inhibition effects on BVOC emissions, which were neglected in Stavrakou et al. (2014). previous studies, YIBs with the MEGAN scheme simulates a decreasing trend of ~1 Tg C 505 a -2 in the past 3 decades. The main cause of the discrepancy in the sign of change is the 506 missing CO 2 inhibition effects in the previous studies. In addition, differences in 507 vegetation models, meteorological forcings, and time frames of investigation also likely 508 contribute. The YIBs result is consistent with a recent study by Sindelarova et al. (2014), 509 who reported a decreasing trend of ~1.2 Tg C a -2 for global isoprene emissions during 510 1980-2010 using the MEGAN scheme and inclusion of a CO 2 inhibition parameterization 511 from Heald et al. (2009). 512 513 4.3 Impacts of CO 2 effects 514 515 Similar to the multi-model ensemble predictions (Sitch et al., 2015), we found that global 516 trends in carbon fluxes are dominantly driven by CO 2 fertilization (Figs. 2 and 5). In the 517 YIBs, the global responses to elevated [CO 2 ] is 0.2% ppm -1 for GPP and 0.27% ppm -1 for 518 NPP, with relatively uniform spatial distribution (Figs. S7a and S7b). The GPP response 519 falls within the range of 0.05-0.21% ppm -1 predicted by 10 terrestrial models (Piao et al., 520 2013) and that of 0.01-0.32% ppm -1 observed from multiple free-air CO2 enrichment 521 (FACE) sites (Ainsworth and Long, 2005). The NPP response is higher than the model 522 ensemble of 0.16% ppm -1 (Piao et al., 2013) and the observed median value of 0.13% 523 ppm -1 (Norby et al., 2005), suggesting that CO 2 fertilization to NPP may be 524 overestimated in the YIBs. One possible cause is the omission of N limitation in the 525 model, which could reduce CO 2 responses by half (Piao et al., 2013). Elevated [CO 2 ] 526 leads to increases of 0.023 Pg C a -1 ppm -1 in NEP, within the multi-model range of 0.003-527 0.06 Pg C a -1 ppm -1 (Piao et al., 2013 (Possell et al., 2005), supporting the simulations with 536 MEGAN. In addition, the magnitude of CO 2 inhibition implemented in MEGAN (-0.25% 537 ppm -1 ) is close to observations (-0.26% ppm -1 ) in Possell et al. (2005). However, most of 538 these experiments are conducted for short-term period and cannot detect LAI changes due 539 to the long-term CO 2 fertilization. In addition, the impacts of CO 2 are dependent on 540 species and environmental conditions (ambient temperature and light availability). For 541 example, Buckley (2001) found almost no responses in isoprene emissions to the elevated 542 [CO 2 ] for oak trees. Furthermore, experiments with high temperature and/or light density 543 show increasing isoprene at elevated [CO 2 ] . These studies suggest that 544 the real responses of isoprene emissions to CO 2 under long-term climate change may not 545 be so linear as predicted in MEGAN scheme. More sensitivity experiments and long-term 546 samplings are required to identify CO 2 -isoprene relationships on broad range of biomes 547 and locations. 548 549

Impacts of meteorology 550 551
Predicted long-term trends show large deviations against observations at tropical areas 552 (Fig. 3), where meteorology plays important and complex roles. Responses of carbon 553 fluxes to temperature are more diverse than to CO 2 (Figs. S8a and S8b). In the YIBs, 554 negative responses of GPP and NPP are predicted in tropical areas, where soil moisture 555 availability limits plant functions (e.g. stomatal conductance) to the increased 556 temperature. Furthermore, for tropical rainforests where ambient temperature is higher 557 than optimal photosynthetic temperature (25-30°C), additional warming decreases carbon 558 assimilation, especially for NPP because of simultaneous increases in plant respiration 559 (Liang et al., 2013). On the contrary, warming leads to enhanced GPP and NPP at wetter 560 and cooler areas in the NH subtropics. Such spatial pattern is consistent with multi-model 561 ensemble predictions (Piao et al., 2013). On the global scale, warming results in changes 562 of -0.7% °C -1 for GPP in YIBs, falling within the range of -1.6-1.4% °C -1 estimated by 10 563 models (Piao et al., 2013). Predicted NPP responses of -15-6% °C -1 (Fig. S8b) is not so 564 positive as the measurements of -8-40% °C -1 , probably because most of current warming 565 experiments are located in subtropics of NH (Wu et al., 2011). Elevated temperature 566 changes NEP by -1.4 Pg C a -1 °C -1 , also within the multi-model range of -5~-1 Pg C a -1 567 °C -1 (Piao et al., 2013). Simulated isoprene emissions with PS_BVOC show similar 568 warming responses as that of carbon fluxes (Fig. S8c), except for tropical rainforests 569 where the former is positive while the latter is negative. Such decoupling is attributed to 570 the differences in optimal temperatures between isoprene (35-40 °C) and photosynthesis 571 (25-30 °C). Simulations with MEGAN scheme show very strong temperature dependence 572 of 6-15% °C -1 (Fig. S8d), consistent with measurements of 5-20% °C -1 for aspen 573 (Niinemets and Sun, 2015) and 9-12% °C -1 for oak (Li et al., 2011). However, 574 experiments with some other species (e.g. spruce in Kivimaenpaa et al. (2013)) show no 575 responses or moderate ones, suggesting that warming sensitivity of isoprene emissions 576 might be dependent on species and ambient conditions. 577 578 Responses to PAR are mostly positive and distributed evenly, with global sensitivity of 579 0.3% W -1 m 2 for GPP and 0.5% W -1 m 2 for NPP (Figs. S9a and S9b). Isoprene emissions 580 from both PS_BVOC and MEGAN schemes show similar responses to PAR, with larger 581 sensitivity in subtropics than that in tropics (Figs. S9c and S9d), likely because the 582 ambient PAR is higher at lower latitude, leading to slower responses of isoprene 583 emissions to the unit changes of light (Guenther et al., 1993). YIBs simulations show that 584 PAR is not the driver of long-term trends in carbon fluxes and BVOC emissions (Fig. 4), 585 likely because changes in solar radiation is limited in the past 3 decades (Figs. S1b). 586 587 Soil moisture dominates climate-driven flux changes in tropical areas (Fig. 4). In YIBs 588 model, changes in soil water availability affect carbon assimilation through the alteration 589 of leaf stomatal conductance and plant phenology (especially for shrublands and 590 grasslands in arid regions). Both GPP and NPP show strong responses to soil wetness 591 variations, especially over tropics where >10% changes are found for every increase of 592 0.01 in soil wetness at 1.5 m (Figs. S10a and S10b). On the global scale, GPP changes by 593 4.7% 0.01 -1 and NPP by 5.5% 0.01 -1 in response to soil wetness. Although experiments 594 also show rapid reductions in carbon assimilation due to drought stress (e.g., Ruehr et al., 595 2012;Xia et al., 2014), the magnitude of such influence is difficult to evaluate because 596 different metrics and depths of soil water are used in measurements. Isoprene emissions 597 from PS_BVOC show similar soil-wetness responses to that of GPP (Fig. S10c), 598 indicating that drought reduces BVOC emissions. However, observations show 599 insignificant changes of isoprene with mild drought stress (e.g., Pegoraro et al., 2006), 600 though such drought tolerance is strongly weakened at severe drought and/or warm 601 conditions (Centritto et al., 2011). Consistent with these experiments, MEGAN scheme 602 does not include drought inhibition on isoprene emissions. Simulations with YIBs show 603 large responses of BVOC to soil wetness in tropical areas (Fig. S10d) and Europe. The afforestation in Europe helps promote regional carbon uptake, resulting 611 in more reasonable trends in LAI compared with remote sensing data (Fig. 1). However, 612 the expansion of crop in China leads to a reduction in LAI, which is not supported by the 613 satellite data. One possible cause is the uncertainty in crop fraction, because data from 614 Hurtt et al. (2011), used by YIBs, show crop expansion while data from Ramankutty and 615 Foley (1999) suggest reductions of the crop fraction over East China over the similar 616 period. The role of land use change in our simulation might be conservative because we 617 consider only land cover changes. Perturbed emissions from land use management, such 618 as forest lodging, cropping practice, use of fertilizer, fire management and so on 619 (Houghton, 2010) may alter regional carbon budget by changing carbon sinks to sources. 620 Studies including gross emissions of land use perturbation estimated a global net land 621 source to atmosphere, which shows decreasing trend in the last 3 decades (Ciais et al., 622 2013). Such change may help strengthen net land carbon sink but is missing in our study. 623 The land biosphere has experienced significant changes in the past 3 decades. At north of 627 30°N, changes in LAI account for 25% of the trends in regional carbon fluxes and 628 isoprene emissions. However, the extension of growing season alone makes insignificant 629 contributions to the increased carbon assimilation. This conclusion is inconsistent with 630 site-level observations that show evident increases in carbon assimilation at early spring 631 and/or late autumn in recent decades (Dragoni et al., 2011;Keenan et al., 2014). The 632 causes for such discrepancies lie in two. First, phenology at specific location may exhibit 633 much more intense changes than that at larger scale. For example, Dragoni et al. (2011)  contributes to the stronger uptake at early spring and late autumn. One difficulty for the 639 observation-based estimate of phenological impacts is that extension of growing season is 640 accompanied by warmer climate, which may stimulate both carbon assimilation and 641 BVOC production. In a recent study, Barlow et al. (2015) found invariant length of land 642 carbon uptake period at high northern latitudes based on the first time differential of 643 atmospheric CO 2 concentrations, suggesting that increased greenness is not necessarily 644 equal to enhanced carbon uptake in shoulder seasons. Furthermore, Barlow et al. (2015) 645 showed that enhanced peak uptake is the main driver for the strengthened carbon sink at 646 high northern latitudes over the past 4 decades. These conclusions are supportive of our 647 simulations for the monthly trends at subtropical regions (North America, Europe, and 648 East Asia) (Fig. S5). Increase of temperature promotes carbon uptake of forest ecosystems at high latitudes 657 (>30°N) while drought tendency dampens GPP and NPP of grasslands and shrublands at 658 low latitudes (30°S-30°N). The widespread increases of LAI at northern lands account for 659 ~25% of the regional GPP trends. Significant changes in phenology are found at north of 660 30°N; however, this temperature-driven phenological change alone is not promoting 661 regional carbon assimilation. Changes in land use show limited influences on global 662 carbon fluxes, except for some regional impacts over Europe (afforestation)   The spatial distribution of GPP and isoprene changes is shown in Figure S2. and meteorology (red), temperature only (magenta), and phenology only (blue). Units of 1187 trends are (a) day a -1 , (b) day a -1 , (c) Pg C a -1 day -1 , and (d) Tg C a -1 day -1 . The spatial 1188 distribution of GPP and isoprene changes is shown in Figure S2.