Land cover change dominates decadal trends of biogenic volatile organic compound (BVOC) emission in China

Land cover change dominates decadal trends of biogenic volatile organic compound (BVOC) emission in China Hui Wang1, 2, Qizhong Wu1, Alex B. Guenther2, Xiaochun Yang1, Lanning Wang1, Tang Xiao3, Jie Li3, Jinming Feng4, Qi Xu1, Huaqiong Cheng1 1College of Global Change and Earth System Science, Joint Center for Global Changes Studies, Beijing Normal University, 5 Beijing 100875, China 2Department of Earth System Science, University of California, Irvine, CA 92697, USA 3State Key Laboratory of Atmospheric Boundary Layer Physics and Atmospheric Chemistry, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China 4Key Laboratory of Regional Climate-Environment for Temperate East Asia, Institute of Atmospheric Physics, Chinese 10 Academy of Sciences, Beijing 100029, China Correspondence to: Qizhong Wu (wqizhong@bnu.edu.cn) & Lanning Wang (wangln@bnu.edu.cn) Abstract. Satellite observations reveal that China has been leading the global greening trend in the past two decades. We assessed the impact of land cover change on total BVOC emission in China during 2001-2016 and found a significant increasing trend of 1.09% yr-1 with increases of 1.35, 1.25 and 1.43 % yr-1 for isoprene, monoterpenes and sesquiterpenes, 15 respectively. Comparison of different scenarios showed that vegetation change is the main driver of BVOC emission change in China. Considerable heterogeneity was observed on regional scales, with the highest increasing trends of BVOC emission found in the Qinling Mountains and in the south of China. The BVOC emission for the year 2016 in these two regions was enhanced by 61.89 and 67.64% compared to that of 2001, respectively. We compared the long-term HCHO vertical columns (VC) from the satellite-based Ozone Monitoring Instrument (OMI) with the estimation of isoprene emission in summer. The 20 results showed statistically significant positive correlation coefficients over the regions with high vegetation cover fractions. In addition, the isoprene emission and HCHO VC both showed statistically significant increasing trends in the south of China where these two variables have high positive correlation coefficients. This result supports our estimation of the variability and trends of BVOC emission in China. Although anthropogenic sources comprise ~63% NMVOC emissions in China, the continued increase of BVOC will enhance the importance of considering BVOC when making policies for controlling ozone 25 pollution in China along with ongoing efforts to reduce anthropogenic emissions.

