Substantial ozone enhancement over the North China Plain from increased biogenic emissions due to heat waves and land cover in summer 2017

In the summer of 2017, heavy ozone pollution swamped most of the North China Plain (NCP), with the maximum regional average of daily maximum 8 h ozone concentration (MDA8) reaching almost 120 ppbv. In light of the continuing reduction of anthropogenic emissions in China, the underlying mechanisms for the occurrences of these regional extreme ozone episodes are elucidated from two perspectives: meteorology and biogenic emissions. The significant positive correlation between MDA8 ozone and temperature, which is amplified during heat waves concomitant with stagnant air and no precipitation, supports the crucial role of meteorology in driving high ozone concentrations. We also find that biogenic emissions are enhanced due to factors previously not considered. During the heavy ozone pollution episodes in June 2017, biogenic emissions driven by high vapor pressure deficit (VPD), land cover change and urban landscape yield an extra mean MDA8 ozone of 3.08, 2.79 and 4.74 ppbv, respectively, over the NCP, which together contribute as much to MDA8 ozone as biogenic emissions simulated using the land cover of 2003 and ignoring VPD and urban landscape. In Beijing, the biogenic emission increase due to urban landscape has a comparable effect on MDA8 ozone to the combined effect of high VPD and land cover change between 2003 and 2016. In light of the large effect of urban landscape on biogenic emission and the subsequent ozone formation, the types of trees may be cautiously selected to take into account of the biogenic volatile organic compound (BVOC) emission during the afforestation of cities. This study highlights the vital contributions of heat waves, land cover change and urbanization to the occurrence of extreme ozone episodes, with significant implications for ozone pollution control in a future when heat wave frequency and intensity are projected to increase under global warming.

1 Introduction

In recent decades, China has been facing severe air pollution issues, particularly for the winter PM_2.5 and summer ozone (Zheng et al., 2015; Cheng et al., 2016; Zhao et al., 2016). It has been noted that the mean concentration of PM_2.5 has generally decreased in the past few years, but the concentration of ozone shows an increasing trend (Li et al., 2017b, 2019; Wang et al., 2017; Chen et al., 2018a), suggesting a greater urgency for ozone pollution control. For instance, Li et al. (2017b) revealed an increase in annual mean ozone in 2016 of 11 µg m⁻³ compared to 2014 in China. Lu et al. (2018) found a 3.7 %–6.2 % increase per year in the mean ozone concentration over 74 cities in China from 2013 to 2017. Since ozone is harmful to both human health (Soriano et al., 2017) and vegetation (Emberson et al., 2009; Avnery et al., 2011), it is vital to investigate the possible mechanisms related to high ozone concentrations. Based on ozone observations from 2013 to 2017, the North China Plain (NCP, an area about 400 000 km² in size, with Beijing located on its northeast edge; 35–42^∘ N, 112–119^∘ E), is identified as the area with the most severe ozone pollution in China compared to other regions such as the Yangtze River Delta and Pearl River Delta, possibly linked to the stimulation effect from enhanced hydroperoxy radicals (HO₂) due to a reduction in aerosol sink resulting from the decrease in PM_2.5 during this period (Li et al., 2019). Chen et al. (2019) investigated the impact of meteorological factors such as temperature, wind speed and solar radiation on ozone pollution from 2006 to 2016 and noted that the severe ozone events in June 2017 around Beijing stand out and suggested a possible connection with the abnormal meteorological conditions. These studies motivated a need for a better understanding of the high-ozone problem over the NCP.

