Anthropogenic drivers of 2013–2017 trends in summer surface ozone in China

Fig. 1 shows the 5-y average (2013–2017) values of the summer mean and maximum MDA8 ozone at the ensemble of sites operated by the China Ministry of Ecology and Environment. The Chinese National Ambient Air Quality Standard for MDA8 ozone is 160 µg m⁻³, corresponding to 82 ppbv at 298 K and 1,013 hPa. This standard is exceeded over much of eastern China. The highest concentrations are in the North China Plain, with values as high as 150 ppbv. Summer mean MDA8 ozone is also highest over the North China Plain, with values of 60–80 ppbv.

Fig. 2 shows the monthly mean MDA8 ozone trends for 2013–2017 in the four megacity clusters highlighted in Fig. 1: Beijing–Tianjin–Hebei (BTH), Yangtze River Delta (YRD), Pearl River Delta (PRD), and Sichuan Basin (SCB). These four megacity clusters are specific target areas in Chinese government plans to decrease air pollution (www.mee.gov.cn/hjzl/dqhj/cskqzlzkyb/). The trends are presented as the anomalies for individual summer months relative to their 2013–2017 means. Also shown is the meteorologically driven variability as described by a multiple linear regression (MLR) model considering a number of meteorological variables from the NASA Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) reanalysis (Methods and Table 1) (23). We use only the top three meteorological predictors for each region (indicated in Fig. 2) to avoid overfitting the data. These include temperature, surface winds, relative humidity, and also surface pressure for PRD. These variables are frequently observed to be correlated with ozone air quality (24) and can be viewed as general indicators of stagnation. Temperature also affects ozone through its control of biogenic VOC emissions and peroxyacetyl nitrate chemistry (25). The coefficients of determination (R²) for the MLR model in fitting the observed ozone anomalies range from 0.60 to 0.86 after removal of the residual linear trends (in black in Fig. 2).

The residual trends in Fig. 2 may be reasonably attributed to the effect of changing anthropogenic emissions. Fig. 3 shows the general trend of this MDA8 ozone residual across China for 2013–2017 after the meteorologically driven variability from the top three variables has been removed for each grid cell with the MLR model. Trends that are statistically significant above the 90% confidence level are marked with black dots. There is a general regional increase in eastern China between Shanghai (YRD) and Beijing (BTH). There are also patterns of decrease in southern and northeastern China away from the major population centers. The average trends for the focus megacity clusters are 3.1 ppbv a⁻¹ for BTH, 2.3 ppbv a⁻¹ for YRD, 0.56 ppbv a⁻¹ for PRD, and 1.6 ppbv a⁻¹ for SCB (SI Appendix, Table S1). The trend in BTH is larger than the earlier 2003–2015 trend of 1.1 ppbv a⁻¹ reported by ref. 14. PRD and SCB show increases even though they are in southern China, indicating some difference between urban centers and the broader region.

Chinese anthropogenic emissions estimated in the MEIC inventory decreased by 21% for NO_x and increased by 2% for VOCs over the 2013–2017 period (17). Emissions of PM_2.5 and its precursors are estimated to have also decreased including by 59% for SO₂ (17). Trends for the four megacity cluster regions are given in SI Appendix, Table S2. Observed average PM_2.5 levels in summer during 2013–2017 decreased by 41% for BTH, 36% for YRD, 12% for PRD, and 39% for SCB. Aerosol optical depth (AOD) decreased by 20% in eastern China (SI Appendix, Fig. S1).

We examined the effects of these changes in NO_x emissions, VOC emissions, and PM_2.5 levels using the nested-grid GEOS-Chem model version 11-02 over Asia (60°–150°E, 10°S–55°N) with a resolution of 0.5° × 0.625°. The GEOS-Chem model includes detailed ozone–NO_x–VOC–aerosol chemistry (26) and has been evaluated in previous studies simulating surface ozone in China (27–30). Our baseline simulation for 2013 is driven by MERRA-2 meteorological data with anthropogenic emissions from the MEIC inventory for China (17) and MIX inventory for other Asian countries (31). SI Appendix, Fig. S2, evaluates the simulation for 2017 with the mean summer MDA8 ozone observations for that year. Observed and simulated concentrations average 58.5 ± 15.4 and 63.0 ± 14.8 ppbv, respectively. Spatial correlation between simulated and observed ozone is high (correlation coefficient R = 0.89).

We then conducted sensitivity simulations with 2013–2017 changes taken together and separately in Chinese NO_x and VOC emissions (SI Appendix, Fig. S3), PM_2.5 affecting aerosol chemistry, and AOD affecting photolysis rates (SI Appendix, Fig. S1) (Methods). All simulations were performed for the same meteorological conditions of JJA 2013 after 1 mo of initialization. Detailed description of the model configuration and the sensitivity simulations is given in SI Appendix.

