Regional to global distributions, trends, and drivers of biogenic volatile organic compound emission from 2001 to 2020

Biogenic volatile organic compounds (BVOCs) are important precursors to ozone and secondary organic aerosols in the atmosphere, affecting air quality, clouds, and climate. However, the trend in BVOC emissions and driving factors for the emission changes in different geographic regions over the past 2 decades has remained unclear. Here, regional to global changes in BVOC emissions during 2001–2020 are simulated using the latest Model of Emission of Gases and Aerosols from Nature (MEGANv3.2) with the input of time-varying satellite-retrieved vegetation and reanalysis meteorology data. Comparison of model simulations with the site observations shows that the model can reasonably reproduce the magnitude of isoprene and monoterpene emission fluxes. The spatial distribution of the modeled isoprene emissions is generally comparable to the satellite retrievals. The estimated annual average global BVOC emissions are 835.4 Tg yr⁻¹ with the emissions from isoprene, monoterpenes, sesquiterpenes, and other BVOC comprised of 347.7, 184.8, 23.3, and 279.6 Tg yr⁻¹, respectively. We find that the decrease in global isoprene emissions (−0.07 % per year) caused by the increase in CO₂ concentrations (−0.20 % per year) is stronger than that caused by changes in vegetation (−0.03 % per year) and meteorological factors (0.15 % per year). However, regional disparities are large. Isoprene emissions increase significantly in Europe, East Asia, and South Asia (0.37 % per year–0.66 % per year). Half of the increasing trend is contributed by increased leaf area index (LAI) (maximum over 0.02 m² m⁻² yr⁻¹) and tree cover. Changes in meteorological factors contribute to another half, with elevated temperature dominating in Europe and increased soil moisture dominating in East and South Asia. In contrast, in South America and Southeast Asia, shifts in vegetation type associated with the BVOC emission capacity, which partly results from the deforestation and agricultural expansion, decrease the BVOC emission and offset nearly half of the emission increase caused by changes in meteorological factors. Overall, isoprene emission increases by 0.35 % per year and 0.25 % per year in South America and Southeast Asia, respectively. In Central Africa, a decrease in temperature dominates the negative emission trend (−0.74 % per year). Global monoterpene emissions show a significantly increasing trend (0.34 % per year, 0.6 Tg yr⁻¹) compared to that of isoprene (−0.07 % per year, −0.2 Tg yr⁻¹), especially in strong greening hotspots. This is mainly because the monoterpene emissions are more sensitive to changes in LAI and are not subject to the inhibition effect of CO₂. The findings highlight the important roles of vegetation cover and biomass, temperature, and soil moisture in modulating the temporal variations of global BVOC emissions in the past 2 decades.

1 Introduction

Emissions of biogenic volatile organic compounds (BVOCs) from the terrestrial vegetation play a pivotal role in atmospheric chemistry and climate due to their large quantity (∼ 1000 Tg yr⁻¹) and high reactivity (Guenther et al., 1995, 2012). Isoprene and monoterpenes (e.g., α-pinene, β-pinene, limonene) are the most prevalent BVOC species, and other species include sesquiterpenes, methanol, and ethanol. These BVOCs account for about 90 % of total non-methane volatile organic compounds (NMVOCs) in the atmosphere (Guenther et al., 2006), which are important precursors for troposphere ozone and secondary organic aerosols (SOAs) through atmospheric oxidation processes and thus influence air pollution, clouds, and Earth's radiative budget. However, BVOC emissions remain highly uncertain as they depend on a diversity of factors, and the relative importance of different factors is still unclear.

The BVOC emissions are determined by many environmental factors such as vegetation, meteorology, and carbon dioxide (CO₂) concentrations. The impact of vegetation on BVOC emissions is primarily reflected in vegetation types (e.g., forest, grassland), tree species, and vegetation biomass density (e.g., land cover, leaf area index). Previous studies have demonstrated that different vegetation types and tree species affect BVOC emissions dramatically (Lathière et al., 2006; Stavrakou et al., 2014; Sindelarova et al., 2022). For meteorological parameters, especially temperature, light, and soil moisture, it was observed that elevated temperature, stronger radiation, and wetter soil significantly promote BVOC emissions (Rinne et al., 2002; Bai et al., 2016; Jiang et al., 2018). In contrast, elevated CO₂ can suppress the emissions of the major BVOC component (e.g., isoprene) (Heald et al., 2009; Wilkinson et al., 2009). The vegetation change has affected nearly half of the global land surface (Hurtt, 2011) with land cover change intensifying in recent decades, particularly in tropical and East Asia (Purves et al., 2004; Pacifico et al., 2012; Piao et al., 2015). Related studies pointed out that one-third of the global vegetation growth area has been greening since the 21st century, with the leaf area index (LAI) increasing by 2.3 % per decade, of which China contributes nearly one-quarter (Chen et al., 2019). In addition, the surge in greenhouse gas emissions since the industrial revolution has led to significant global warming and meteorological changes. Thus, due to the large spatio-temporal variations in these factors, there may exist significant regional and global differences in BVOC emissions as well as emission trends.

Ground-based measurements can sample BVOC fluxes from the leaf to the canopy scale (Müller et al., 2010; Bai et al., 2016; Sarkar et al., 2020). However, the spatial and temporal coverages of such observations are limited and cannot be extended to represent BVOC emissions in a larger domain. The BVOC emission models, such as the Model of Emissions of Gases and Aerosols from Nature (MEGAN, Guenther et al., 1995, 2006, 2012) and the BVOC Photosynthesis-Dependent Scheme (PS_BVOC, Arneth et al., 2007; Unger et al., 2013), consider the main factors influencing BVOC emissions. Compared to the widely used MEGAN version 2.1, the recently released MEGAN version 3.2 has addressed several gaps in BVOC emission modeling (Guenther et al., 2020), including (i) more refined BVOC emission factors, where vegetation emission factors are calculated based on tree species rather than the original fixed plant functional type (PFT) emission factor, and (ii) consideration of environmental stress caused by extreme weather and air pollution. These models have been applied to simulate the BVOC emissions from regional to global scales and to identify the impacts of various factors on the variations in BVOC emissions (Fu and Liao, 2012; Stavrakou et al., 2014; Chen et al., 2018; Wang et al., 2021; Li et al., 2022).

Previous studies have explored the long-term emission trends in BVOCs. Chen et al. (2018) paid attention to changes in vegetation from 2000 to 2015 on the global scale, and reported that the global total isoprene emission declined by only 1.5 % for this period. Opacka et al. (2021) complemented the work of Chen et al. (2018) by incorporating different satellite-retrieved land cover datasets and pointed out that land cover changes from 2001 to 2016 mitigate isoprene emissions ranging from −0.33 % to −0.04 % per year, while temperature and radiation changes enhance isoprene emissions by 0.94 % per year. Sindelarova et al. (2022) argued that although vegetation change from 2000 to 2019 exerted a small effect (−0.11 % per year) on the overall change in global isoprene emissions, the changes were significant in some hotspots. Purves et al. (2004) found that changes in vegetation in the eastern United States from 1980 to 1990 could lead to an increase in BVOC emissions by about 17 %. Fu and Liao (2012) analyzed the changes in BVOC emissions in China from 2001 to 2006 and found that the interannual variability of isoprene emissions was dominated by changes in the meteorological fields, while the variability of monoterpene emissions was more sensitive to changes in vegetation. Stavrakou et al. (2014) explored the factors influencing isoprene emissions in Asia from 1979 to 2012 and found that enhanced temperature and light led to a 0.52 % per year increase in emissions in China, and oil palm expansion in Southeast Asia increased the isoprene emissions by more than 1 % per year. Chen et al. (2018) showed that from 2000 to 2015, afforestation caused a 5 %–10 % increase in isoprene emissions in northeastern China and India, while deforestation led to about a 10 % reduction in isoprene emissions in the Amazon Basin, West Africa, and Southeast Asia. Wang et al. (2021) indicated that the greening of China from 2001–2016 led to a significant increase (up to 11.7 %) in BVOC emissions, and the regional accumulated BVOC emissions in 2018 could be 26 % higher than those in 2001, mainly due to the increase in vegetation cover and LAI (Li et al., 2022). However, these previous studies were conducted at various spatio-temporal scales using fixed vegetation or meteorological datasets for a given year, leading to difficulties in comparing their magnitude and even sign of the BVOC emission trends, and it is still unclear which factor dominates the BVOC emission trends in different hotspot regions.