Tropospheric ozone is closely related to both anthropogenic emissions and biogenic emissions, including volatile organic compounds (VOCs) and nitrogen oxides (NO_x) (Sillman, 1995, 1999; Tonnesen and Dennis, 2000; Xing et al., 2011; Fu et al., 2012). In the past few years (i.e., 2012–2017), anthropogenic emissions such as NO_X continued to decrease (Liu et al., 2016) and anthropogenic VOCs changed little (Zhao et al., 2018; Zheng et al., 2018; Li et al., 2019). Biogenic VOC (BVOC) was reported to enhance hourly ozone by 3–5 ppbv in the NCP, especially in areas north of Beijing, based on a 2 d simulation from 31 July to 1 August 1999 (Wang et al., 2008). The annual BVOC emission in this area increased by 1 %–1.5 % per year from 1979 to 2012 (Stavrakou et al., 2014) due to changes in land use and climate. Broadleaf trees in general have a higher emission rate of BVOC than grass, shrubs and crops (Guenther et al., 2012). A dramatic increase in forest (trees) coverage is evident in the last 20 years over the NCP (Chen et al., 2018b), partly attributable to the “Three-north Forest Protection Project”. For example, trees planted before the 2008 Olympic Games doubled the BVOC emissions in Beijing from 2005 to 2010 (Ghirardo et al., 2016). Urban landscape may even emit more BVOC than natural forest because of favorable conditions such as lower tree densities and better light illumination (Ren et al., 2017). Ren et al. (2017) found that BVOC emitted by urban landscape accounted for 15 % of total BVOC emissions in Beijing in 2015. Over highly polluted urban areas of the NCP, ozone production is highly sensitive to VOC emissions (Liu et al., 2012; Han et al., 2018). Therefore, elevated BVOC emissions can greatly enhance ozone formation in the NCP.

Besides emissions, tropospheric ozone is also closely related to meteorological conditions, such as heat waves (Gao et al., 2013; Fiore et al., 2015; Otero et al., 2016), low wind speed and stagnant weather (Jacob and Winner, 2009; Sun et al., 2017; Zhang et al., 2018). Weather conditions concomitant with heat waves, including high temperature, low wind speed and little cloud coverage, may enhance ozone production (Jaffe and Zhang, 2017; Pu et al., 2017; Sun et al., 2019). At the same time, such meteorological conditions also promote emissions of BVOC and ozone formation (Zhang and Wang, 2016). Using a global model, Fu and Liao (2014) suggested a slight to moderate increase in biogenic isoprene west and north of Beijing due to land cover and land use alone and an even more obvious increase when meteorological changes are considered. In the summer of 2017, heat waves swept over the majority of the NCP, providing an excellent opportunity to investigate how the heat wave may have modulated BVOC emissions and subsequent severe ozone events in the NCP. Observation data and modeling are used to delineate various factors contributing to enhanced biogenic emissions and elevated ozone concentrations. More details of the data and model are provided in the “Methods” section.

2 Methods

Data and model configuration