air quality and the climate system. The emission of BVOC is controlled by multiple environmental factors like temperature, radiation, concentration of CO2 and other stresses, and is affected by climate changes (Guenther et al., 1995;Arneth et al., 2007;Penuelas and Staudt, 2010). Besides the climatic factors, the land cover change also plays a key role in the variability of BVOC emission (Stavrakou et al., 2014;Unger, 2014;Chen et al., 2018), e.g., cropland expansion has been estimated to dominate the reduction of 5 isoprene, the dominant BVOC species, in last century (Lathière et al., 2010;Unger, 2013) although there are large uncertainties associated with these estimates.
China has been greening in recent decades (Piao et al., 2015). A recent study points out that China accounts for 25% of the net increase of global leaf area during 2000-2017 (Chen et al., 2019). The increase of forest area plays a dominant role in greening in China with multiple programs to maintain and expand 10 forests (Zhang et al., 2016;Bryan et al., 2018;Chen et al., 2019). The enhancement of vegetation cover rate and biomass can lead to the increase of BVOC emission in China and induce a corresponding impact on local air quality and the climate system. Previous studies have investigated the long-term emission trend of dominant BVOC species like isoprene in China (Fu and Liao, 2012;Li and Xie, 2014;Stavrakou et al., 2014;Chen et al., 2019). Li and Xie (2014)  There are no long-term in-situ observation of BVOC in China to validate our estimation of interannual 5 variability of BVOC emission, however, satellite formaldehyde (HCHO) observations provide an opportunity to validate the interannual variability of isoprene, the dominant compound among BVOC species that accounts for almost half of total BVOC emission in China . Since HCHO is an important proxy of isoprene in forest regions with no significant anthropogenic impact, satellite observed HCHO columns are widely used to derive regional ecosystem isoprene emission (Palmer et al., 10 2003;Marais et al., 2012;Stavrakou et al., 2015;Kaiser et al. 2018). Zhu et al. (2017b) reported the increasing trend of HCHO vertical columns (VC) detected by Ozone Monitoring Instrument (OMI) driven by increasing cover rate of local forest in the northwestern US. Stavrakou et al. (2018) also used the longterm HCHO VC to investigate the annual variability of BVOC induced by climate variability. We have used the long-term HCHO record from 2005-2016 by OMI to assess our estimation of annual isoprene 15 variability.
where Fi, ei and gi represent the emission amount, standard emission factor and emission activity factor of chemical species i. The standard emission factor in this study is based on the plant functional type (PFT) distribution, and the PFT scheme in MEGAN v2.1 is the scheme adopted in Community Land Model 4.0 (Lawrence et al., 2011). The emission activity factor gi accounts for the impact of multiple environmental factors and expresses it as: 5 where gp, i, gT, i gA, i gSM, i and gC, i represent the activity factors for light, temperature, leaf age, soil moisture and CO2 inhibition impact. The Cce (=0.57) is a factor to set the gi equal to 1 at the standard conditions. The LAI is the leaf area index and defines the amount of foliage and the leaf age in MEGAN. The light and temperature response algorithms in MEGAN v2.1 are from Guenther et al. (1991), Guenther et al. (1993) 10 andGuenther et al. (2012), which described enzymatic activities controlled by temperature and light conditions. The CO2 inhibition algorithm is from Heald et al. (2009), and only the estimation of isoprene emission considers the impacts of soil moisture and CO2 concentration. The detailed descriptions of these factor algorithms can be found in Guenther et al. (2012) and Sakulyanontvittaya et al. (2008).