This study provides a comprehensive analysis of BVOC emission trends from 2001 to 2020 on a regional to global scale using the latest BVOC emission model, MEGANv3.2, combined with time-varying satellite vegetation retrievals and meteorological reanalysis data. More importantly, this study further identifies the contribution of various driving factors to these trends. The findings of this research shed light on the importance of BVOC emissions in air quality and aerosol radiative forcing in hotspot regions. This paper is organized as follows. In Sect. 2, we introduce the MEGANv3.2 model and its input data, model experimental design, and BVOC observation data used in model evaluations. In Sect. 3, the spatio-temporal distributions and trends in BVOC emissions in different regions are analyzed. The contributions of various driving factors (i.e., vegetation, meteorology, and CO₂) are quantified. In Sect. 4, we discuss the uncertainty of the model results and the sources of biases in the BVOC emission trends by comparing with previous studies. Conclusions of this study are given in Sect. 5.

2 Method and data

2.1 MEGANv3.2

The MEGAN emission model has been widely used to simulate BVOC emissions at global and regional scales (Guenther et al., 2006, 2012). It has also been incorporated into various Earth system models and chemical transport models (Müller et al., 2008; Li et al., 2013; Sindelarova et al., 2014; Messina et al., 2016; Bauwens et al., 2018; Chen et al., 2018). Here, the latest version of MEGANv3.2 is applied to estimate the BVOC emissions from 2001 to 2020. Compared to the earlier version MEGANv2.1 (Guenther et al., 2012), MEGANv3.2 estimates vegetation emission factors based on variable plant species measurements instead of on fixed plant functional type (PFT; Guenther et al., 2020). Specifically, while MEGANv2.1 uses a look-up table of emission factors for the 15 PFTs corresponding to the biological emission classes (see Table 2 in Guenther et al., 2012), MEGANv3.2 uses the so-called emission factor processor, to estimate the landscape average emission factors, which are based on the following three databases: (1) growth form datasets for four PFTs (tree, shrub, grass, and crops); (2) ecotype dataset, composed of a mix of emission-specific tree species and grass associated with specific emission capacities; and (3) updated tree species and grass datasets corresponding to the biogenic emission classes. These updates can distinguish the differences in vegetation emission factors in regions with the same PFT but with varying plant species. The new version also considers the additional stress factors of emissions by using the simple threshold function, including high and low temperature, strong wind, and heavy O₃ pollution. Additionally, the number of BVOC components in MEGANv3.2 is expanded from the original 148 to over 200. MEGANv3.2 calculates the BVOC emissions rate (ER) as follows:

$$\begin{array}{}\text{(1)}& \mathrm{ER}=\mathrm{EF}\cdot \mathrm{EA},\end{array}$$

where EF and EA represent the standard emission factor (i.e., ER at “standard” conditions) and the nondimensional emission activity factor, respectively.

The EF map can be obtained by running the MEGANv3.2 emission factor processor, which combines growth form and ecotype data with plant species community and species emission factor datasets to generate the mean emission factor. A brief algorithm is shown below:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{EF}& ={\mathrm{EF}}_{\mathrm{tree}}\cdot f_{\mathrm{tree}}+{\mathrm{EF}}_{\mathrm{shrub}}\cdot f_{\mathrm{shrub}}+{\mathrm{EF}}_{\mathrm{grass}}\cdot f_{\mathrm{grass}}\\ & +{\mathrm{EF}}_{\mathrm{crop}}\cdot f_{\mathrm{crop}},\end{array}\end{array}$$

where EF_tree, EF_shrub, EF_grass, and EF_crop represent the species emission factors for the four types of growth forms (i.e., PFTs) and f is the fraction of the specific growth form in a model grid cell.

The emission activity factor considers the effect of various environmental factors and is calculated as

$$\begin{array}{}\text{(3)}& \mathrm{EA}=\mathrm{LAIv}\cdot {\mathit{\gamma }}_p\cdot {\mathit{\gamma }}_T\cdot {\mathit{\gamma }}_{\mathrm{HT}}\cdot {\mathit{\gamma }}_{\mathrm{LT}}\cdot {\mathit{\gamma }}_{\mathrm{SW}}\cdot {\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{3}}}\cdot {\mathit{\gamma }}_A\cdot {\mathit{\gamma }}_{\mathrm{SM}}\cdot {\mathit{\gamma }}_C,\end{array}$$