The distribution of observed data is shown in Fig. 1. For instance, the meteorological observations used in this study such as daily maximum temperature, daily mean wind speed and daily total precipitation were obtained from the China Meteorological Data Service Center (CMA, http://data.cma.cn, last access: 23 September 2019), shown as blue dots in Fig. 1. Observed surface ozone data are obtained from the China National Environmental Monitoring Centre (http://www.pm25.in, last access: 23 September 2019), shown as red dots in Fig. 1. Meteorological Assimilation Data Ingest System (MADIS) hourly 2 m temperature, specific humidity, 10 m wind speed and direction are available from MADIS (https://madis.ncep.noaa.gov, last access: 23 September 2019), shown as hexagons in Fig. 1.

Figure 1 Distribution of observational sites over the NCP. Blue dots: daily maximum temperature daily mean wind speed at 10 m and daily total precipitation from the China Meteorological Administration (CMA); red dots: ozone monitoring sites from the China National Environmental Monitoring Centre; black hexagon: hourly temperature at 2 m (T2), specific humidity at 2 m (Q2), wind speed (WS10) and direction (WD10) at 10 m from MADIS; green box: urban area of Beijing. B, H and T represent Beijing, Hebei Province and Tianjin, respectively.

For modeling the meteorological conditions, WRF V3.8.1 is used in this study. The domain is centered at 34^∘ N, 110^∘ E, with a total of 34 vertical layers and top pressure at 50 hPa. The spatial resolution is 36 km. The physics parameterizations used in this study are the same as in our previous studies (Gao et al., 2017; Zhang et al., 2019), including the Morrison double moment microphysics (Morrison et al., 2009), the Rapid Radiative Transfer Model for General Circulation Models (RRTMG) longwave and shortwave radiation (Iacono et al., 2008; Morcrette et al., 2008), the unified Noah land surface model (Chen and Dudhia, 2001), the Mellor–Yamada–Janjic planetary boundary layer (PBL) scheme (Janjić, 1990, 1994; Mellor and Yamada, 1982), and the Grell–Freitas cumulus scheme (Grell and Freitas, 2014). The initial and boundary conditions were generated from the NCEP Climate Forecast System Reanalysis (CFSR) version 2 (Saha et al., 2013), with a spatial resolution of $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$.

For modeling atmospheric chemistry, the widely used Community Multi-scale Air Quality (CMAQ) model (Byun and Ching, 1999; Byun and Schere, 2006), with the latest version 5.2, was used in this study. The major gas phase chemistry was represented by Carbon Bond version 6 (CB6) and Aerosol Module Version 6 (AERO6) aerosol module. Initial and boundary conditions were from the Model for Ozone and Related chemical Tracers, version 4 (MOZART-4) (Emmons et al., 2010). A dynamical downscaling tool was developed in this study to link the Mozart output to CMAQ, based upon the package of Mozart to WRF-Chem (mozbc: https://www2.acom.ucar.edu/wrf-chem/wrf-chem-tools-community; last access: 23 September 2019). With this tool, the default clean air profile provided by the CMAQ 5.2 package was replaced by more realistic boundary variations at both the surface and different vertical levels. A continuous run from 1 June to 4 July was performed, with the first week discarded as spinup.

The anthropogenic emissions of air pollutants in China were estimated by Tsinghua University, detailed in previous studies (Wang et al., 2014; Zhao et al., 2013, 2017, 2018) and updated based on the Multiresolution Emission Inventory for China (MEIC, $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$; http://www.meicmodel.org/, last access: 23 September 2019) (Li et al., 2017a).

The biogenic emissions were calculated by the Model of Emissions of Gases and Aerosols from Nature version 2.1 (MEGAN; Guenther et al., 2006, 2012). MEGAN input data includes three components: plant functional type (PFT), leaf area index (LAI) and emission factors (EFs). There is a total of 19 emission species including isoprene and terpenes derived from more than 100 emission compounds. For each of the 19 species, the emission rates F_i (µg m⁻² h⁻¹) for a certain grid are defined in Eq. (1), with i denoting the species.

$$\begin{array}{}\text{(1)}& F_i={\mathit{\gamma }}_i\sum {\mathit{\epsilon }}_{i,j}{\mathit{\chi }}_j,\end{array}$$

where ε_i,j and χ_j are the emission factor and fractional coverage of plant functional type (j) in each grid, respectively. γ_i is the emission activity defined based on light (denoted as L), temperature (T), leaf age (LA), soil moisture (SM), leaf area index (LAI) and CO₂ inhibition (denoted as CI), following Eq. (2).

$$\begin{array}{}\text{(2)}& {\mathit{\gamma }}_i=C_{\mathrm{CE}}\mathrm{LAI}{\mathit{\gamma }}_{L,i}{\mathit{\gamma }}_{T,i}{\mathit{\gamma }}_{\mathrm{LA},i}{\mathit{\gamma }}_{\mathrm{SM},i}{\mathit{\gamma }}_{\mathrm{CI},i},\end{array}$$

where C_CE is the canopy environment coefficient and 0.57 was used following Guenther et al. (2012).

Compared with the previous version 2.0 with only four PFTs, there are 16 types of PFTs represented in the new MEGAN version (Guenther et al., 2006, 2012), allowing for more accurate estimations of PFT-differentiated emission factors. PFT and LAI data were from the MODIS MCD12Q1 (Friedl et al., 2010) and MCD15A2H datasets (Myneni et al., 2015), respectively. The eight vegetation types in MODIS were apportioned to the 16 PFT types in MEGAN 2.1 based on the temperature zone. For example, MODIS has only one type of broadleaf deciduous trees, while MEGAN 2.1 has three: broadleaf deciduous tropical, temperate and boreal trees. The broadleaf deciduous trees in MODIS are mapped onto the three MEGAN types based on the latitudinal boundaries of the tropical, temperate and boreal zones, with detailed mapping information provided in Table S4 in the Supplement. Monthly mean LAIs were used in this study. The meteorological conditions used to generate biogenic emission in MEGAN were provided by the WRF simulation.

3 Results

3.1 Observed ozone features

The Technical Regulation on Ambient Air Quality Index (HJ633-2012) defines six classes of ozone-related pollution based on the daily maximum 8 h ozone concentration (MDA8). Classes I and II are clean conditions (MDA8 less than 82 ppbv), class III (82–110 ppbv) indicates slight pollution, class IV (110–135 ppbv) represents medium pollution, and classes V and VI are severe pollution conditions with MDA8 higher than 135 ppbv. Utilizing the observed MDA8 from the China National Environmental Monitoring Centre (http://www.pm25.in, last access: 23 September 2019), we first analyze the severe ozone pollution events, considering their large impact on human health. The observed MDA8 was interpolated to a $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid. Figure 2 shows the number of severe ozone pollution days (MDA8 greater than 110 ppbv) during the summer of 2014–2017. The number of severe ozone pollution days in 2017 is larger than 9 in most areas, which is substantially higher than that of the other 3 years when most areas have fewer than 6 d. Frequent occurrence of severe ozone pollution happens in southern Beijing and south of Hebei Province (the area marked with letter H in Fig. 1).

Figure 2 The number of severe ozone pollution days (MDA8 greater than 110 ppbv) during the summer of 2014–2017 over the NCP.

3.2 Meteorological factors modulating the high-ozone events

Correlation between MDA8 ozone and daily maximum 2 m temperature (Fig. 3) shows statistically significant values for all 4 years, confirming the significant impact of temperature on ozone. However, the correlation in 2017 is obviously higher than in the other 3 years, and the regression slope of 4.21 ppbv ^∘C⁻¹ is about 1.07 to 1.84 ppbv ^∘C⁻¹ higher than in the other 3 years, demonstrating the larger impact of temperature in 2017. Both the higher correlation (0.74) and the larger slope in 2017 are contributed mainly by days with ozone above the top 10 % (104 ppbv), which are related to the long-lasting high-ozone periods (see Table S1 and Fig. S1 in the Supplement) during 14–21 June and 26 June–3 July. Removing data above the top 10 % brings the correlation (0.63) and slope closer to those of the other 3 years (Fig. S2). Furthermore, the mean temperature in 2017 is not statistically different from that of the other 3 years, suggesting that the higher temperature period has disproportionate effects on ozone. Jaffe and Zhang (2017) also found a larger regression slope between ozone and temperature during the abnormally warm month of June 2015 in the western US compared to the previous 5 years with more normal temperatures.

Figure 3 The correlation between summer MDA8 ozone and daily maximum 2 m temperature (T_max) for 2014–2017 over the NCP. Regional mean was calculated from the observational sites over the NCP, so each data point corresponds to a regional mean value of MDA8.

To further delve into the meteorological factors modulating the ozone variations in the summer of 2014–2017, the time series of summer MDA8 ozone is shown in Fig. 4, along with daily maximum temperature, wind speed and daily total precipitation. From Fig. 4d, the two long-lasting ozone episodic events (event 1: 14–21 June; event 2: 26 June–3 July) occur during heat waves concomitant with stagnant (calm or low wind speed), dry (little or no precipitation) air and strong solar radiation (not shown), conducive to ozone formation and accumulation. During the first 3 d of these two high-ozone episodic events, the regional mean daily maximum temperature is 32.3 ^∘C, accounting for the 90th percentile relative to a 30-year period during 1987–2016. Moreover, almost half of the stations with at least 3 continuous days exceeding their respective 95th percentile from 1987 to 2016, satisfying the definitions of heat waves. This feature during the heat wave period was illuminated in Table S2 as well, showing that among all the observational stations with MDA8 ozone exceeding 110 ppbv, 87 % (62 %) and 96 % (81 %) occur with daily precipitation less than 1 mm (daily precipitation less than 1 mm and daily mean wind speed lower than 3 m s⁻¹). Long-lasting hot and stagnant weather conditions were not clearly observed during 2014–2016 (Fig. 4a–c).

Figure 4 Time series of observed MDA8 ozone (red lines; based on sites from the China National Environmental Monitoring Centre; red points in Fig. 1), daily maximum temperature at 2 m (blue lines), daily mean wind speed at 10 m (green lines) and daily total precipitation (yellow bars) over the NCP (based on sites from CMA; blue dots in Fig. 1) during the summer from 2014 to 2017. The regional precipitation was set to zero for a certain day if less than 15 % (9 sites) of the total sites (58 sites) with daily total precipitation greater than 1 mm. The horizontal blue dashed lines in each panel denote 31 ^∘C.

Figure 5 The correlation between summer MDA8 ozone and daily maximum VPD during 2014–2017 over the NCP. Regional mean was calculated from the observational sites over the NCP, so each data point corresponds to a regional mean value of MDA8.

3.3 Effect of land use and biogenic emission on ozone

Biogenic emissions contribute importantly to ozone formation. The MEGAN model has been widely used to simulate biogenic emissions in air quality modeling studies (Guenther et al., 2012), but recent research suggested that biogenic emissions may be underestimated in the model for several reasons:

• a.

  water-stressed impact on biogenic emissions. Zhang and Wang (2016) found that two high-ozone events in the US were associated with excess isoprene release due to dry and hot weather conditions that induced water stress in plants. The increased vapor pressure deficit (VPD; the pressure difference between saturation vapor and ambient vapor) drives the release of more isoprene, but the VPD effect on biogenic emissions has not been taken into consideration in MEGAN 2.1, so the subsequent influence of biogenic emissions on ozone may be largely underestimated. Zhang and Wang (2016) suggested a doubling of daily biogenic isoprene when the daily VPD reaches 1.7 kPa or greater. It should be noted that this parameterization was based upon the observed information over the US; more tests may be needed in future when applying it to areas besides the US. The monthly mean VPD spatial distribution in June 2017 (Fig. S3) as well as the high correlation between observed MDA8 ozone and VPD (Fig. 5; with time series shown in Fig. S4) suggest enhanced isoprene emission in the NCP, so we will test this VPD mechanism using model simulations. Please note that in the latest version MEGAN 3 (Jiang et al., 2018), a new approach was developed to quantify the drought effect on the isoprene emissions based on both photosynthesis and water stress, yielding a general reduction of monthly mean isoprene emission across the globe, including northern China. The impact of changes in isoprene emissions, based on the new method, on ozone formation deserves further evaluation in future.

• b.

  Effect of changes in land cover on biogenic emissions. As reflected by the much higher emission factor, biogenic isoprene emission is enhanced in broadleaf forest relative to other land cover types such as needleleaf forest, shrub, grass or crop (Table 2 in Guenther et al., 2012). In the NCP, broadleaf tree is the dominant land cover type, and its coverage has been increasing dramatically since the 1970s, primarily as a result of the Three-north Protection Forest System project. For example, based on Moderate Resolution Imagine Spectroradiometer (MODIS) land use data (Friedl et al., 2010), the coverage of broadleaf deciduous temperate trees nearly doubled from 2003 to 2016 over the NCP (Fig. 6a–c). This has resulted in a substantial increase in isoprene emissions between 2003 and 2016 (Fig. 6), particularly north of Beijing, Hebei and Tianjin, where the increase is more than 200 %. Combining point (a) described above, the underestimation of biogenic emission due to changes in land cover may be exaggerated in years with high temperatures and high VPD. It is vital to quantify the effect of land cover changes on biogenic emissions such as isoprene and the subsequent impact on ozone formation.

• c.

  Impact of urban landscape on biogenic emission. The land use type cataloged in the MODIS MCD12Q1 product (Friedl et al., 2010) does not take into consideration urban green spaces, which may lead to a 15 % underestimation of total BVOC emissions in 2015 over Beijing (Ren et al., 2017). Generally, urban ozone production is highly sensitive to VOC emissions (Xing et al., 2011; Liu et al., 2012). Bell and Ellis (2004) found a doubling of ozone in urban areas relative to rural areas for the same percentage increase in biogenic emissions. The impact of biogenic emission from urban landscape on urban ozone formation has not been considered in previous studies. For a sensitivity analysis, we added a 15 % increase in the total BVOC emissions in Beijing to investigate their impact on urban ozone formation. These emissions were distributed evenly in the urban core area of Beijing as the increase in biogenic emissions from urban landscape were only available for Beijing.

Figure 6 Spatial distribution of broadleaf deciduous trees in 2003 (a), 2016 (b) and their differences (2016–2003; c) and the biogenic isoprene emissions during the heat wave periods (14–21 June; 26 June–3 July 2017) based on the land cover in 2003 (d), 2016 (e) and their differences (2016–2003; f).

Figure 7 The comparison of isoprene concentrations between model simulations and observations in Beijing. The black dots represent the observed data from various contributions to the literature, whereas the hollow triangles (in black, red, green and blue) represent the model simulations for the four cases described above (cases 2–5). For each observational dataset, the corresponding reference number is given after the site name in (a) and (b), with site locations shown in (c). One exception is the unpublished work in THU* which is from the observations using a proton-transfer reaction time-of-flight mass spectrometer (PTR-ToF-MS) conducted by Tsinghua University (paper in preparation). Please note that no observation period matches our simulation time exactly, making the comparison more qualitative rather than quantitative. However, the model evaluation did match the respective location and time (i.e., daytime or selected hour) among different observations. The model simulation period used in the comparison is from 8 June to 4 July 2017. For observations, in (a), the dots represent the mean isoprene concentrations during daytime in August from 2005 to 2011 at Peking University (PKU; Zhang et al., 2014; left side of a) and from 16 July to 18 August 2008 at the Chinese Research Academy of Environmental Science (CRAES; Yang et al., 2018; right side of a). In (b), the dots on the left represent the mean isoprene concentration from 08:00 and 16:00 (local standard time) in August from 2004 to 2006 (with detailed measurement time shown in Table 1 of Shao et al., 2009) in PKU. The observational data on the right side of (b) are on a daily mean scale during a certain period (with one site of CY showing minimal and maximal daily mean values during the period) from four sources. The two leftmost dots are located at the campus of Tsinghua University (THU), with one from 15 to 20 August 2006 (Duan et al., 2008) and the other from 14 July to 5 August 2017 (paper in preparation as explained above). The third dot represents data measured at PKU from 24 July to 27 August 2008 (Liu et al., 2015), and the fourth dot indicates data observed at Chaoyang District (CY; Gu et al., 2019).

To elucidate the mechanism modulating the ozone events discussed above, the regional meteorology and air quality model WRF/CMAQ was used to conduct simulations from 8 June to 4 July 2017. The WRF simulations generally meet the benchmark standard for meteorological variables (Table S3). For air quality simulations, five scenarios were designed, with biogenic emissions ignored in the base case. Compared to the base case, case 2 adds biogenic emission associated with the land cover of 2003 and cases 3, 4 and 5 are the same as case 2 except for the inclusion of the VPD effect, both VPD and land cover of 2016, and VPD and land cover of 2016 combined with the effect of urban green spaces, respectively. To validate the reasonableness of adding the biogenic emission, we first evaluate the simulated isoprene concentration, one of the most important species closely related to ozone formation, from WRF/CMAQ among different cases. Since there is a lack of observed ambient isoprene concentration during this study period, the data available (mostly over Beijing) from the literature were retrieved and used as cross comparison with the model results (Fig. 7). From Fig. 7a, b, the observed mean isoprene concentration ranges from 0.4 to 1.6 ppbv at various sites of Beijing. By taking into consideration isoprene emissions from VPD, land cover of 2016 and urban green spaces (case 5) the model simulations yield the best performance, with isoprene concentration of 0.8 ppbv to 1.4 ppbv. However, the other cases (with isoprene concentrations of 0.1 ppbv to 0.2 ppbv) substantially underestimate the isoprene concentrations. Therefore, the isoprene emissions from urban green spaces (comparing case 5 and case 4) in Beijing play a vital role in the isoprene concentrations, which subsequently affect the ozone formation which will be further evaluated and discussed below.

Since the effect of urban landscape was only applied to Beijing in case 5, we use case 4 (combination of VPD and land cover change effects) (referred to as B_MDA8) as the reference. Therefore, we first compare MDA8 ozone in case 4 with observations. To facilitate the comparison, observational data were interpolated to the model grids, and reasonable performance is achieved with an MFB/MFE (mean fractional bias and mean fractional error) of −7 %/16 % (Fig. 8). Considering the mean bias likely attributed to the factors such as emission uncertainty or model inherent biases, a bias correction was applied to each case by adding 7 % of mean observed MDA8 ozone from 8 June to 4 July 2017.

Figure 8 MDA8 ozone evaluation over the NCP from 8 June to 4 July in 2017. NMB, NME, MFB and MFE represent normalized mean bias, normalized mean error, mean fractional bias and mean fractional error, respectively.

Zooming into the two ozone episodic events (14–21 June, 26 June–3 July), the mean MDA8 values of case 4 are 98.02, 108.89, 95.75 and 98.98 ppbv for the NCP, Beijing, Hebei and Tianjin, respectively, during the heat wave periods (14–21 June 2017; 26 June–3 July 2017), whereas the MDA8 ozone value for the case (case 1) without biogenic emission are 87.15, 93.06, 84.78 and 89.65 ppbv for the corresponding region. The ozone increment from case 2 to case 5 (as well as observations; magenta stars in Fig. 9a) relative to case 1 is shown in Fig. 9a for these regions. Including biogenic emission based on the land cover of 2003 (case 2) yields an extra mean MDA8 ozone of 7.84 ppbv (8 % of B_MDA8), 9.96 ppbv (9 % of B_MDA8), 7.86 ppbv (8 % of B_MDA8) and 6.99 ppbv (7 % of B_MDA8) for the NCP, Beijing, Hebei and Tianjin, respectively (yellow bars in Fig. 9a), compared to case 1. Including the VPD effect (case 3) adds an extra mean MDA8 of 1.71 ppbv in the NCP compared to case 2, and the enhancement is highest in Beijing (3.08 ppbv) (green bars in Fig. 9a). Additional MDA8 ozone enhancement is simulated by including the effect of land cover change (increase in natural broadleaf forest; Fig. 6a–c; case 4), i.e., an extra MDA8 of 1.32 ppbv in the NCP relative to case 3, with the highest contribution of 2.79 ppbv in Beijing (blue bars in Fig. 9a). The urban landscape (case 5) in Beijing yields an extra 4.74 ppbv or 4 % of MDA8 compared to case 4, almost doubling the effect of VPD and land cover change in Beijing. The larger percentage increase in MDA8 ozone (41 % from Fig. 9a, which is shown in Fig. 9b as well) due to urban landscape relative to the prescribed 15 % increase in BVOC emission in Beijing supports the notion of an amplified MDA8 ozone response in urban areas because of the high sensitivity of ozone to VOC emissions, which well matches observational data (magenta star).

Figure 9 Biogenic contribution to MDA8 ozone during the heat wave periods (14–21 June; 26 June–3 July), shown by the individual (a) and percentage contribution (b) of standard biogenic emissions (Bio-emis) using MEGAN 2.1 with the land cover of 2003 (Bio-emis), VPD effect, land cover (LC) change and urban green spaces. The color bars (a) represent the simulated contributions of biogenic emissions (yellow), VPD (green), land use changes (blue) and urban green (red) to the MDA8 ozone concentrations in the NCP, Beijing, Hebei and Tianjin, respectively. The magenta stars in (a) represent the observed biogenic emissions calculated by subtracting the contribution to MDA8 ozone simulated in the base case from the observed total MDA8 ozone. The box-and-whisker plot shows the contribution of biogenic emissions, VPD, land cover change and urban green spaces to the total MDA8 ozone in BTH (black) and Beijing (brown) (y axis on the left), and the percentage increment (right y axis) of VPD, land cover change, and urban green relative to MDA8 induced by Bio-emis for BTH (blue) and Beijing (purple). Please note that urban green spaces are only available for Beijing. The top and bottom edges of the boxes represent the 75 and 25 percentiles, with the centered line and red dot showing the median and mean, respectively.

Figure 10 A schematic diagram of the impact of biogenic emission on ozone formation. N-BVOC refers to natural biogenic emission, P-BVOC refers to the biogenic emission from planted forest and in this study represents the increase in forest coverage. U-BVOC refers to urban BVOC generated from urban green spaces. The red thick upward arrows indicate that extra VOCs may be induced by the heat waves.

To further illustrate the contributions of BVOC to MDA8, Fig. 9b shows the contribution of biogenic emissions (Bio_emis, based on land cover of 2003), VPD, land cover change and urban landscape (or urban green) to MDA8 as a fraction of the MDA8 of B_MDA8 (left y axis in Fig. 9b) and as a percentage increment relative to the MDA8 contributed by biogenic emissions in case 2 (right y axis in Fig. 9b) in BTH (Beijing, Tianjin, Hebei; with letters B, T and H marked in Fig. 1) and Beijing. For BTH, the mean contribution to B_MDA8 is 9 %, 2 % and 2 % for Bio_emis, VPD and land cover change (red dots in the black bars in Fig. 9b), respectively, with maximum contributions of 22 %, 10 % and 10 %. For Beijing, the contributions of Bio_emis, VPD, land cover change and urban landscape are 9 %, 3 %, 3 % and 4 %, respectively (red dots in the brown bars in Fig. 9b). Urban landscape (19 %) contributes more than Bio_emis (17 %) in the urban area of Beijing in terms of the maximal contribution (maximum value of the brown box in Fig. 9b). Compared with Bio_emis, the mean increments are 19 % and 17 % for VPD and land cover change (red dots in the blue bars in Fig. 9b). For Beijing, the mean additional enhancements are 30 %, 28 % and 41 % for VPD, land cover change and urban landscape relative to Bio_emis (red dots in the purple bars in Fig. 9b), with a combined increment of 99 % compared to the MDA8 ozone contributed by biogenic emission based on the land cover of 2003. Although only grid cells with both simulations and observations available are used in Fig. 9b, the results are similar if all model grids points are used (not shown).

In order to demonstrate whether changes in land cover and VPD play any role during normal ozone conditions, we conducted another sets of simulations (the same as cases 2–4 discussed above) from 8 June to mid-July in 2016, a similar period as in 2017. The mean MDA8 ozone concentrations over the NCP during this entire period in 2017 for case 2 is 79.03 ppbv, and statistical significant enhancement (1.34 ppbv) was achieved in case 3. In comparison to case 3, the land cover change in case 4 shows a statistically significant increase as well (1.13 ppbv). As expected, looking at the entire period in 2016 (8 June–4 July), a statistically significant (and even higher in comparison to 2016) increase was achieved in case 3 (1.55 ppbv) compared to case 2 (90.11 ppbv) and case 4 (1.23 ppbv) compared to case 3. Therefore, the land cover and VPD may be applied in both episodic events and conditions with normal ozone concentrations.

Herein the mechanisms for ozone enhancement are summarized in the schematic of Fig. 10. Both natural and anthropogenic emissions contribute to ozone formation. Because of the Three-north Protection Forest System project, natural forest north of Beijing has more than tripled in area coverage compared to 2003, leading to an increasing trend in biogenic emissions. Under heat wave conditions, biogenic emissions may be further enhanced through the effect of VPD in addition to the effect of temperature. For urban areas, even more biogenic emissions may be emitted from urban landscapes. All these mechanisms for increasing biogenic emissions could enhance ozone formation, particularly over urban areas such as Beijing.

4 Discussion

The mechanisms contributing to the severe ozone pollution events in the summer of 2017 in the NCP were investigated. Two severe tropospheric ozone pollution events occurred in the NCP during the periods of 14 to 21 June and 26 June to 3 July. We provided support for the roles of the observed meteorological conditions including high temperature and stagnant dry weather, which favor high ozone concentrations. More importantly, the influence of biogenic emissions on ozone formation was investigated in more detail by incorporating important biogenic emission factors that are typically ignored in regional model simulations. Biogenic emissions based on the land cover of 2003 yields an extra mean MDA8 ozone of 7.84 ppbv for the NCP. Including the VPD effect and land cover change adds 1.71 and 1.32 ppbv of ozone in the NCP. These contributions are even larger in Beijing, with VPD adding 3.08 ppbv and land cover change adding 2.79 ppbv. Most notably, biogenic emissions from urban landscape (i.e., green spaces) have so far not been considered in ozone regional modeling studies to our knowledge. By adding this source in the urban area of Beijing, substantial ozone enhancement was simulated, bringing the WRF/CMAQ simulation of MDA8 closer to observations. The urban landscape in Beijing yields an extra 4.74 ppbv of MDA8, comparable to the combined effect of VPD and land cover change in Beijing. Together, the combined effect of VPD, land cover change and urban landscape doubles the effect of biogenic emission calculated based on the land cover of 2003 and not including the VPD and urban landscape effects. Please note that although the urban isoprene emission from landscape in Beijing only accounts for 15 % (Ren et al., 2017), the location of the emissions may play a much larger role in contributing to the urban isoprene concentration. As was shown in Fig. 6, most of the isoprene emissions from the forest in Beijing are located in the rural area, which is relatively far from the urban area. Considering the short lifetime of isoprene, it may not be as efficient as the urban isoprene emission resulting from urban landscape directly in modulating the isoprene concentrations. Therefore, the urban isoprene emission may play a much more significant role in urban photochemical reactions compared to the isoprene emissions from the forest over the rural areas.

The BVOC emissions from urban green spaces are projected to increase by more than twice in 2050 due to urban area expansion (Ren et al., 2017). Together with the more frequent heat waves projected for the future (Gao et al., 2012; Zhang et al., 2018), the impact of biogenic emissions on ozone pollution in the NCP will likely play an increasingly important role in ozone pollution and should be taken into considerations in future air quality management plans to address issues of air quality and health. The effect of urban green spaces was only considered in Beijing in this study as we lack the data to parameterize this effect in other regions. Considering the substantial effect of urban green spaces on urban ozone formation, it is vital to evaluate similar effects in other cities where ozone pollution is a concern.