Fig. 4 shows the differences in MDA8 ozone resulting from these 2013–2017 anthropogenic changes. Changes in NO_x and VOC emissions (mainly due to decreased NO_x emissions; SI Appendix, Fig. S4) increase ozone in the urban areas of BTH, YRD, and PRD and in the broader urban region around Beijing, while decreasing ozone elsewhere, following expected patterns of VOC-limited and NO_x-limited conditions. Ozone production in urban areas is expected to be VOC-limited because NO_x concentrations are very high, but ozone production on a more regional scale in summer is expected to be NO_x-limited. The modeled ozone sensitivity is generally consistent with previous measurement-based, satellite-retrieved, and model inferences of NO_x- vs. VOC-limited conditions for ozone production in China (4, 12, 22).

However, we find that changes in PM_2.5 are more important than changes in NO_x or VOC emissions in driving ozone trends, particularly in the North China Plain, and this is mainly due to aerosol chemistry rather than photolysis (Fig. 4). The relevant aerosol chemistry involves reactive uptake of the gaseous precursors to ozone formation, as described in GEOS-Chem by first-order reactive uptake coefficients γ (32). This includes reactive uptake of the hydroperoxy radical (HO₂) with coefficient γ = 0.2 and conversion to H₂O or H₂O₂ (32–34) and reactive uptake of nitrogen oxides (NO₂, NO₃, and N₂O₅) with conversion to HNO₃ (32, 35). Uptake of HO₂ is by far the dominant effect (Fig. 4). It accounts in the model for most of the sink of hydrogen oxide radicals (HO_x ≡ OH + peroxy) in eastern China (SI Appendix, Fig. S5). This suppresses the HO₂ + NO reaction by which ozone is produced. The effect is particularly important in the North China Plain where PM_2.5 concentrations are highest.

The importance of aerosol chemistry as a sink for ozone precursors in China has been previously pointed out in model studies (21, 22), which found ozone decreases of 6–12 and 10–20 ppb, respectively, over eastern China as a result of this chemistry. Ref. 21 found the dominant effect to be the reactive uptake of nitrogen oxides, but we find that effect to be small in part because of VOC-limited conditions and in part because summertime conditions are not conducive to nighttime NO₃/N₂O₅ chemistry.

The HO₂ uptake coefficient γ = 0.2 used in our simulation is consistent with a large body of experimental and modeling literature (32). It is specifically consistent with laboratory measurements of HO₂ uptake by aerosol particles collected at two mountain sites in eastern China (34), which showed γ values averaging 0.23 ± 0.07 and 0.25 ± 0.09 at each site. Ref. 34 attributed this reactive uptake to aerosol-phase reactions of HO₂ with transition metal ions (TMI) and organics. In our standard simulation, we assume that the product of HO₂ uptake is H₂O, as for example through Cu/Fe TMI catalysis (33):However, if Fe(II) reacts with HO₂ instead, then the product becomes H₂O₂:We conducted a sensitivity simulation assuming the product to be H₂O₂ instead of H₂O, and this showed no significant difference in results because the recycling of HO_x radicals from H₂O₂ is inefficient (SI Appendix, Fig. S6).

Overall, the pattern of simulated 2013–2017 ozone trends from the combined changes in emissions and PM_2.5 (Fig. 4) is roughly consistent with the observed pattern of residual (presumed anthropogenic) trends in Fig. 3. The largest increases extend from Shanghai (YRD) to the North China Plain. Ozone decreases over most of southern China except in urban regions (as in PRD and SCB). There are some discrepancies between model and observed trends. The model underestimates the observed trend in BTH, possibly because the 50-km grid is too coarse to resolve strongly VOC-limited conditions in urban cores. Observations show ozone increases in western China, whereas the model suggests that emission controls should have produced decreases. Terrain is high in that region so that ozone has a large background component (30), and the increasing trend could reflect the more general trend of increasing background ozone at northern midlatitudes (36). Anthropogenic emissions in western China may also be underestimated (31). Observations show mixed trends in the eastern peninsula of Shandong province as well as decreases in northeastern China that are not captured by the model. Eastern Shandong may be difficult to model due to marine influence. For northeastern China, the model simulates an ozone increase because of the PM_2.5 decrease, but it may overestimate the low PM_2.5 concentrations in that region (SI Appendix, Fig. S2).

There is a pressing need to continue to decrease PM_2.5 levels in China because of the benefit for public health. Our finding that decreasing PM_2.5 causes an increase in ozone calls for decreasing NO_x and VOC emissions to overcome that effect. Model sensitivity simulations decreasing either NO_x or VOC emissions relative to 2017 levels show ozone benefits from both in the four megacity clusters (Fig. 5), consistent with ozone production being in the transitional regime between NO_x- and VOC-limited (12). The larger gains are from NO_x emission reductions as the chemistry becomes increasingly NO_x-limited, but VOC emission reductions are important to decrease ozone in urban cores (SI Appendix, Fig. S7). Gains from decreasing NO_x and VOC emissions are additive (37); thus, there is benefit in decreasing both.