where LAIv represents the leaf area index of vegetation-covered surfaces and is obtained by dividing LAI with VCF (vegetation cover fraction). γ_p, γ_T, γ_HT, γ_LT, γ_SW, ${\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{3}}}$, γ_A, γ_SM, and γ_C represent the activity factors for downward shortwave radiation, 2 m air temperature, high temperature, low temperature, strong wind, O₃ pollution, leaf age, soil moisture, and CO₂ concentration, respectively. In MEGANv3.2 the increases in temperature, radiation, and soil moisture favor the BVOC emissions, while high/low temperature (> 40 or < 10 °C) and strong wind (> 12 m s⁻¹) as well as heavy O₃ pollution and high CO₂ concentration suppress the BVOC emissions. Currently, the inhibition effect of CO₂ on BVOC emissions in the model is only available for isoprene. In this study, we consider the effect of all above activity factors on BVOC emissions except for the O₃ pollution factor (i.e., ${\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{1}$). More details of these algorithms are described in Guenther et al. (2006, 2012, 2020).

2.2 Data

2.2.1 Vegetation datasets

The vegetation parameters driving MEGANv3.2 include LAI, VCF, and PFT. In this study, the Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation retrievals from 2001 to 2020 were used. LAI data were obtained from Yuan et al. (2011), which improved the MODIS Version 6 product MCD15A2H (Myneni et al., 2015) with a temporal resolution of 8 d and a spatial resolution of 0.5° × 0.5°. The LAIv calculated in MEGANv3.2 is defined as LAI divided by VCF, representing the leaf area index per unit vegetation area. The VCF was from the yearly MODIS MOD44B Version 6 dataset (DiMiceli et al., 2015), which contains three ground cover components (i.e., tree cover, non-tree cover, and bare soil cover). We summed the first two observed variables as vegetation cover. The raw VCF product (250 m pixel size, sinusoidal grid) was further converted to a 0.5° × 0.5° latitude/longitude grid by a conservation interpolation before being used in MEGANv3.2. The PFT was obtained from the yearly MODIS MCD12C1 product with a spatial resolution of 0.05° (Friedl and Sulla-Menashe, 2015). To calculate the growth form fractions used by the emission factor processor in Eq. (2), the selected 17 MODIS IGBP (International Geosphere Biosphere Programme) global vegetation classification types from the above MODIS PFT product were mapped to four main PFT classification types (i.e., tree, shrub, grass, and crop) in MEGANv3.2 based on methods from Sulla-Menashe and Friedl (2018). The ecotype dataset mentioned in Sect. 2.1 is based on satellite imagery and ground surveys and comes from the MEGAN development group (https://bai.ess.uci.edu/megan/data-and-code/growth-form-and-ecotypes, last access: 21 November 2022) The reprocessed datasets were conservatively interpolated to a spatial resolution of 0.5° × 0.5° as model inputs.

2.2.2 Meteorological datasets

The meteorological parameters driving the MEGANv3.2 model were from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) (Gelaro et al., 2017). To reduce the biases in the simulated radiation flux trends from aerosols, MERRA-2 has assimilated aerosol optical depths from space-based observations in the long-term and considered the interaction of aerosols with the climate system. The selected variables used in MEGANv3.2 include 2 m temperature, surface downward shortwave radiation, surface soil moisture, water vapor mixing ratio, 10 m wind speed, precipitation, surface air pressure, low-level wind speed, cloud cover, and snow cover. Photosynthetically active radiation (PAR) in MEGANv3.2 was obtained by dividing the surface downward shortwave radiation by 2. The temporal resolution of these variables is either 1-hourly or 3-hourly, and the 3-hourly data are linearly interpolated to the uniform 1-hourly data. All selected parameters were further interpolated from the original 0.5° × 0.625° to a spatial resolution of 0.5° × 0.5° (consistent with the resolution of the vegetation datasets) for driving the MEGANv3.2 model. In our study, MERRA-2 data from 2001 to 2020 were used.

2.2.3 Observations

We used in situ observations of BVOC emission fluxes collected from the literature for the comparison with our model estimates. In total, isoprene emission fluxes at 26 observation sites and monoterpene emission fluxes at 11 observation sites were collected and listed in Tables S1 and S2 in the Supplement. Other BVOC emission fluxes were not collected, mainly due to few observations of them and their small contributions to the total BVOC burdens. The units of all collected data were converted to mg C m⁻² d⁻¹ for easy comparison.

In addition, we also employed global isoprene burden data from space-based observations for the comparison with the simulation. The space-based observations from the Cross-track Infrared Sounder (CrIS) include direct retrievals of the global isoprene column burden with an optimal estimate for January, April, July, and October of 2013 (Fu et al., 2019; Wells et al., 2020). Although isoprene column burdens are different from emission fluxes, there is a strong positive correlation between them, and thus isoprene burden data can provide a good reference for qualitative analysis of the spatial distribution of isoprene emissions.

2.3 Simulations

To isolate the contribution of different influencing factors (vegetation, meteorology, and CO₂) to BVOC emission trends from 2001 to 2020, we perform nine sensitivity experiments (Table 1). These experiments consist of two groups. The first group contains four experiments: EMIT_ALL is the control experiment that considers the historical changes in all factors. EMIT_VEG, EMIT_MET, and EMIT_CO₂ consider only the historical changes in vegetation parameters, meteorological factors, and CO₂ concentration, respectively, while the other factors are fixed as those in 2001. These three experiments were used to quantify the contributions of vegetation, meteorology, and CO₂ concentrations to the BVOC emission trends, respectively. In the second group, five experiments were conducted to isolate the contributions of individual vegetation parameters (i.e., PFT, LAIv) and meteorological factors (i.e., temperature, light, and soil moisture). For vegetation parameters, the experimental setup is the same as EMIT_VEG but with PFT (EMIT_VEG_FIX_PFT) or LAIv (EMIT_ VEG_FIX_LAIv) fixed as that in 2001. The difference between EMIT_VEG and EMIT_VEG_FIX_PFT/EMIT_VEG_FIX_LAIv represents the contribution of LAIv and PFT historical changes to the BVOC emission trends. For meteorological factors, the experimental setup is the same as EMIT_MET but with temperature (EMIT_MET_FIX_T2m), light (EMIT_MET_FIX_RAD), or soil moisture (EMIT_MET_FIX_SM) fixed as that in 2001. The difference between EMIT_MET and EMIT_MET_FIX_T2m, EMIT_MET_FIX_RAD, and EMIT_MET_FIX_SM represents the impact of temperature, light, and soil moisture changes to the BVOC emission trends, respectively. The model horizontal resolution is 0.5° × 0.5°, the temporal resolution is 1 h, and the simulation period is 2001–2020. The input variables include the satellite-retrieved vegetation parameters and MERRA-2 reanalysis data as described above.

Table 1 Description of model experiments driven with vegetation parameters, meteorological parameters, and CO₂.

LAIv: leaf area index of vegetation-covered surfaces; PFT: plant functional type; T2m: 2 m temperature; RAD: surface solar radiation; SM: soil moisture; CO₂: CO₂ concentrations.

3 Results

3.1 Spatio-temporal distribution of BVOC emissions

3.1.1 Spatial distribution of BVOC emissions

The latest version of the MEGANv3.2 model used in this study has not been fully evaluated in previous studies; thus we first compared the simulation results with the in situ isoprene and monoterpene emission fluxes collected from the literature. For a fair comparison, the simulated isoprene and monoterpene fluxes were interpolated to the sample-specific locations and averaged over the same measurement period.

Figure 1a shows the locations of the observation sites, which are mainly distributed in North America, Europe, and Asia. Four main vegetation types are included: (1) evergreen broadleaf forest (EBF), (2) deciduous broadleaf forest (DBF), (3) evergreen needleleaf forest (ENF), and (4) grassland (grass). Observations show that there is a large range of isoprene emission fluxes (10–120 mg C m⁻² d⁻¹) and monoterpene emission fluxes (0.1–15 mg C m⁻² d⁻¹) for different PFTs (Guenther et al., 2012). Isoprene emission flux is generally larger in EBF and DBF than in ENF and grass. Compared to the observations of isoprene emission fluxes, MEGANv3.2 can simulate the magnitude of isoprene emission fluxes with a correlation coefficient of 0.48 and a mean bias of −3.25 mg C m⁻² d⁻¹ (Fig. 1b). The model can also simulate larger isoprene emission fluxes for EBF and DBF than for ENF and grass. However, the model tends to underestimate isoprene fluxes from grassland by a factor of about 10, possibly because the prescribed grass emission factors are too low. The model also tends to mostly overestimate the isoprene emission fluxes from EBF while underestimating them from DBF. For ENF, although there is only one station, the model significantly overestimates the isoprene emission fluxes by a factor of 6.

Fewer monoterpene samples (11) were collected than isoprene samples (26). As shown in Fig. 1c, the model overestimates the monoterpene emission with a correlation coefficient of 0.34 and a mean bias of 3.65 mg C m⁻² d⁻¹. Although the overestimation is mostly within a factor of 10, it largely overestimates the monoterpene emission flux at two sites located in eastern China and northwestern South America by 2 orders of magnitude. Comparison of long-term observations at the K34 tower site (vegetation type: EBF) in 2013 in the Amazon reveals that the simulated seasonal variations of isoprene emission fluxes are similar to the observation (Fig. 1d). The model also captures the increase in emission during the dry season with a correlation coefficient of 0.48, although there is a smaller contrast between dry season and other seasons in the model.

Not surprisingly, there are still large discrepancies between the simulations and observations, which may be ascribed to the deficiency in emission parameterizations such as vegetation emission factors. It may be also ascribed to input parameters such as vegetation and meteorology. In addition, note that the comparison of in situ isoprene measurements with the model is not always representative. Isoprene has a very short lifetime (minutes to hours) in the atmosphere, which implies that its measured fluxes depend on the local tree or plant species near the observation site. The model results at a horizontal resolution of 50 km represent a regional mean for the grid cell where the specific observation site is located, which can also partly explain the difference between the simulations and observations because of the large spatial variability in BVOC emissions. The large bias in the seasonal variation of isoprene fluxes in MEGANv3.2 may be due to a lack of representation of the isoprene emission capacity of tree species at different leaf ages (Alves et al., 2018). Additionally, the model bias arises from a lack of realistic representations of leaf phenology, canopy structure, soil moisture feedbacks, and variation in isoprene emissions due to the complex biodiversity in the Amazon region.

Figure 1 (a) The location distribution of isoprene and monoterpene observation sites based on the literature collection (the five-pointed star represents the location of the K34 tower site; EBF: evergreen broadleaf forest; DBF: deciduous broadleaf forest; ENF: evergreen needleleaf forest; grass: grassland) and comparison with simulated (b) isoprene fluxes, (c) monoterpene fluxes, and (d) monthly variations of isoprene fluxes measured at the K34 tower site in 2013 (the upper and lower limits of the observations represent 1 standard deviation). R: correlation coefficient; bias: absolute bias; RMSE: root mean square error; Num: total number of observation samples. In (b) and (c), the 1:1 (solid) and $\mathrm{1}:\mathrm{10}/\mathrm{10}:\mathrm{1}$ (dashed) lines are plotted for references. In (d), gray shade denotes the dry season.

Figure 2 shows the comparison of simulated isoprene emission fluxes from this study, the IASB-TD-OMI (top–down isoprene emissions constrained by HCHO columns from the Ozone Monitoring Instrument) dataset, and CrIS. IASB-TD-OMI datasets employ a “top–down” approach to constrain isoprene emission fluxes simulated by a chemical transport model with formaldehyde observations from the Ozone Monitoring Instrument (OMI) (Stavrakou et al., 2014, 2015). CrIS is directly retrieved from satellite observations (Sect. 2.2). Note that CrIS provides the isoprene burden, which cannot be directly compared with the modeled emission flux. Here we use it to indicate the relative intensity of isoprene emission (Sect. 2.2) and validate the spatial patterns of simulated isoprene emissions.

The simulated results and IASB-TD-OMI show similar spatial patterns of isoprene emission, with the correlation coefficients of 0.64, 0.77, 0.74, and 0.77 for January, April, July, and October, respectively. However, except for January, the simulated emissions from this study are systematically higher than those obtained from IASB-TD-OMI by about 20 %. The differences are mainly concentrated in South America, Central Africa, and Southeast Asia, which may be partly due to differences in estimation methods (i.e., top–down vs. bottom–up), emission model parameters (e.g., vegetation emission factors), and the meteorological datasets used (MERRA-2 vs. ERA-Interim). In addition, the simulated isoprene emissions are compared to the isoprene column burdens retrieved from the CrIS. The spatial distributions of simulated isoprene emission and the CrIS-retrieved column burden are moderately correlated, with the correlation coefficient varying between 0.43 and 0.58. The isoprene emission and the column burden show similar spatial patterns (e.g., large values in eastern North America and Central Africa) and seasonal variations (e.g., the lowest emission flux and column burden in July), which suggests that the model estimates reasonable isoprene emissions.

Figure 2 The comparison of (left column) simulated isoprene emission fluxes (0.5° × 0.5°) with (middle column) IASB-TD-OMI (0.5° × 0.5°) and (right column) CrIS data (2° × 2.5°) in January, April, July, and October of 2013. The annual global emission or burden is given at the bottom left of each panel. The correlation coefficients R of the simulation results with IASB-TD-OMI and CrIS are marked in the panels in the middle and right columns, respectively. Gray shade denotes the regions where data are not available.

Figure 3 shows the spatial distribution of annual emission fluxes for four BVOC categories, including isoprene, monoterpenes, sesquiterpenes, and other BVOCs. The contribution of each category to total BVOC emission is also shown. The spatial patterns of the emission fluxes for the four BVOC categories are relatively similar (Fig. 3a–d). The strongest emissions (> 20 g m⁻² yr⁻¹ for total BVOC emission) are mainly located in the tropical regions such as South America, Central Africa, and Southeast Asia. The latitudinal distribution of emissions shows that the peaks of all four categories are located at 0–15 °S, decreasing gradually toward the poles (Fig. 3f). Isoprene emission accounts for the largest fraction (>40 %) of total BVOC emission at 0–15° S, and the fraction decreases at higher latitudes. In contrast, monoterpene emission accounts for the smallest fraction (around 20 %) at 0–15° S, and the fraction increases to ∼40 % at high altitudes. The latitudinal variations of different BVOC categories are mainly caused by different kinds of PFT cover at different latitudes.

Figure 3 The spatial pattern of the annual average BVOC emission fluxes of (a) isoprene, (b) monoterpenes, (c) sesquiterpenes, (d) other BVOCs, (e) total BVOC, and (f) the latitudinal distribution (left: total emission; right: fraction of each category; pie chart: the fraction of each category in global total emissions) for the period 2001–2020. The annual global emissions are given at the bottom left of each panel.

In total, annual global BVOC emission is 835.4 Tg yr⁻¹ during 2001–2020. The total emission for isoprene, monoterpenes, sesquiterpenes, and other BVOC is 347.7, 184.8, 23.3, and 279.6 Tg yr⁻¹, respectively, and accounts for 41.6 %, 22.1 %, 2.8 %, and 33.5 %, respectively, of the total BVOC emission. Overall, isoprene is the dominant component for total BVOC emissions. The contribution from sesquiterpenes is relatively minor compared to isoprene and monoterpene emissions. Table 2 further compares the BVOC emissions estimated in this study with others. The annual global BVOC emission (835.4 Tg yr⁻¹) from this study is within the range (558–1005 Tg yr⁻¹) of previous estimations with MEGANv2.1 using different drivers. The annual global isoprene emission estimated in this study (347.7 Tg yr⁻¹) is smaller than most previous bottom–up estimates but is close to the top–down estimate based on GOME2 (344.7 Tg yr⁻¹). Interestingly, the isoprene emission estimated in this study is 12 % lower than that estimated by Weng et al. (2020), while the monoterpene emission estimated in this study is 26 % higher, although the present study used the same driving fields (i.e., MERRA-2 meteorological field and MODIS vegetation parameters) as Weng et al. (2020). These discrepancies are mainly ascribed to the differences in vegetation emission factors between the two versions of MEGAN. As compared to MEGANv2.1, isoprene emission factors are smaller and monoterpene emission factors are larger in MEGANv3.2. Note that MEGANv2.1 only utilizes fixed emission factors corresponding to the PFTs, but the PFT is insufficient to characterize the emission factors; e.g., tree species with the same PFT may have very different BVOC emission rates. MEGANv3.2 further considers differences in emission factors for tree species with the same PFTs. Thus, the vegetation emission factors in MEGANv3.2 are more accurately represented. However, we note that the uncertainties associated with emission factors are still large due to the limited observational data (Guenther et al., 2020). The estimated sesquiterpene emission in this study (23.3 Tg yr⁻¹) is consistent with previous estimates based on the MERRA meteorological field (20.0–21.6 Tg yr⁻¹) but is higher than those based on the ERA (the European Center for Medium-Range Weather Forecasting reanalysis products) meteorological field (11.9–16.6 Tg yr⁻¹). This is partly due to the higher values of 2 m temperature and downward shortwave radiation field in MERRA than in ERA in the tropical regions, resulting in higher emission (Sindelarova et al., 2022). In addition, the other BVOC emission estimated in this study (279.6 Tg yr⁻¹) is higher than in previous studies (115.5–278.8 Tg yr⁻¹), which can be partly attributed to more BVOC components (e.g., butane, butanenitrile, acetophenone, benzene cinnamaldehyde, cinnamic acid) considered in MEGANv3.2 than in MEGANv2.1.

Table 2 Comparison of annual global BVOC emission rate with previous studies (Tg yr⁻¹).

Note: MEGAN – Model of Emission of Gases and Aerosols from Nature; HEMCO – Harvard-NASA Emissions Component; MACC – Monitoring Atmospheric Composition and Climate project; MOHYCAN – Model of Hydrocarbon emissions by the CANopy; LPJ-GUESS – Lund-Potsdam-Jena General Ecosystem Simulator; BVOCEM – Biogenic Volatile Organic Compound Emission Model; OMI – Ozone Monitoring Instrument; GOME – Global Ozone Monitoring Experiment instrument; GEOS-Chem – the global 3-D model of atmospheric chemistry driven by meteorological input from the Goddard Earth Observing System (GEOS) of the NASA Global Modeling and Assimilation Office (GMAO); MERRA – the Modern-Era Retrospective analysis for Research and Applications; ERA5 – the European Center for Medium-Range Weather Forecasting fifth generation of atmospheric reanalysis products; ERA-Interim – the European Center for Medium-Range Weather Forecasting interim reanalysis products; GEOS-FP – GEOS-Chem met field archive of the GMAO “forward processing” product; Qian 2006 – Qian et al. (2006) atmospheric forcing; CRU – the Climatic Research Unit of the University of East Anglia; NCEP – the National Center for Environmental Prediction reanalysis product; UM – climate model output from the UK Met Office Unified Model; GEOS-STRAT – met field data product compatible with GEOS-Chem from GMAO; MODIS – Moderate Resolution Imaging Spectroradiometer; CLM4 – Community Land Model; ESA-CCI – Climate Change Initiative of the European Space Agency; CLM-SP – standard global simulation constrained by observed land cover; GFW – Global Forest Watch; LPJV – Lund-Potsdam-Jena vegetation; SDGVMV – Sheffield Dynamic Vegetation Model; AVHRR – Advanced Very High Resolution Radiometer.

3.1.2 Seasonal variation of isoprene emissions

The seasonal cycle of simulated isoprene emissions over the nine regions (i.e., NAM: North America; EUR: Europe; NAS: North Asia; EAS: East Asia; SAS: South Asia; SEAS: Southeast Asia; SAM: South America; CAF: Central Africa; AUS: Australia) is presented in Fig. 4. Overall, different regions show different seasonal cycles. Isoprene emissions mainly depend on vegetation and meteorological conditions. At midlatitudes to high latitudes, due to the densest vegetation and highest temperatures in the summer months (Fig. S1 in the Supplement), isoprene emissions peak in the summer months. In contrast, isoprene emissions are lowest during winter months in these regions. For example, in North America, Europe, North Asia, and East Asia (Northern Hemisphere), isoprene emissions are highest in July and lowest in January. We note that the seasonal variations of isoprene emissions in Europe are consistent with the results from the EMEP (European Monitoring and Evaluation Programme) model, which are derived based on the cover fractions and emission factors of detailed tree species and other vegetation (categorized into six PFTs) (Sindelarova et al., 2022). In Australia (Southern Hemisphere), isoprene emissions peak in austral summer. At low latitudes of northern China, South Asia, and South America isoprene emissions peak at the end of the dry season in May and October, respectively, due to the higher temperatures (Fig. S1). In the tropical regions and low latitudes of the Southern Hemisphere, including Southeast Asia and Central Africa, isoprene emissions exhibit a typical bimodal distribution with two peaks in April and October, corresponding to the peak LAI and temperature in these 2 months (Fig. S1).

In total, the global total isoprene emission peaks in July and reaches its minimum in January. This is different from the study of Sindelarova et al. (2014), which showed the peak of global total isoprene emission in October and the lowest value in June, mainly because Sindelarova et al. (2014) simulated greater isoprene emissions over South America and Central Africa than our study.

3.2 Trends in BVOC emissions

Figure 5 shows the interannual variation of total isoprene emissions over the nine regions for the last 20 years from the four experiments in the first group (ALL: EMIT_ALL; VEG: EMIT_VEG; MET: EMIT_MET; CO₂: EMIT_CO₂). By comparing these experiments, the individual impact of three major factors (i.e., vegetation, meteorology, and CO₂) on isoprene emissions can be identified. There are strong interannual variations in isoprene emissions in all the regions. As shown in Table 3, the standard deviation of the regionally summed isoprene emissions is 0.1–8.1 Tg yr⁻¹ in the nine regions, and the ratio of standard deviation to the climatology mean is 0.01–0.08. For the trends, although the trend in global isoprene emissions between 2001 and 2020 is weak (∼0.07 % per year; p=0.67), significant trends in regional emissions can be found in some specific regions. In Europe, East Asia, and South Asia, regionally summed isoprene emissions exhibit significantly increasing trends of about 0.37 % per year–0.66 % per year. In contrast, there is a significantly decreasing trend in Central Africa (−0.74 % per year). The trends are weak in the other regions (North America, South America, Southeast Asia, and Australia) except in North Asia, where a decreasing trend of −0.32 % per year is modest but not statistically significant due to the large interannual variability of isoprene emission in this region.

For the impacts of various factors, the interannual variations of isoprene emissions are mainly determined by meteorological factors in all the regions, as the time series of isoprene emissions in MET follows closely those in ALL. In comparison, vegetation factors play a much smaller role in the interannual variations of isoprene emissions. As CO₂ concentrations increase constantly, there is little impact from CO₂ concentrations on the interannual variations of isoprene emissions. The most remarkable impacts of vegetation and CO₂ concentrations is on the trends in isoprene emissions during the last 20 years.

Vegetation factors lead to significant trends in isoprene emissions (p<0.1 for AUS; p<0.05 for other regions) in all regions. In Europe, North Asia, East Asia, South Asia, and Australia, vegetation changes (increases in vegetation cover or shift in vegetation type) strongly promote isoprene emissions, with the trends of 0.28 %–0.41 % per year. In contrast, vegetation changes (decreases in vegetation cover or shift in vegetation type) reduce isoprene emissions in South America and Southeast Asia (with the trends of about −0.15 % per year), which may be ascribed to the local deforestation.

Figure 4 The seasonal variation of isoprene emissions for each region outlined in black on the map averaged during 2001–2020. The shaded area represents 1 standard deviation. The nine regions are as follows: NAM – North America; EUR – Europe; NAS – North Asia; EAS – East Asia; SAS – South Asia; SEAS – Southeast Asia; SAM – South America; CAF – Central Africa; AUS – Australia.

Figure 5 Interannual variation and trends in regional isoprene emissions for each region outlined by black rectangles on the map. Trends are expressed by the relative change in percentage (i.e., the linear change during 2001–2020 divided by the mean value). p values denote statistically significant levels using the Mann–Kendall test. ALL represents simulated results considering interannual variability of all drivers (i.e., vegetation, meteorology, and CO₂), while VEG, MET, and CO₂ represent simulated results considering only interannual variability of vegetation, meteorology, and CO₂ concentrations, respectively.

Table 3 The statistical parameters (avg: average value, Tg yr⁻¹; SD: standard deviation, Tg yr⁻¹; SD / avg: ratio of standard deviation to average value) of total isoprene emissions over the nine regions during 2001–2020 from the four experiments in the first group (ALL, VEG, MET, and CO₂).

There are also considerable contributions from meteorological factors to the trends in isoprene emissions. Changes in meteorological factors lead to significantly increasing trends in isoprene emissions (0.28 % to 0.43 % per year; p<0.1) in Europe, East Asia, South Asia, and South America. In contrast, meteorological factors cause a strongly declining trend in isoprene emissions in Central Africa (−0.61 % per year; p<0.05). There are also increasing trends in isoprene emissions (0.25 % to 0.35 % per year) in North America and Southeast Asia and decreasing trends (−0.15% to −0.38 % per year) in North Asia and Australia due to changes in meteorological factors, although the trends are not statistically significant (p>0.1).

As expected, the CO₂ inhibition effect results in a significant decrease of −0.20 % per year in isoprene emissions (Fig. 5). Note that this study uses a globally uniform and yearly mean CO₂ concentration without considering spatial and seasonal variations of CO₂ concentration. Additionally, the CO₂ concentration can also indirectly affect isoprene emissions by changing meteorological and vegetation factors. Specifically, the increase in CO₂ concentration is partly responsible for global warming and thus higher temperatures (e.g., Monson et al., 2007), and the increased CO₂ concentration can potentially lead to a larger LAIv (Monson et al., 2007). These indirect effects are not explicitly considered in this study.

Different from isoprene emissions, there is no statistically significant effect of CO₂ concentration on monoterpene emissions as suggested by previous studies (Malik et al., 2019, 2023). Therefore, monoterpene emissions in MEGANv3.2 only consider the effects of vegetation and meteorological factors and show a significantly positive trend of 0.34 % per year globally (Fig. S2). In total, the increasing trends caused by vegetation changes and meteorological factors are comparable globally. In Europe, East Asia, and South Asia, the increasing trends induced by vegetation changes (increase in vegetation cover or shift in vegetation type) are larger than those induced by meteorology changes, indicating the dominant role of vegetation change in these regions. In North Asia and Australia, vegetation changes lead to increasing trends, while meteorological changes lead to decreasing trends. The overall impacts are the increase in monoterpene emissions, due to the dominant impacts from vegetation changes. In the other regions, changes in meteorological factors dominate vegetation changes in determining the trends in monoterpene emissions. In particular, in South America and Southeast Asia, changes in meteorological factors lead to strongly increasing trends in monoterpene emissions, which are much larger than those induced by vegetation changes.

Vegetation changes affect isoprene and monoterpene emissions differently, mainly due to differences in the emission capacity (i.e., emission factors, Eq. 2) of vegetation types for these two components. Vegetation changes lead to a stronger trend in the monoterpene emission than that of the isoprene emission in some regions. For example, the isoprene-emission trends are 0.08 % per year, 0.40 % per year, 0.28 % per year, 0.30 % per year, and 0.41 % per year, while the monoterpene-emission trends are up to 0.19 % per year, 0.54 % per year, 0.43 % per year, 0.47 % per year, and 0.56 % per year in North America, Europe, North Asia, East Asia, and South Asia, respectively. However, in Southeast Asia and South America, the effect of vegetation changes on the trends in monoterpene emissions is weaker. In addition, the effect of meteorological changes on monoterpene emissions is found to be weaker than that on isoprene emissions in most regions, especially in Central Africa, where the trend for isoprene emission is −0.61 % per year, while the trend for monoterpene emission is −0.26 % per year. Overall, compared to isoprene emissions, a stronger increasing trend and a weaker decreasing trend in monoterpene emissions are found in these hotspot regions when all influencing factors are considered.

3.3 Drivers of BVOC emission trends

3.3.1 Changes in vegetation factors, meteorology, and CO₂ concentrations

BVOC emissions depend on various factors including vegetation parameters, meteorological conditions, and CO₂ concentration. Figure 6 shows the global distribution of the trends in these influencing factors of BVOC emissions from 2001 to 2020, including vegetation parameters (i.e., VCF, LAI), meteorological parameters (i.e., surface 2 m temperature, surface solar radiation, and soil moisture), and CO₂ concentrations. Different regions show different VCF trends. A moderate increase of 0.1 % to 1 % per year in VCF can be found in central North America, East Asia, and India, while a decrease of −0.1 % to −1 % per year exists in central South America, Central and southwestern Africa, the west of Australia, and Central Asia. The decrease in VCF is mainly related to local deforestation or wildfire burning. Specifically, the four main PFTs used in MEGANv3.2 (i.e., tree, shrub, grass, and crop; Eq. 2) have shifted significantly in some regions (Fig. 7a–d). Tree cover increases in East Asia and Europe by about 0.3–0.5 % per year. Grass cover increases in central South America, with a maximum increasing trend of more than 1 % per year, corresponding to the decrease in tree and shrub cover in these regions. Crop cover also increases in some regions such as eastern South America and northern India.

Figure 6 Trends in (a, b) vegetation factors, (c–e) meteorological factors, and (f) CO₂ concentrations from 2001 to 2020. (a) Vegetation cover fraction (VCF), (b) leaf area index (LAI), (c) surface 2 m temperature (T2m), (d) surface solar radiation, (e) soil moisture, and (f) CO₂ concentration. Stippling denotes regions where the trend is statistically significant (p<0.1) using the Mann–Kendall test.

Most vegetated areas are becoming greener (i.e., have a higher LAI) during 2001–2020, especially in Europe, East Asia, and South Asia, with a positive trend exceeding 0.02 m² m⁻² yr⁻¹. The decreasing trends in LAI only exist in some isolated regions in eastern South America, Central Africa, and the west of Australia. Note that VCF and LAI reflect two different aspects of vegetation information using different retrieval methods. VCF represents the amount of ground covered by the vertical projection of vegetation, and its variation focuses on the extension or shrinkage of vegetation on the ground surface. While LAI reflects the amount of vegetation biomass, its variation emphasizes the changes in vegetation biomass content per unit of ground surface area. Therefore, the variations of these two parameters may be different. LAIv in MEGANv3.2 (Eq. 3) is calculated by dividing LAI by VCF, representing the leaf area per unit of the canopy area (vegetation-covered area), which implicitly reflects the overall growth condition of vegetation (e.g., the amount of plant foliage, canopy structure) per unit of vegetation-covered area. Due to the rapid shrinking of vegetation cover, LAIv in South America shows a more remarkable increasing trend than LAI does (Fig. 7e).

For meteorological factors, trends in surface 2 m temperature (T2m) vary significantly in different regions. T2m increases with a small trend (0.01 to 0.05 °C yr⁻¹) in most regions during 2001–2020. A larger increase of 0.05 to 0.2 °C yr⁻¹ occurs in the Arctic and Europe. Despite increasing trends in T2m in most regions, there are also cooling trends of −0.05 to −0.2 °C yr⁻¹ in some regions such as northeastern North America, Central Africa, and Central Asia. Large trends in surface solar radiation exist in tropical regions, with significant dimming trends in central South America and Central Africa and brightening trends in some parts of Southeast Asia. Soil moisture increases significantly in large parts of Central Africa, South Asia, and northeastern Asia, while it decreases significantly in South America and Australia. Finally, global CO₂ levels have increased dramatically from 370.57 ppm in 2001 to 412.44 ppm in 2020. Overall, these factors have changed significantly over the past 20 years, which has led to significant changes in BVOC emissions as shown in Sect. 3.2.

Figure 7 Trends in (a)–(d) cover fractions of the four plant functional types (PFTs) and (e) LAIv (leaf area index of vegetation-covered surfaces) from 2001 to 2020. (a) Tree, (b) shrub, (c) grass, and (d) crop. Stippling denotes regions where the trend is statistically significant (p<0.1) using the Mann–Kendall test.

3.3.2 Contribution of drivers to BVOC emission trends

We further decompose the main drivers (vegetation, meteorology, and CO₂ concentrations) and quantify the contribution of individual drivers to isoprene emission trends. Figure 8 shows the spatial patterns of isoprene emission trends estimated from the nine experiments listed in Table 1. The pattern in Fig. 8a is affected by a combination of vegetation change (Fig. 8b), meteorological variability (Fig. 8c), and CO₂ concentration change (Fig. 8d). The statistically significant and positive trends are mainly located in the midlatitudes of Eurasia, with the strongest trends of 2 % to 5 % per year, while the negative emission trends are distributed in Central Africa and central Australia, with the largest magnitudes ranging from −5 % to −8 % per year. Vegetation changes increase isoprene emissions in most regions except in the Amazon, while CO₂ concentration changes decrease isoprene emissions uniformly across the globe (relative trend: −0.20 % per year; absolute trend: −0.7 Tg yr⁻¹; p<0.01). For meteorological variability effects, the areas with statistically significant trends are smaller than those considering the effects of vegetation or CO₂ concentration changes and show different signs in different regions. For instance, the positive trends are mainly in Europe, Central Asia, and northern South America, while the negative trends are in Central Africa, central Australia, and the area around Lake Baikal.

In terms of the effects of individual factors, the trends in isoprene emissions illustrated in Fig. 8e and f are both statistically significant in most areas, which result from the changes in LAIv (Fig. 7e) and PFT (Fig. 7a–d, especially for tree cover). Note that changes in LAIv and PFT lead to opposite trends for the total global isoprene emissions (LAIv effect: 0.07 % per year, p<0.01; PFT effect: −0.1 % per year, p<0.01), resulting in a weak effect of vegetation changes on the global isoprene emission trend (−0.03 % per year; p=0.16). The trends in isoprene emissions illustrated in Fig. 8g–i result from the changes in 2 m temperature (Fig. 6c), surface solar radiation (Fig. 6d), and soil moisture (Fig. 6e). Temperature changes play a dominant role in the isoprene emission trends, with statistically significant positive trends in Europe and South America of up to 1 % to 2 % per year, while negative trends are found in Central Africa, northwestern South Asia, and southern North Asia, with maximum trends of −2 % to −5 % per year. Radiation changes exert a weaker effect on the isoprene emission trends (−0.2 % to −0.5 % per year) but show statistically significant negative trends over most regions (global scale; relative trend: −0.06 % per year; absolute trend: −0.2 Tg yr⁻¹; p<0.05). In some specific regions, changes in soil moisture significantly affect the trends in isoprene emissions. Positive trends are mainly distributed in Central Africa and Central Asia (2 % to 5 % per year), while negative trends are located in eastern Amazonia and Australia (up to −5 % to −10 %).

Figure 8 Spatial distribution of isoprene emission trends from 2001 to 2020. (a) ALL represents simulated results considering interannual variability of all drivers. (b, c, and d) VEG, MET, and CO₂ represent simulated trends when considering only interannual variability of vegetation, meteorology, and CO₂ concentrations, respectively. (e, f) Contributions from individual vegetation parameters including LAIv (leaf area index of vegetation-covered surfaces) and PFT (plant functional type). (g, h, and i) Contributions from individual meteorological factors including T2m (surface 2 m temperature), RAD (surface solar radiation), and SM (soil moisture). Stippling denotes regions where the trend is statistically significant (p<0.1) using the Mann–Kendall test. Trends are expressed by the relative change in percentage (i.e., the linear change during 2001–2020 divided by the mean value) and absolute change. p values represent statistically significant levels using the Mann–Kendall test.

To understand the discrepancies in isoprene emission trends across regions and to identify the contribution of individual influencing factors, nine hotspot regions (same as Fig. 4) are selected for further analysis. As shown in Fig. 9a, in some regions of the Northern Hemisphere, such as Europe, East Asia, and South Asia, both vegetation and meteorology changes strongly boost the isoprene emissions. In contrast, isoprene emissions in Central Africa decline sharply, mainly due to the meteorological changes. In other regions, especially Southeast Asia, South America, and Australia, the effects of the three factors offset each other, resulting in overall small trends.

The vegetation parameters can be decomposed into LAIv and PFT cover (Fig. 9b, Eqs. 2 and 3; the difference between EMIT_VEG and EMIT_VEG_FIX_PFT/EMIT_VEG_FIX_LAIv in Table 1). The LAIv calculated in MEGANv3.2 combines the information from LAI and VCF, which implies the growth of vegetation per unit of vegetation-covered area, while PFT cover reflects the amount of extension for different vegetation types on the ground surface. Large increases in the isoprene emissions mainly exist in Europe, North Asia, East Asia, South Asia, and Australia, where changes in LAIv and PFT jointly promote the emissions. These regions have experienced a dramatic increase in LAIv, along with an obvious expansion of tree cover. In particular, crop cover in South Asia and shrub cover in Australia increase rapidly (Fig. 9d). In other regions, however, changes in the PFT cover cancel out the contribution from increased LAIv, resulting in weaker positive trends in the isoprene emissions. Southeast Asia and South America even have experienced a moderate decrease of −0.14 % to −0.16 % per year in the emissions, mainly due to shifts in local vegetation functional types induced by the deforestation (Fig. 9d). In these regions, primary broadleaf evergreen forests are converted to some economic trees and crops (e.g., rubber, oil palm, and sugar cane).

Three meteorological variables namely temperature, radiation, and soil moisture, which are the main influencing factors from meteorology, are selected for quantifying their effects on the trends in isoprene emissions (Fig. 9c, d; the difference between EMIT_MET and EMIT_MET_FIX_T2m, EMIT_MET_FIX_RAD, and EMIT_MET_FIX_SM in Table 1). The emission trends caused by meteorological factors are large in Europe, East Asia, South Asia, South America, and Central Africa, although the emission trends are not statistically significant in some regions (p>0.1). Elevated temperature is found to dominate the rise in emissions in Europe and South America, and temperature cooling dominates the falling in emissions in Central Africa. The enhancement of isoprene emissions in East and South Asia is dominated by the increase in soil moisture. Overall, for meteorological parameters, changes in temperature and soil moisture exert the largest influence on the trends in isoprene emissions, with little effects from changes in radiation except in some regions of South America and Central Africa.

The spatial pattern of monoterpene emission trends is similar to that of isoprene emission trends but shows stronger positive trends in larger areas, especially in greening hotspots (Fig. S3). The discrepancies are mainly owing to the fact that monoterpene emissions are more sensitive to changes in LAIv (Fig. S3d; LAIv effect: 0.15 % per year, p<0.01) and not sensitive to the inhibition effect of CO₂. Monoterpene emission trends differ significantly in some regions from isoprene emission trends, which is mainly due to the different impacts of vegetation and meteorological factors on their emissions (Fig. S4). In North America, Southeast Asia, and South America, monoterpene emissions exhibit significantly increasing trends, while isoprene emission trends are not significant. For vegetation parameters, changes in LAIv increase the monoterpene emissions more significantly than the isoprene emissions, while changes in PFT exert a weaker effect on the monoterpene emission trends. The combined effects cancel out negative emission trends in Southeast Asia and South America and promote the increase in emissions in other regions. The effect of meteorological factors on monoterpene emission trends is similar to that on the isoprene emission trends, although the effect is weaker.

Figure 9 Trends in isoprene emissions and associated influencing factors from 2001–2020 in nine regions. (a) ALL represents simulated results considering interannual variability of all drivers, while VEG, MET, and CO₂ represent simulated results considering only interannual variability of vegetation, meteorology, and CO₂ concentrations, respectively. (b) Contributions from individual vegetation parameters including LAIv (leaf area index of vegetation-covered surfaces) and PFT (plant functional type). (c) Contributions from individual meteorological factors including T2m (surface 2 m temperature), RAD (surface solar radiation), and SM (soil moisture). (d) Trends in individual influencing factors including vegetation parameters, meteorological factors, and CO₂ concentrations. Tree, shrub, grass, and crop represent the cover for the four plant functional types (PFTs). The single, double, and triple asterisks denote 90 %, 95 %, and 99 % confidence intervals (CIs) using the Mann–Kendall test, respectively.

4 Discussion

In this study, we simulate the global BVOC emissions over the last 20 years using the latest version of the BVOC emission model MEGANv3.2 with time-varying vegetation and meteorological parameters and CO₂ concentrations. We quantify the impacts of the various driving factors on the BVOC emission trends. Isoprene emissions are found to increase significantly in Europe, East Asia, and South Asia, with comparable contributions from vegetation and meteorological factors. Warming in Europe and wetting of soil in East Asia and South Asia lead to an increase in the isoprene emissions, while cooling in Central Africa leads to a significant decrease in the isoprene emissions.

The results presented here demonstrate the heterogeneous spatial and temporal variabilities of vegetation and meteorological factors, leading to different trends in BVOC emissions in different regions of the world. In some regions (e.g., Europe, East Asia, and South Asia), vegetation and meteorological factors combine to promote the BVOC emissions, resulting in significant positive trends, while in other regions, such as South America and Southeast Asia, the effects of vegetation and meteorological factors on the BVOC emissions offset each other, resulting in weak emission trends. For vegetation parameters, the biomass-related parameter LAIv and type-related parameter PFT cover are both proportional to the intensity of BVOC emissions. Although multiple satellite data reveal a global trend towards greening (i.e., increased LAI), there is a substantial decrease in PFT cover (i.e., tree cover) in some local regions of Southeast Asia and South America, leading to a decreasing trend in BVOC emissions there.

Note that the selection of the reference year (i.e., the year 2001 in Table 1) may cause variations in simulated BVOC emissions, mostly affecting the magnitude rather than the sign of the absolute trends. Since this study focuses on the relative trends in BVOC emissions (i.e., the ratio of absolute trend to multi-year means), differences in the reference year have little effect on the magnitude and sign of our estimation results. We show that during 2001–2020, global isoprene emissions decrease by a trend of −0.07 % per year, with changes in the meteorological factors, vegetation, and CO₂ concentrations contributing to a trend of 0.15 % per year, −0.03 % per year, and −0.20 % per year, respectively. Some previous studies (e.g., Chen et al., 2018; Opacka et al., 2021; Sindelarova et al., 2022) also investigated the trends in isoprene emissions during 2000s and 2010s and obtained similar results regarding the positive trends due to meteorological changes and negative trends due to vegetation change. Chen et al. (2018) showed that the isoprene emission trend was −0.1 % per year between 2000 and 2015. Sindelarova et al. (2022) used time-varying ERA5 meteorological data and a static CLM4 land cover map to calculate the isoprene emission for 2000–2019 and suggested an isoprene emission trend of 0.35 % per year. When replacing the static land cover map with the annually varying ESA-CCI data, Sindelarova et al. (2022) also found that the emission trend decreases to 0.24 % per year. Opacka et al. (2021) found that the MEGAN-MOHYCAN model driven by the time-varying ERA-Interim meteorological field yielded an isoprene emission trend of 0.94 % per year for 2001–2016, while time-varying vegetation parameters can offset the positive trend by 0.04 % per year (based on MODIS land cover data) or 0.33 % per year (based on modified MODIS land cover data). The differences in the magnitude of trends due to meteorological changes or vegetation changes may be ascribed to the difference in the input data (meteorology and vegetation) and emission parameterizations. In fact, driven by the ERA5, ERA-Interim, and MERRA-2 meteorological fields, the model simulates a trend in isoprene emission from 0.15 % to 0.94 % per year (a factor of 6 difference). Different vegetation data lead to the simulated isoprene emission trends differing from −0.03 % to −0.33 % per year (a factor of 11 difference). In addition, different model parameterizations can lead to large discrepancies in emission changes. For example, based on the same vegetation data, Li et al. (2022) showed that isoprene emission induced by vegetation changes increased by about 20.6 % over China between 2001 to 2018, which is larger than that in our study (12.0 %). The difference can be explained by the fact that we use different vegetation emission factors and vegetation activity factors from Li et al. (2022).

Vegetation cover tends to increase in a large part of global land areas and decrease in smaller areas such as South America (i.e., the Amazon) and Southeast Asia (Fig. 6a). The trends in global isoprene emission due to vegetation cover change are largely determined by the decreasing trends in the Amazon and Southeast Asia due to the large contributions from these regions to global total isoprene emissions (Figs. 3a, 8f). The decrease in isoprene emissions in the Amazon and Southeast Asia is mainly due to the change in PFT, such as from the original tropical broadleaf evergreen forest to broadleaf deciduous trees with agricultural economic benefits (Figs. 7, 9). In some areas, such as Europe, East Asia, and South Asia, the increased forest coverage due to afforestation contributes significantly to the increase in isoprene emissions, which is similar to the results from Chen et al. (2018). However, Sindelarova et al. (2022) simulated little changes in isoprene emissions in Europe, mainly due to differences in the satellite PFT data used in their and our studies.

In this study, we use the latest version of MEGAN, which considers more refined factors compared to its predecessors. However, there are still large uncertainties in simulated BVOC emissions and their trends, such as in the treatment of drought and heat wave stress and the emission factors of certain tree species. Previous observations have shown that isoprene and monoterpene emissions are affected differently by drought severity (Brilli et al., 2007; Kaser et al., 2022; Otu-Larbi et al., 2020; Šimpraga et al., 2011). Isoprene emissions remain unchanged under mild drought but increase under moderate drought with increased leaf temperature due to changes in stomatal conductance (Kaser et al., 2022). However, under severe drought, isoprene emissions drop due to reduced substrate supply (Brilli et al., 2007). The effects of drought on monoterpenes are similar to those of isoprene (Lavoir et al., 2009; Ormeño et al., 2007). Therefore, these processes need to be better described in the model (Wang et al., 2022). In addition, more BVOC flux observations are urgently needed for model validations. In particular, the vegetation emission factors for tree species may be largely biased due to a scarcity of available observations and need to be further refined.

5 Conclusions

In this study, the time-varying meteorological data and satellite observations are used to drive the latest BVOC emission model MEGANv3.2 to simulate global BVOC emissions for the past 20 years. The contributions from different factors to the trends in BVOC emissions from global to regional scales are quantified.

Compared with site observations, the model can simulate isoprene emissions within a factor of 10 at most stations but systematically overestimates the monoterpene emissions. Compared to space-borne isoprene retrievals, the model can capture high isoprene emission regions such as South America and Central Africa. There are large seasonal variations in isoprene emissions, which are mainly determined by temperature and vegetation variations. The relative contribution of isoprene (monoterpene) emissions to BVOC emissions tends to decrease (increase) with latitudes, which is mainly ascribed to the meridional variations of PFT cover and corresponding emission factors.

Isoprene emissions increase significantly in Europe, East Asia, and South Asia (at rates of 0.37 % per year–0.66 % per year), with changes in both vegetation and meteorological factors contributing almost equally to the trends. For different meteorological factors, isoprene emission trends are mainly driven by the increase in temperature in Europe and by the increase in soil moisture in East and South Asia. In South America and Southeast Asia, shifts in PFT cover lead to a significant decrease in the BVOC emissions, which cancels out nearly half of the increasing trends induced by the changes in meteorological parameters. In addition, despite the increase in global mean temperature, there is a decrease in temperature in Central Africa, resulting in a significantly decreasing trend in isoprene emission in this region (−0.74 % per year). The dominant factors of monoterpene emission trends are similar to those of isoprene emissions, while monoterpene emissions show a stronger increasing trend or a weaker decreasing trend in most regions. In addition, monoterpene emissions are more sensitive to changes in LAIv, resulting in more pronounced increasing trends in greening hotspots.

Overall, our study highlights the significant BVOC emission trends both globally and regionally. More importantly, the results from this study clarify the contributions from different drivers and deepen our understanding of long-term BVOC emission trends at regional to global scales. Changes in BVOC emissions may have important impacts on ozone and atmospheric particle formation, which consequently impact the atmospheric chemistry, radiation, and climate. These interactions involving BVOCs will be investigated by using a coupled meteorology–chemistry with the BVOC emission model.