15
The land cover parameters for driving MEGAN including LAI, PFT and vegetation cover fraction (VCF) were provided by satellite datasets . The MODIS MOD15A2H for 2001 and MCD15A2H for 2002-2016 LAI datasets were adopted in this study. The parameter LAIv in MEGAN is calculated as: where VCF is provided by MODIS MOD44B datasets. 20 The PFT was used to determine the canopy structure and standard emission factors in MEGAN (Guenther et al., 2012). The PFT data in this study is obtained from the MODIS MCD12C1 land cover product, and MODIS PFT classification were converted to MEGAN PFT classification based on the climatic criteria described by Bonan Gordon et al. (2002) using the climatology of ERA-interim dataset (Berrisford et al., https://doi.org/10.5194/acp-2020-28 Preprint.

Meteorological Datasets
The hourly meteorological fields were provided by the Weather Research and Forecast (WRF) Model V3.9 (Skamarock et al., 2008) simulations. The meteorological simulation is driven by ERA-Interim 5 reanalysis data (Berrisford et al., 2011) with 27 km horizontal spatial resolution and 39 vertical layers.
The physical schemes were presented in supplemental Table S1.
Since light and temperature conditions are the main environmental drivers of BVOC emission (Guenther et al., 1993;Sakulyanontvittaya et al., 2008), we assessed the reliability of the WRF simulated downward shortwave radiation (DSW) and 2-meter temperature (T2) using the in-situ observations from 98 radiation 10 observation sites and 697 meteorology observation sites in China. The in-situ observations used in this study are from the National Meteorological Information Center (http://data.cma.cn/). We converted the hourly model outputs and daily observations to monthly average values from 2001 to 2016 for comparison.
For DSW, the average mean bias (MB), mean error (ME) and root mean square error (RMSE) are 40.37 (± 20.81), 43.55 (± 17.52) and 49.79 (± 17.70) W m-2 among 98 sites, and the overestimation of DSW 15 simulation is a common issue in multiple simulation studies and may be induced by the lack of aerosol radiation affect and cloud simulation (Wang et al., 2011;Situ et al., 2013;Wang et al., 2018). For T2, the average MB, ME and RMSE are -1.19 (± 2.87), 2.40 (± 2.14) and 2.65 (± 2.11) °C among 697 sites over China. We also compared the monthly anomalies of DSW and T2 from the model simulation and observation to validate the interannual variability of meteorological fields simulated by WRF. As shown 20 in Figure 1, the results indicate that the model accurately reproduced the interannual variability of DSW and T2, and the correlation coefficients of DSW and T2 anomaly between the simulation and observation reached 0.77 and 0.88, respectively. Our WRF simulation successfully captured the long-term meteorological variabilities and is reasonable to use for estimating the impact of climatic variability on BVOC emission in China for this study.

Satellite Formaldehyde Observation
The satellite HCHO VC used in this study is from the Belgian Institute for Space Aeronomy (BIRA-IASB) and retrieved using the differential optical absorption spectroscopy (DOAS) algorithm (De Smedt et al., 2012;De Smedt et al., 2015). The detailed description of the BIRA-IASB OMI HCHO product can be found in De Smedt et al. (2015), and we used the monthly HCHO VC with 0.25° ´ 0.25° spatial resolution. 5 Since the OMI instrument is temporally stable (Dobber et al., 2008;De Smedt et al., 2015), the OMI HCHO VC product is suitable for long-term analysis (Jin and Holloway, 2015) and was used to primarily validate our estimation of isoprene emission variability. The major sources of tropospheric HCHO are biogenic VOC, anthropogenic source and open fires (Zhu et al., 2017a). Since biogenic isoprene is the dominant source of HCHO in the forest regions without obvious anthropogenic impact (Palmer et al., 10 2003), we used HCHO as the proxy of isoprene to validate the interannual variability of isoprene estimates.

Scenarios and Analysis Method
We designed four scenarios (S1-S4) to investigate the impact of land cover change and climatic conditions on BVOC emission. The configurations of the four scenarios are shown in Table 1, and S1 was considered as the standard or "full" scenario with both annually updated land cover parameters (LAIv and PFT) and The climatic variability can affect the growth of vegetation and then affect LAI values (Piao et al., 2015).
In this study, the interaction between climate and ecosystem is not considered in the offline MEGAN model, which means that the meteorological conditions, e.g. precipitation, will not affect the LAI values.
The LAI input for MEGAN model in this study was obtained from the remote sensed LAI products. 25 Therefore, when the time of the meteorological condition is inconsistent with that of the LAI input, the indirect impact of meteorological conditions on BVOC emission through biomass and phenology were https://doi.org/10.5194/acp-2020-28 Preprint. Discussion started: 16 March 2020 c Author(s) 2020. CC BY 4.0 License. neglected. We used the experiments with the inconsistent LAI and meteorological conditions to investigate the direct effect of climatic variability on BVOC emission.
The chemical species emissions estimated by MEGAN were grouped into four major categories including isoprene, monoterpene, sesquiterpene and other VOCs since the terpenoids account for the majority of total BVOC emission and have known impacts on atmospheric oxidants and SOA (Wang et al., 2011). 5 The trend analysis in this study was done following the Theil-Sen trend estimation method and the results were tested by the Mann-Kendall non-parametric trend test, and the corresponding results in this study were calculated using the trend_manken (https://www.ncl.ucar.edu/Document/Functions/Builtin/trend_manken.shtml) function of the NCAR Command Language (NCL, https://www.ncl.ucar.edu/).

The Variability of BVOC Emission in China During 2001-2016
As shown in Table 2, the average annual emission during 2001-2016 of isoprene, monoterpene, sesquiterpene and other VOCs estimated from S1 are 7.56 (±0.74), 1.37 (±0.12), 0.16 (±0.02) and 6.73 (±0.46) Tg, respectively. Isoprene is the dominant species and accounts for about half of the total BVOC emission in China. S1 is the standard scenario that includes both annually updated meteorological fields 15 and vegetation conditions. In comparison with previous studies (Table 3), our estimation of isoprene emission is very close to the results by Stavrakou et al. (2014) and Tie et al. (2006) while our estimation of monoterpene emission is considerably (57 to 72%) lower than other estimations. Multiple factors including interannual variations, horizontal resolution, meteorological and land cover inputs can lead to the discrepancy of these estimations. 20 The total estimated BVOC emission has a statistically significant increasing trend with rates of 1.09 and 1.19% yr-1 for the S1 and S2 scenarios (d in Figure 2), respectively. The increasing rate of isoprene, monoterpene and sesquiterpene are 1.35, 1.25 and 1.43% yr-1 respectively for the S1 scenario. In comparison, the increasing rates of these species in the S2 scenario are higher than those in the S1 scenario with 1.58, 1.51 and 1.61 % yr-1 for isoprene, monoterpene and sesquiterpene despite the direct impact of 25 meteorological conditions. Although the S1 scenario considers the impact of annual meteorological variability, the BVOC emission is still in a significant upward trend driven by the increasing forest area meteorological variability. Furthermore, as shown in Figure 2, the largest discrepancy in total BVOC emission in S3 and S4 appears between the year 2014 and 2016, with the total BVOC emission in the year 2016 being 21.1 % and 22.5 % higher than that of the year 2014 in S3 and S4, respectively. In general, the interannual variability of meteorological conditions leads to ~20% difference in BVOC emission during our study period. In contrast, the CV of total BVOC in S2 is 5.74%, which is close to that of 5% 5 in S3 and S4, showing that interannual variability is dominated by meteorology even though the trend is dominated by landcover. Moreover, the largest discrepancy of total BVOC in S2 (23.9%) occurred between the year 2002 and 2016 and is very close to that estimated solely for meteorological conditions. However, as mentioned above, the comparison here only considered the direct impact of meteorological conditions, and the meteorological conditions also can affect the growing process and phenology which 10 can influence BVOC emission indirectly (Peñuelas et al., 2009). Considering the direct and indirect impact of climatic conditions as well as land cover change, the CV of total BVOC in the "full" scenario S1 is 8.15%, which is higher than the other scenarios. The highest and lowest total BVOC emission in S1 are in year 2016 and 2010, respectively, with 2016 being 34.60% (5.00 Tg) higher than 2010. These results show that both landcover and meteorology can individually contribute ~20%, and together over 15 30%, to the estimated annual variability in China BVOC emissions within this 16-year time period.

The Regional Variability of BVOC Emission in China
The hotspots of BVOC emission are mainly located in the northeast, central and south of China where the forest is widely distributed and the climate is warm and favorable for emitting BVOC as shown in Figure   3. The Changbai Mountains, Qinling Mountains, the southeast and southwest China forest regions, 20 southeast Tibet and Taiwan are the regions with highest BVOC emission in 2001, which is broadly consistent with the previous estimations (Tie et al., 2006;Li et al., 2013).
The spatial distributions of statistically significant (p>0.9) changing trends in S1-S4 are also presented for individual categories in Figure 3. In general, the spatial distributions of trends of different species in S1 and S2 are highly consistent since the vegetation development is the main driver of the increasing strong negative trend is found at the boundary of Jiangxi and Hunan provinces. While a positive increasing trend induced by meteorology is also found in Tibet, western Sichuan and southeastern Yunnan province in S3 and S4, and it is clear that most of the trend hotspots in S2 do not overlap with those in S3 and S4, which further indicates that the vegetation development is the main driver of BVOC increasing trend in S1, the "full" scenario, rather than meteorological conditions. average emission rate during the study period appears in 2010 because of the lowest cover fraction of broadleaf trees. However, in S1 with the impact of meteorological variability, the average emission rate 20 in 2016 is 5.89 g m-2 yr-1. This is 7.48% higher than the average emission rate in 2001 of 5.48 g m-2 yr-1, but the lowest average emission rate of 4.49 g m-2 yr-1 is still in 2010. In general, it is clear that land cover change is a dominate factor impacting the interannual variability of BVOC emission on a decadal scale in regions undergoing rapid landcover change as has been suggested by previous studies (Unger, 2014;Chen et al., 2018). 25 The estimated increase of BVOC in regions like the Qinling Mountains and southern China are expected to affect regional air quality. For the Qinling Mountains and surrounding areas, as estimated by Li et al. (2018) using the WRF-chem model, the average contribution of BVOC to O3 could reach 16.8 ppb for the https://doi.org /10.5194/acp-2020-28 Preprint. Discussion started: 16 March 2020 c Author(s) 2020. CC BY 4.0 License. daily peak concentration and 8.2 ppb for the 24h concentration in the urban region of Xi'an, one of the biggest cities near the Qinling Mountains suffering poor air quality in recent years (Yang et al., 2019). For southern China, Situ et al. (2013) reported that BVOC emission could contribute an average 7.9 ppb surface peak O3 concentration for the urban area in the Pearl River Delta region, and the contribution from BVOC even reached 24.8 ppb over PRD in November. Since BVOC plays an important role in local air 5 quality, the change of BVOC emission may have an even greater effect on the local ozone pollution. For instance, the simulation study by Li et al. (2018) also found that the urban region of Xi'an is VOC-limited because of the abundant NOx emission there. Therefore, the increase of BVOC emission in the Qinling Mountains would further favor the formation of O3 in the urban region of Xi'an.

10
The lack of long-term in-situ observations of BVOC in China makes it difficult to validate the variability and trend estimation of BVOC emission in this study. However, since biogenic isoprene is the dominant precursor of formaldehyde in rural regions with minimal anthropogenic influence (Palmer et al., 2003),  The grid level correlation coefficients between the average summer HCHO VC and isoprene emission estimated in our study are shown in Figure 5d, and the grids with statistically significant correlations (p>0.9, N=12) grids are marked with black dots. A positive correlation can be found in the northeast, Plain have an increasing trend because of the increase of human activities (Smedt et al., 2010), there is also an obvious increasing trend of HCHO VC in the developed Yunnan and Guangxi provinces in the south of China. Moreover, these regions, especially Guangxi province also show a statistically significant positive correlation between isoprene emission and HCHO VC as presented in Figure 5d. This indicates that biogenic emissions might be the main driver of the increased HCHO in Guangxi province. 25

Comparison of BVOC Emission with Anthropogenic Emission in China.
China has initiated a series of pollution control policies in recent years (Zheng et al., 2018;Ma et al., 2019) and achieved success in controlling some air pollutants, especially PM2.5 Yu et al., 2019;Xu et al., 2019). However, the ozone pollution is still severe in China especially the mega-city areas (Li et al., 2019;Xu et al., 2019). It is widely known that NOx and VOC are the precursors of ground NOx/VOC emission would correspondingly promote the formation of ozone (Tang et al., 2010;Li et al., 2018;Li et al., 2019). Both AVOC and BVOC show the increasing trend since 2010, as shown in Figure   6, with an increase of 9.7% for AVOC and 34.6% for BVOC from 2010 to 2016, respectively. Combined with the decreasing trend of NOx emission, the overall changes of precursors are expected to make it more difficult to control ozone pollution. In addition, there are multiple studies pointing out the 20 interactions between anthropogenic emission and biogenic emission on ozone pollution in mega cities including Beijing (Pang et al., 2009;Shao et al., 2009), Shanghai (Geng et al., 2011), Guangzhou (Situ et al., 2013) and Xi'an . The current trends in biogenic as well as anthropogenic emissions suggest that biogenic emissions will play an even more important role in future air pollution in China.

25
Satellite observations have shown that China has led the global greening trend in recent decades (Chen et al., 2019). In this study, we investigated the impact of this greening trend on BVOC emission in China The increase of BVOC reported by this study is expected to lead to a more complex situation for making the policies for controlling ozone pollution in China. The recent pollution control policies in China have effectively initiated the control of PM2.5 pollution, but the ozone pollution is still severe especially in 5 urban areas Yu et al., 2019;Xu et al., 2019;Li et al., 2019). Although anthropogenic emission is still the dominant source of NMVOC in China and is ~70 % higher than the average biogenic emission in China, the BVOC still makes an important contribution to ozone pollution in mega cities including Beijing (Pang et al., 2009;Shao et al., 2009), Shang Hai (Geng et al., 2011), Guang Zhou (Situ et al., 2013) and Xi'an  and may further increase in importance considering the continuing 10 greening trend over China in the future.

Author Contribution
QW, LW and HW planned and organized the project. HW, JF and QX prepared the input datasets. HW modelled and analyzed the data. HW and QW wrote the manuscript. HW, AG and QW revised the manuscript. AG, XY, LN, XT, JL, JF and HC reviewed and provided key comments on the paper. 15

Data Availability
The source code of MEGAN v2.