In summary, we analyzed the factors driving 2013–2017 trends in summertime surface ozone pollution across China, taking advantage of the extensive network data available since 2013. We removed the effect of meteorological variability by using a multiple linear regression model fitting surface ozone to meteorological variables. The residual shows an increasing trend of 1–3 ppbv a⁻¹ in urban areas of eastern China that we attribute to changes in anthropogenic emissions. Decrease in anthropogenic NO_x emissions can increase ozone in urban areas where ozone production is expected to be VOC-limited. However, we find that a more important and pervasive factor for the increase in ozone in the North China Plain is the rapid decrease in PM_2.5, slowing down the reactive uptake of HO₂ radicals by aerosol particles and thus stimulating ozone production. Decreasing ozone in the future will require a combination of NO_x and VOC emission controls to overcome the effect of decreasing PM_2.5. There is a need to better understand HO₂ aerosol chemistry and its implications for ozone trends in China. Extending the observational record beyond the relatively short 5-y period will also provide more insights into the factors driving ozone trends in China.

All of the measurements, reanalysis data, and GEOS-Chem model code are openly available for download from the websites given below. The anthropogenic emission inventory is available from www.meicmodel.org, and for more information, please contact Q.Z. ([email protected]).

Hourly surface ozone concentrations for JJA 2013–2017 were obtained from the public website of the China Ministry of Ecology and Environment (MEE): beijingair.sinaapp.com/. The network had 450 monitoring stations in 2013 summer, growing to 1,500 stations by 2017 and including about 330 cities. We average the hourly data on the 0.5° latitude × 0.625° longitude MERRA-2 grid and compute daily MDA8 ozone on that grid. Trend analyses use all available data for a given year. Only using sites with 5-y records does not change the results. Most sites in the four focused megacity clusters were already operational in 2013.

Meteorological fields for 2013–2017 were obtained from the MERRA-2 reanalysis produced by the GEOS of the NASA Global Modeling and Assimilation Office (accessible online through https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2) (23). The MERRA-2 data have a spatial resolution of 0.5° × 0.625°. They match well with observed daily maximum temperature and relative humidity at Chinese weather stations (SI Appendix, Fig. S8) (38) and provide us with a full gridded ensemble of meteorological variables. We average them over either 24 h or daytime hours (8–20 local time), depending on the variable (Table 1). All data are normalized for use in the MLR model (see below) by subtracting their 2013–2017 mean for that day of the year and dividing by the standard deviation.

A number of previous studies have examined meteorological influences on ozone variability in China (4, 9, 39, 40), On the basis of these studies we considered the correlation of MDA8 ozone across China with a large number of candidate meteorological variables from the MERRA-2 archive (SI Appendix, Table S3 and Fig. S9). This led us to adopt nine variables as featuring the strongest correlations (Table 1). We applied a stepwise MLR model for each 0.5° × 0.625° grid cell:

where y is the normalized daily MDA8 ozone concentration and ($x_1$, …, $x_9$) are the nine meteorological variables. The interaction terms are up to second order. The regression coefficients ${\beta }_k$are determined by a stepwise method adding and deleting terms based on Akaike information criterion statistics to obtain the best model fit (41). Similar MLR models have been successfully applied to quantify the effect of meteorological variability on air pollutants in North America, Europe, and China (42–44).

We first apply the MLR model to identify the key meteorological variables driving the variability of daily surface ozone for each grid cell. Only the three locally dominant meteorological variables are regressed onto deseasonalized monthly MDA8 ozone to fit the effect of 2013–2017 meteorological variability on ozone within a 0.5° × 0.625° grid cell. This is done to avoid overfitting. We find that the dominant meteorological variables driving ozone variability are consistent across grid cells on a regional scale.

The ozone simulations use the nested-grid version of the GEOS-Chem chemical transport model with detailed oxidant–aerosol chemistry, driven by MERRA-2 assimilated meteorological data and with a horizontal resolution of 0.5° × 0.625° over East Asia (version 11-02; acmg.seas.harvard.edu/geos/). Anthropogenic emissions in China are from the MEIC inventory (see below). The base simulation is for the summer of 2013, and sensitivity simulations examine the effects of 2013–2017 changes in Chinese anthropogenic emissions, PM_2.5, and AOD, as described below. Additional sensitivity simulations isolate the effects of PM_2.5 and AOD changes on photolysis rates, NO_x aerosol chemistry, and HO₂ aerosol chemistry. Results presented in Fig. 4 are differences between the sensitivity simulations and the base simulation. Further details on the GEOS-Chem simulations are in SI Appendix.

The MEIC (www.meicmodel.org) is used to estimate China’s anthropogenic emissions and their trends from 2013 to 2017 (17, 31). MEIC is a widely used bottom-up emission inventory framework that follows a technology-based methodology to calculate emissions from more than 700 anthropogenic source types in China.

Observed PM_2.5 concentrations during 2013–2017 are from the same MEE observation network as ozone. Local changes in PM_2.5 concentrations from 2013 to 2017 affecting aerosol chemistry are applied as scaling factors to GEOS-Chem aerosol surface areas in the boundary layer below 1.3 km. AOD trends for 2013–2017 are from the monthly level 3 product of the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument aboard the Aqua satellite, reported at 550-nm wavelength with a resolution of 1° × 1° (https://ladsweb.modaps.eosdis.nasa.gov/). These trends in AOD are applied as scaling factors to simulated AOD in the GEOS-Chem calculation of photolysis rates (see details in SI Appendix, sections 1 and 2).