Rapid increase in summer surface ozone over the North China Plain during 2013–2019: a side effect of particulate matter reduction control?

While the elevated ambient levels of particulate matters with aerodynamic diameter of 2.5 µm or less (PM_2.5) are alleviated largely with the implementation of effective emission control measures, an opposite trend with a rapid increase has been seen in surface ozone (O₃) in the North China Plain (NCP) region over the past several years. It is critical to determine the real culprit causing such a large increase in surface O₃. In this study, 7-year surface observations and satellite retrieval data are analyzed to determine the long-term change in surface O₃ as well as driving factors. Results indicate that anthropogenic emission control strategies and changes in aerosol concentrations as well as aerosol optical properties such as single-scattering albedo (SSA) are the most important factors driving such a large increase in surface O₃. Numerical simulations with the National Center for Atmospheric Research (NCAR) Master Mechanism (MM) model suggest that reduction of O₃ precursor emissions and aerosol radiative effect accounted for 45 % and 23 % of the total change in surface O₃ in summertime during 2013–2019, respectively. Planetary boundary layer (PBL) height with an increase of 0.21 km and surface air temperature with an increase of 2.1 ^∘C contributed 18 % and 12 % to the total change in surface O₃, respectively. The combined effect of these factors was responsible for the rest of the change. Decrease in SSA or strengthened absorption property of aerosols may offset the impact of aerosol optical depth (AOD) reduction on surface O₃ substantially. While the MM model enables quantification of an individual factor's percentage contributions, it requires further refinement with aerosol chemistry included in the future investigation. The study indicates an important role of aerosol radiative effect in development of more effective emission control strategies on reduction of ambient levels of O₃ as well as alleviation of national air quality standard exceedance events.

1 Introduction

Elevated ambient levels of ozone (O₃) are of great concern due to their important impact on human health, ecosystem productivity, atmospheric chemistry and climate change (Monks et al., 2015; Tai et al., 2014; Tan et al., 2019). O₃ is produced by a series of photochemical reactions involving nitrogen oxides (NO_x= NO + NO₂) and volatile organic compounds (VOCs) in the presence of solar radiation. Ambient levels of O₃ are highly dependent on emissions of O₃ precursors, solar radiation and other physical processes such as regional and vertical transport (Sun et al., 2019; Ni et al., 2018; Liu et al., 2019; Wang et al., 2016b). While O₃ concentrations show a steady decreasing trend in Europe and North America, an opposite trend with an accelerating increase rate is observed in China (Lu et al., 2018; Li et al., 2019a). Due to high nonlinearity of the O₃–NO_x–VOC relationship and complexity of processes governing ambient levels of O₃, a large uncertainty remains in the determination of the impact of different driving factors on changes in surface O₃ concentrations under the polluted atmospheric conditions. Thus, accurate quantification of relative contributions of individual factors to the large increase in surface O₃ concentrations over the heavily polluted regions such as China continues to represent one of the major challenges to research communities and government policy makers.

Anthropogenic emissions are the key in driving change in surface O₃. With rapid development of industrialization and urbanization, anthropogenic emissions of NO_x and VOCs, two major precursors of O₃ formation, have been increasing significantly in China over the past several decades (Zeng et al., 2019). For instance, tropospheric columns of NO₂ (TCNO₂), an indicator of anthropogenic emission intensity of NO_x were increased by 307 % in Beijing from 1996–2011 (Huang et al., 2013), which caused a strong increasing trend of O₃ in the lower troposphere. Meanwhile, an increase in surface O₃ at a rate of 2 % a⁻¹ was observed in Beijing from 1995–2005 (Ding et al., 2008), and a similar increase with 1–2 ppb a⁻¹ was monitored at urban and remote sites in eastern China (Sun et al., 2016; Gao et al., 2017; Ma et al., 2016; Tang et al., 2009). However, little light was shed on change in surface O₃ compared to its counterpart PM_2.5 which was elevated to the severe pollution level in eastern China, especially over the North China Plain (NCP) region (Zeng et al., 2019; Zhai et al., 2019). The severity of PM_2.5 pollution has been largely alleviated after the stringent emission control strategies were implemented by the Chinese government at a national level in 2013 (Zeng et al., 2019). According to the estimate by Multi-resolution Emission Inventory in China (MEIC), anthropogenic emissions of PM_2.5 decreased by approximately 60 % and NO_x emissions decreased by 21 %. Significant reductions were also seen in other air pollutants such as SO₂ but not for VOCs, which showed an increase of 2 % instead over the period of 2013–2017 (Zheng et al., 2018). As a result, monthly mean PM_2.5 concentrations decreased by 41 % for the Beijing–Tianjin–Hebei (BTH) region, which is similar to the NCP region presented in this study, and aerosol optical depth (AOD) was reduced by 20 % in eastern China (Li et al., 2019a). However, an opposite trend with an accelerating increase rate of O₃ was observed in the NCP region during this period (Lu et al., 2018; Cooper et al., 2014). The fact that O₃ formation was dominated by a VOC-sensitive regime may partly account for such an increase in the NCP region, but it is not clear how much the change of surface O₃ is attributed to anthropogenic emission control efforts.

Aerosol radiative effect is another factor imposing a large constraint on change in surface O₃. Aerosols attenuate surface-reaching solar near-ultraviolet (UV) radiation effectively and reduce photolysis rate of NO₂, a key parameter in determining O₃ formation. The impact of aerosol radiative effect on photolysis rate of NO₂ or O₃ photochemical production is highly dependent on aerosol optical properties as described by AOD, single-scattering albedo (SSA) and asymmetry factor. AOD is a measure of extinction of solar beam by aerosols (e.g., dust and haze), used as a proxy of representing severity of fine particulate matter pollution or aerosol mass concentrations. SSA denotes the relative contributions of scattering versus absorption effect to total aerosol extinction efficiency with “0” for pure absorption and “1” for pure scattering effect. Both numerical simulations and observations showed that aerosols with UV-scattering effect may accelerate photochemical production of O₃ but aerosols with a strong absorption property (e.g., mineral dust and soot) may inhibit O₃ production in the atmospheric boundary layer (Dickerson et al., 1997; Mok et al., 2016). The lowest photolysis rate coefficient was closely linked with the highest AOD (Liu et al., 2019; Dickerson et al., 1997). It was observed that surface PM_2.5 concentrations decreased by 41 % whereas surface O₃ increased at a rate of 3.1 ppb a⁻¹ over the BTH region from 2013–2017 (Li et al., 2019a). Decrease in PM_2.5 was considered to be one of the important causes leading to such an increase in surface O₃ due to additional O₃ production associated with reduced sink of hydroperoxy radicals (HO₂) (Li et al., 2019a). They pointed out that increase in surface O₃ associated with decrease in PM_2.5 was more prominent than that with reduction of NO_x emissions over the NCP region where O₃ formation was dominated by a VOC-limited regime. Liu and Wang (2020a, b) found the reduction of PM emissions increased the O₃ levels by enhancing the photolysis rates and reducing heterogeneous uptake of reactive gases (mainly HO₂ and O₃), of which the latter is more important than the former. A similar impact associated with aerosol radiative properties on O₃ production was observed in other regions over the world. For instance, the combined effect associated with optical properties of BrC and black carbon (BC) reduced the net change in O₃ production by up to 18 % compared to BC alone in the Amazon Basin (Mok et al., 2016). Thus, surface O₃ changes are dependent on not only aerosol concentrations (AOD used as a proxy) but also aerosol optical properties such as SSA. Relative importance of different aerosol optical property parameters to change in surface O₃ needs to be addressed.

Figure 1 (a) Locations of North China Plain (NCP, 32–40^∘ N, 114–121^∘ E) and three other major air pollution regions, Beijing–Tianjin–Hebei (BTH, 37–41^∘ N, 114–118^∘ E), Yangtze River Delta (YRD, 30–33^∘ N, 118–122^∘ E) and Pearl River Delta (PRD, 21.5–24^∘ N, 112–115.5^∘ E) in China with terrain heights included and (b) the location of an ultraviolet (UV) radiation observational site, YCA (Yucheng site), areas covered by the NCP region and surrounded by mountains.

NCP, the largest alluvial plain of China, is surrounded by the Yan Mountains with a main peak of 2116 m in the north, the Taihang Mountains with the highest elevation of 2882 m in the west, and the Dabie and Tianmu mountains in the south, bordered by the Yellow Sea in the east (see Fig. 1). Such a complex terrain is not conducive to dispersion and dilution of air pollutants and makes them easily trapped. Meanwhile, the total energy consumption was increased by more than a factor of 5 from 1985–2016 (Zeng et al., 2019). The NCP has become one of the most polluted regions in China. The highest PM_2.5 concentration reached 900 µg m⁻³ during winter, and heavy PM_2.5 pollution events were the major concern to air quality during that period (Gu, 2013; An et al., 2019), but surface PM_2.5 concentrations have reduced substantially. Meanwhile, O₃ exceedance events became more frequent and more serious in the NCP region (Zhang et al., 2015; Lang et al., 2017; Zhai et al., 2019). Hourly surface O₃ concentrations went up to 150.0 ppb and the increase rate reached 3.1 ppb a⁻¹, much higher than those observed in other polluted regions such as the Yangtze River Delta (YRD) and Pearl River Delta (PRD) in China (Li et al., 2019a; Lyu et al., 2019). The elevated surface O₃ has become an emerging critical air quality issue in this region (Wang et al., 2006; Shi et al., 2015). Understanding of the driving factors such as a rise in surface O₃ has become a very hot topic (e.g., Li et al., 2019a, b). However, most of the related studies are limited to the contributions of atmospheric chemistry and changes in O₃ precursors' emissions. Relative importance of aerosol radiative effect associated with a substantial decrease in aerosols or PM_2.5 and meteorological variability to the enhancement of surface O₃ is not well qualified.

In this study, seven-year air quality observational data provided by the China National Environmental Monitoring Center (CNEMC) network are examined to determine the temporal and spatial variations in surface O₃ over the NCP region over the period of 2013–2019. A series of analyses are presented to investigate the long-term change trend of surface O₃ and the statistical relationships with NO_x and VOC emissions, meteorological variables, and aerosol radiative optical property parameters. A box model with Master Mechanism (MM) is then utilized to determine the response of surface O₃ to the key driving factors. The specific objectives include (1) to identify the key factors driving the increase in surface O₃ over the NCP, the most polluted region in China, and (2) to quantify the relative contributions of anthropogenic emissions (e.g., NO_x and VOCs), aerosol concentrations, aerosol optical properties and meteorological variability to the increase in surface O₃ in summertime during 2013–2019.

2 Data and methods

2.1 Observational data

Data used in this study include hourly-averaged surface observations of O₃ and PM_2.5 from 2013 to 2019 provided by the CNEMC (http://106.37.208.233:20035/, last access: 20 December 2020). UV data measured at the Yucheng site (i.e., YCA, 36.87^∘ N, 116.57^∘ E) in the NCP region are obtained from the Chinese Ecosystem Research Network (http://www.cern.ac.cn/, last access: 20 December 2020) from the years 2013–2016. AOD is derived from the monthly level 3 product of the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument aboard the Aqua satellite, reported at the 550 nm wavelength with a resolution of 1^∘ × 1^∘ (Platnick, 2015). TCNO₂ data are retrieved from the daily level 3 products of the Ozone Monitoring Instrument (OMI) aboard the Aura satellite with a resolution of 0.25^∘ × 0.25^∘ (Krotkov, 2019). Shortwave radiation data are provided by the Land Data Assimilation System (FLDAS) (NASA, 2018) at a resolution of 0.1^∘ × 0.1^∘. SSA retrieved from OMI/Aura Near UV Aerosol Optical Depth and Single Scattering Albedo V003 (OMAERUV) (Torres, 2006) at 388 nm are used to evaluate the impact of aerosol scatting/absorption properties on change in surface O₃. Daily max temperature at 2 m (T_2 max), 10 m wind speed and the planetary boundary layer height (PBLH) are derived from the Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) reanalysis data at a horizontal resolution of 0.5^∘ × 0.625^∘ (Global Modeling and Assimilation Office, 2015).

2.2 Model description and configurations

The MM model is utilized to quantify the relative contributions of anthropogenic emissions and aerosol optical and radiative properties to the change in surface O₃. The MM is a chemistry box model, originally developed and updated by the scientists at the National Center for Atmospheric Research (NCAR). It includes a detailed and flexible gas phase chemical mechanism consisting of approximately 5000 reactions for simulating temporal variations in chemical species of interest. The hydrocarbon chemistry in the MM is treated explicitly with photo-oxidation of partly oxygenated organic species included. Alkanes, alkenes and aromatics are considered to be initial hydrocarbon reagents in the gas-phase mechanism. The Gear-type solver is used in the MM model to handle large numbers of chemical reactions and species and the integration time steps varied during the simulations (Madronich and Calvert, 1989, 1990). The TUV (Tropospheric Ultraviolet and Visible radiation) model is called by the MM model for update of chemical reaction rates every 15 min. This model computes time-dependent chemical evolution of an air parcel initialized with a known composition and additional emissions. It is assumed that no dilution is included in the simulations given the difficulty of getting inputs to calculate the dilution rate. The transport in and out of air pollutants reached a quasi-equilibrium state over the study domain, and no heterogeneous processes were included in the MM model. The MM model has been widely used to investigate the impact of different factors such as emissions, chemistry and meteorological conditions on simulations of O₃ and other chemical species (e.g., Liu et al., 2019; Geng et al., 2007).

Photolysis rate j(NO₂) is calculated by using the Tropospheric Ultraviolet and Visible (TUV) radiation model which is embedded into the NCAR MM (Madronich and Flocke, 1999). In the fully coupled system, the TUV is called by the MM model for update of photolysis rates of NO₂ and other chemical species (e.g., H₂O₂, O₃, NO₃, N₂O₅) every 15 min dynamically. The TUV model is initialized with the monthly means of AOD, SSA and total columns of O₃ retrieved from satellite measurements as well as other meteorological parameters such as cloud fractions at the central point of the NCP (36^∘ N, 117.5^∘ E) in June.

HO₂ radicals are important to O₃ formation. HONO photolysis as the primary production of OH radicals and formaldehyde (HCHO) photolysis as the net radical source of HO₂ can lead to major changes in the HO_x and NO_x budget that may have an important effect on O₃ production and loss (e.g., Aumont et al., 2003; Brasseur and Solomon, 2006; Lin et al., 2012). The role of HO₂ radicals can be determined by the following reactions.

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\mathrm{HO}}_{\mathrm{2}}+\mathrm{NO}\to {\mathrm{NO}}_{\mathrm{2}}+\mathrm{OH},\text{(R2)}& {\displaystyle }{\mathrm{NO}}_{\mathrm{2}}+hv\to \mathrm{NO}+\mathrm{O}\left({}^{\mathrm{3}}\mathrm{P}\right),\text{(R3)}& {\displaystyle }{\mathrm{O}}_{\mathrm{2}}+\mathrm{O}\left({}^{\mathrm{3}}\mathrm{P}\right)\to {\mathrm{O}}_{\mathrm{3}},\end{array}$$

where hv represents ultraviolet radiation at the wavelengths of 200–400 nm. The MM model has a capability of quantifying the role of radicals in O₃ formations under different pollution conditions.

The MM simulations are conducted for the predefined box as shown in Fig. 1 to represent ensemble mean behaviors and responses of the model to changes of different model inputs over the NCP region. The 24 h simulations are conducted with the initial hour at 00:00 local time (LT). The inputs of the simulations include meteorological data (e.g., air temperature, cloud and PBLH), aerosol radiative properties (i.e., AOD and SSA) and emissions. While all the meteorological inputs are generated from observational data, the initial values of chemical species such as VOC species (e.g., acrylic, 2-methylbutane, toluene, p-xylene, isoprene), N₂, O₂, H₂O, NO₂, O₃, etc., are obtained from climatology or background values.

Emissions (NO_x and VOCs) are calculated from the MEIC emission inventory. Aerosol radiative property parameters from MODIS and OMAERUV are assumed as constants for all the simulations. All the simulations are driven by the monthly means averaged over the entire NCP region. The temporal variations at an interval of 4 h are included in the model inputs to represent the diurnal variations in different meteorological variables such as T_2 max and PBLH from MERRA-2 reanalysis.

Six groups with a total of 16 numerical experiments with the MM model are designed to quantify the roles of different factors in driving change in O₃ concentrations (Table 1). Case A stands for the base case in which the emissions were generated from MEIC in the base year 2012 (http://www.meicmodel.org/, last access: 20 December 2020) with an adjustment for 2013 use, and the spatial distributions of NO_x and VOCs are presented in Fig. S1. Case B represents a scenario for the year 2019 with NO_x and VOC emission changes by −35 % and +10 % with respect to the case in the year 2013, respectively. The changes in NO_x emissions (−35 %) and VOC emissions (+10 %) in 2019 were obtained by extrapolating their respective changes during the period from 2013 to 2017 (Li et al., 2019a). Case C1 denotes a scenario with a decrease in AOD from 1.0 (i.e., the case for the year 2013) to 0.75 (i.e., for the year 2019) according to MODIS measurements, and six other members in group C are used to examine the impact of varying AOD on the change in surface O₃. Case D1 is the one with a change of SSA from 0.95 (2013) to 0.93 (2019). Case E is a scenario with T_2 max increase from 29.9 ^∘C in 2013 to 32.0 ^∘C in 2019 based on the regional average calculated with the MERRA-2 reanalysis in the NCP region. Case F is designed to assess the impact of the increased PBLH (i.e., increase from 0.76 km in 2013 to 0.97 km in 2019) on surface O₃ change in the NCP region. Case G mimics the year 2019 representing the combined effect of changes in emissions, AOD, SSA, T_2 max and PBLH. The 24 h simulations are completed for each case to quantify the contributions of individual factors to the changes in surface O₃ from 2013 to 2019. More details of the numerical experiments are presented in Table 1.

Table 1 A summary of numerical experiments with the NCAR MM model.

^a The year 2013: NO_x emission is 2.0×10¹² mol cm⁻² s⁻¹, and VOC emission is 7.3×10⁹ mol cm⁻² s⁻¹. ^b The year 2019: NO_x emission is 1.3×10¹² mol cm⁻² s⁻¹, and VOC emission is 8.0×10⁹ mol cm⁻² s⁻¹.

3 Results and discussion

3.1 Spatiotemporal variations in surface O₃, PM_2.5, AOD and TCNO₂

Figure 2 shows a comparison of spatially distributed monthly means of the maximum daily hour average (MDA8) O₃, 24 h average PM_2.5, AOD and TCNO₂ over eastern China in June between 2013 and 2019 derived from in situ and satellite observations. It is clear that NCP was the most polluted region with the highest values of MDA8 O₃, PM_2.5, AOD and TCNO₂ over the past decade. The 24 h average PM_2.5 concentrations were higher than 75.0 µg m⁻³ (the Grade II National Ambient Air Quality Standard, NAAQS defined for residential areas) at most of the monitoring stations across the NCP region in June 2013. The highest 24 h average PM_2.5 reached to 766.0 µg m⁻³ and the corresponding AOD was 1.0. Compared to the well-established monitoring network of PM_2.5, observational sites for O₃ were pretty sparse except for the BTH, YRD and PRD across eastern China in 2013. While the TCNO₂ was 2 times higher than that observed in North America (Stavrakou et al., 2008), exceedance events of the MDA8 O₃ were not frequently observed across eastern China in 2013. PM_2.5 was the major air pollutant in the NCP region during that time period.

PM_2.5 concentrations, AOD and TCNO₂ have been reduced substantially as a result of the implementation of a strict anthropogenic emission reduction policy in 2013. For instance, the monthly mean of PM_2.5 concentrations decreased from 95.5 to 33.2 µg m⁻³ with a percentage reduction of 65 %. Monthly mean AOD was reduced from 1.0 in 2013 to 0.75 in 2019, indicating that PM_2.5 continued to decrease at a rate of −10 % a⁻¹ to −11 % a⁻¹ which was similar to that during 2013–2017 (Li et al., 2019a). A similar decrease trend was seen in both TNO₂ (Fig. 2g, h) and in situ NO₂ measurements (Fig. S2). On the other hand, a rapid increase in surface O₃ concentrations was observed in the NCP region over the past several years. The hot spots with MDA8 O₃ higher than 75.0 ppb were extended to the entire NCP as well as the neighboring regions in 2019 (Fig. 2b). The highest MDA8 O₃ reached 112.8 ppb in 2018, which was even higher than the level (110.0 ppb) observed in Los Angeles (Lin et al., 2017). Compared to the cases observed in 2017 (Li et al., 2019a, b), air pollution events with higher surface O₃ became more severe and more frequent. The frequency of NAAQS exceedance events for surface MDA8 O₃ (i.e., greater than 160 µg m⁻³) in June increased from 30 % in 2013 to 63 % in 2019. Here percentage represents the proportion of MDA8 exceedance days to a total of 30 d (i.e., June).

Figure 2 A comparison of spatial distributions of monthly mean of MDA8 O₃ (ppb) (a, b) and PM_2.5 (µg m⁻³) (c, d) obtained from in situ observations, AOD (e, f) and tropospheric column of NO₂ (TCNO₂, 10¹⁶ cm⁻²) (g, h) derived from satellite observations between 2013 (in left column) and 2019 (in right column) in eastern China (NCP indicated by the box).

Reduction in NO_x emissions and a slight increase in VOC emissions could be part of the reasons causing such an increase over the NCP region where O₃ formation was dominated by a VOC-limited regime. To better understand the relationship of increase in surface O₃ with the decrease in NO₂, the change in monthly mean O_x (a sum of O₃ and NO₂) was plotted in Fig. S3. It is clear that O_x showed an increasing trend over the past 7 years during daytime and nighttime in both urban Beijing and the NCP region. Meanwhile, Li et al. (2019a) attributed the increase to the fact that removal of HO₂ radicals was reduced and more O₃ production was promoted. On the other hand, attenuation of UV radiation became less evident as PM_2.5 or AOD continually decreased. Strengthening UV radiation may accelerate photolysis of NO₂ and eventually lead to more O₃ production. Importance of aerosol radiative effect in the increase in surface O₃ via acceleration of photolysis of NO₂ can be further evaluated through numerical experiments.

Figure 3 A comparison of spatial distributions of monthly mean of T_2 max (^∘, a, b), wind speed (m s⁻¹, c, d), surface-reaching shortwave radiation (W m⁻², e, f) and PBLH (m, g, h) between 2013 (in left column) and 2019 (in right column) in eastern China (the NCP indicated by the box).

Meteorological conditions are another critical factor affecting O₃ production. Typically, higher air temperature is responsible for higher photochemical reaction rates and more O₃ photochemical production (Porter and Heald, 2019). As shown in Fig. 3, the NCP was the hottest spot region with T_2 max which was about 4.0 ^∘C higher than that in the neighboring regions. In addition, the increase rate of T_2 max in the NCP was larger than that observed in other regions in eastern China. T_2 max and surface-reaching shortwave radiation increased by 3 % and 7 %, respectively, over the past several years. In addition to human factors such as urbanization and industrialization, the decrease in aerosols (e.g., PM_2.5 and AOD) could be an important factor driving such a rise in air temperature due to weakening aerosol radiative effect.

3.2 Yearly changes in surface O₃ during 2013–2019 and driving factors

As presented above, NCP was the most polluted region with extremely high ambient levels of air pollutants. Surface O₃ showed a rapid increase over the period from 2013–2019 while PM_2.5 and other pollutants such as NO_x experienced a significant reduction. O₃ has become a major air quality concern in summer. June was the month with the highest monthly mean MDA8 O₃ concentrations (Fig. S4). In this section, we attempt to investigate the yearly change rate and to identify the factors that drove such a large increase in surface O₃ over the NCP region throughout the period of 2013–2019.

Figure 4 shows the yearly changes in monthly means of MDA8 O₃, PM_2.5, AOD, SSA, TCNO₂, T_2 max and PBLH over the NCP region in June from 2013–2019. The change in monthly mean of surface MDA8 O₃ showed an opposite trend to that of PM_2.5 concentrations and other air pollutants. A similar large change trend was seen in diurnal variation patterns (Fig. S5). The increase rate of monthly mean MDA8 O₃ (4.6 ppb a⁻¹) during 2013–2019 was much higher than that observed in the same region during the period of 2005–2015 (1.1 ppb a⁻¹) (Ma et al., 2016) and other regions such as Mount Tai, the YRD, Hong Kong and North America where the changes were less than 2.1 ppb a⁻¹ during the same time period (e.g., Sun et al., 2016; Gao et al., 2017; Wang et al., 2017; Xu et al., 2019). At the same time, a large decrease can be found from the time series of PM_2.5, AOD, TCNO₂ (Fig. 4b–d) and in situ NO₂ measurements (Fig. S6). It is noted that SSA also showed a decreasing trend (Fig. 4e). Decrease in SSA was likely due to the fact that reduction of inorganic aerosols (e.g., sulfate and nitrate) was larger than that of carbonaceous ones (Zhang et al., 2020). Another noticed feature is that MDA8 O₃ showed a decreasing trend in 2019 relative to 2018, which was opposite to that during 2013–2018 (Fig. 4a). It is worth paying attention to the change trend in the coming years.

Figure 4 Long-term changes in monthly mean of (a) MDA8 O₃, (b) PM_2.5, (c) AOD, (d) TCNO₂, (e) SSA, (f) T_2 max and (g) PBLH averaged over the North China Plain in June over the period of 2013–2019.

To understand the factors driving the change in surface O₃, a series of scatter plots are presented to examine the relationships between the surface MAD8 O₃ and individual factors such as aerosol optical properties (i.e., AOD and SSA), TCNO₂, T_2 max and surface-reaching shortwave radiation over the past 7 years in June (Fig. 5). The values discussed here represent the monthly means. MAD8 O₃ showed two different regimes with an opposite dependence of O₃ formation on PM_2.5 concentrations. The first regime showed a decrease trend with increasing surface PM_2.5 when PM_2.5 concentrations were less than approximately 140.0 µg m⁻³ whereas the second one showed no trend with increasing PM_2.5 when PM_2.5 concentrations were higher than 140.0 µg m⁻³. The first regime was highly related to aerosol radiative effect, which has been discussed above. For the second regime, the impact of aerosol radiative effect on surface O₃ photochemical production seemed very minor or even negligible. Instead, O₃ production was suppressed significantly and MAD8 O₃ concentrations were less than 20.0 ppb. In this case, removal of surface O₃ through titration of NO was not effective and surface O₃ showed an increase rather than a decrease trend with increasing NO_x concentrations under the strong NO_x conditions as indicated by TCNO₂ higher than 40–45×10¹⁵ (cm⁻²) in the troposphere. Here the threshold value of 140.0 µg m⁻³ represents an observed reality in this region but it needs to investigate whether such a threshold value exists in other regions.

Figure 5d–e further demonstrate the critical role of meteorological factors in change of surface O₃. MDA8 O₃ showed a near-linear increasing trend with increasing T_2 max and surface-reaching shortwave radiation with respective linear regression correlation coefficients of 0.88 and 0.93. Increase in T_2 max and strengthening shortwave radiation caused by decrease in PM_2.5 (a proxy of aerosols) played a positive role in driving the increase in surface O₃ in the NCP region. On the other hand, MDA8 O₃ showed a decreasing trend with 10 m wind speed (Fig. 5f). That may explain why improvement of stagnation atmospheric conditions may alleviate severity of surface O₃ pollution to some extent. The positive correlation between the PBLH and O₃ shown in Fig. 5g represents one case when radiation is stronger and temperature is higher, which is favorable for O₃ formation. Meanwhile, higher PBLH could enhance the transport down of O₃-enriched air aloft, resulting in an increase in surface O₃ (Reddy et al., 2012). On the other hand, some studies found a negative correlation between the PBLH and O₃. They claimed that a shallower PBL may suppress the dispersion of pollutants and lead to higher O₃ (Yan et al., 2018; Jiang et al., 2015; Wei et al., 2016; Huang et al., 2005)

Enhancement of UV radiation resulting from reduction in surface PM_2.5 represents an important mechanism in driving the increase in surface O₃ concentrations. It can be further illustrated by Fig. 6. While UV radiation displays a nonlinear decreasing trend with surface PM_2.5 concentrations, surface O₃ (hourly) shows a near-linear increasing trend with surface-reaching UV radiation. UV radiation attenuation approaches a constant with a value of 0.1–0.3 MJ m⁻² when surface PM_2.5 concentrations reach around 300 µg m⁻³ or above.

Analyses presented above demonstrate that all the exceedance events of MDA8 are observed under conditions with PM_2.5 less than 60 µg m⁻³, TCNO₂ of equal to or less than 5.0×10¹⁵ (cm⁻²), T_2 max higher than 28.0 ^∘C and surface-reaching shortwave radiation stronger than 250.0 W m⁻². Reduction in aerosols (e.g., surface PM_2.5 as a proxy) concentrations may strengthen UV radiation, increase T_2 max and eventually promote more surface O₃ production.

Figure 5 Response of MDA8 O₃ to (a) PM_2.5, (b) TCNO₂, (c) AOD, (d) T_2 max, (e) shortwave radiation, (f) wind speed and (g) PBLH observed in the NCP region, China, during 2013–2019.

3.3 Relative contributions of different driving factors to increase in surface O₃

In this section, the box model MM is utilized to quantify the relative contributions of individual driving factors to the increase in surface O₃ over the NCP region during 2013–2019. A simulation–observation comparison is presented to evaluate the performance of the MM model on simulations of surface O₃ (Fig. 7), of which the O₃ observations averaged over all the stations in NCP are considered to be the standard observed concentrations. The simulated O₃ peak was about 1 h later than the observation, which was likely due to uncertainty of emission inventory and other meteorological factors. Overall, the MM model was able to mimic the observed variation pattern and peak value as indicated by the correlation coefficient of 0.95 between simulated and observed O₃.

Figure 6 (a) The relationships of surface O₃ concentrations (hourly) with (a) UV radiation and (b) UV radiation with PM_2.5 concentrations based on the observations at the Yucheng site during the time period of 08:00–17:00 LT in June 2013–2016.

Figure 7 Comparisons of (a) regional averaged surface O₃ observations in the NCP and simulated surface O₃ (A, control case) and (b) simulated net changes in O₃ among different driving-factor conditions.

A series of numerical experiments were then completed with the MM model to quantify the relative contributions of anthropogenic emissions (i.e., NO_x and VOCs), AOD, SSA, air temperature and PBLH to the change in surface O₃ over the NCP region during 2013–2019. The results are presented in Table 2. The changes in emissions of O₃ precursors (i.e., NO_x and VOC_s) (i.e., Case B) and decrease in AOD (i.e., Case C1) were the two major contributors with their respective positive contributions of 45 % and 70 % to the increment in surface O₃. But increase in surface O₃ associated with AOD reduction was largely offset by the reduction in SSA. Moreover, air temperature played a non-negligible role, and the increase in T_2 max accounted for 12 % of surface O₃ enhancement (Case E). Meanwhile, the increase in PBLHs also contributed about 18 % to the increment in surface O₃ (Case F). As indicated by Case G, the combined effect by multiple factors was larger than the simple summation of individual factors' contributions, or the total percentage of contribution by an individual factor was less than 100 %. This is likely due to the fact that O₃ production is not the linear function of an individual factor's contribution. Complex interplay among different factors may account for the rest of the increase (i.e., 2 %).

It is not surprising that reduction in NO_x emissions brought about an increase in surface O₃ since O₃ formation was dominated by a VOC-limited regime in most parts of the NCP region. Several numerical experiments were conducted to understand the mechanism of reduced PM_2.5 or AOD facilitating the increase in surface O₃. It is known that photolysis rate of NO₂, j(NO₂), plays a critical role in O₃ formation. The parameter j(NO₂) was highly dependent on aerosol optical properties such as AOD and SSA, as well as solar zenith angle (θ) (Dickerson et al., 1997). As shown in Fig. 8a, decreasing AOD was conducive to photolysis of NO₂ due to reduction of attenuated UV radiation entering the PBL. However, weakened scattering or strengthened absorption property of aerosols (i.e., reduced SSA) may attenuate the UV entering the PBL, deaccelerating photolysis of NO₂. Thus, decrease in SSA may counteract the impact associated with decrease in AOD, which may slow down the increase in surface O₃ to some extent. In addition, j(NO₂) showed the highest value at noontime ($\mathit{\theta }=\mathrm{0}{}^{\circ }$ or sec (θ)=1) and tended to decrease when θ became larger (i.e., early morning or late afternoon). Figure 8b further demonstrates that O₃ formation or MDA8 O₃ showed a near-linear increasing trend with j(NO₂). While decrease in PM_2.5 concentrations or AOD strengthened the UV amount entering the PBL, reduction in SSA may counteract the impact of decreased AOD partially. But impact of AOD outpaced that of SSA. Thus, surface O₃ (e.g., MDA8 O₃) still showed a large increase with the combined effect of AOD and SSA over the past several years.

Figure 8 (a) Response of photolysis rate of NO₂, j(NO₂), to different values of aerosol optical depth (AOD) and single-scatter factor (SSA) and (b) change in MDA8 O₃ with j(NO₂) simulated by the MM model for the cases with SSA = 0.95 and AOD varying from 0.5 to 1.25.

Now let us turn our attention to O₃ chemistry in the varying polluted regions. As illustrated in Fig. 9, HO₂ radicals were sensitive to aerosol properties (i.e., AOD and SSA) but the sensitivity was highly reliant on the solar zenith angle (θ). The HO₂ radical was more sensitive to AOD or SSA in the afternoon than in the morning while photolysis rate of HO₂ is more sensitive to AOD or SSA. It is noted that higher net O₃ production is highly associated with the faster decrease in j(O₃) than j(NO₂) in the afternoon (Gerasopoulos et al., 2006). HO₂ radical abundance reduced as aerosol optical property became more absorptive. This indicates that decrease in SSA may cause reduction of HO₂, less NO₂ and then less O₃ production. The HO₂ peak hour was matched well with that of the O₃ peak (around 15:00 LT), further confirming its important role in O₃ formation. Decrease in AOD may accelerate production of HO₂ radicals or slow down their sink, which was conducive to production of NO₂ (Li et al., 2019a) but decrease in SSA may offset its impact if aerosols show strong absorption properties. Meanwhile, strengthened UV associated with weakened aerosol radiative effect was conducive to photolysis of NO₂. As a result, more O₃ is produced. This accounted for a substantial increase in surface O₃ while PM_2.5 decreased over the past several years (2013–2019). The results are consistent with the findings by Li et al. (2019a).

Table 2 Relative percentage contributions of emissions (Case B), AOD (Case C1), SSA (Case D1), air temperature (Case E) and PBLH (Case F) to the change in MDA8 O₃ over the NCP region during 2013–2019.

Figure 9 Response of concentrations of HO₂ to different values of aerosol optical depth (AOD) and single-scatter factor (SSA).

4 Discussions

In this study, a box model NCAR MM with the detailed NO_x–VOC–O₃ chemistry included is utilized to quantify percentage contributions of emissions, aerosol optical properties and meteorological variabilities to increase in surface O₃ over the NCP region during 2013–2019. The findings may provide more scientific evidence to policy makers on developing more effective control strategies on reduction in ambient levels of O₃ as well as exceedance events. However, several points deserve further discussions.

First, the impact of aerosol radiative effect on surface O₃ formation is dependent on not only aerosol abundance (i.e., AOD) but also aerosol scattering and absorption properties (i.e., SSA). Their impacts can be offset to some extent when AOD and SSA show the same change trend (either increase or decrease) or can be strengthened substantially when both AOD and SSA show an opposite change trend. Here the study on the NCP region represents the first case since both AOD and SSA showed a decrease trend over the past several years. Even so, the combined impact of aerosol radiative effect due to reductions in AOD and SSA still contributed 23 % of the total change in surface O₃ in the NCP over the past several years. This reminds us that the impact of aerosol radiative effect could be more substantial if both AOD and SSA show an opposite change trend. Moreover, compared to impact of change in AOD on surface O₃ formation (e.g., Dickerson et al., 1997; Wang et al., 2016a; Xing et al., 2015, 2017), studies on impact of change in SSA on surface O₃ formation are fewer (Dickerson et al., 1997; Mok et al., 2016). Thus, changes of individual aerosol radiative property parameters must be addressed carefully in order to present more accurate quantification of impact of aerosol radiative effect on change in surface O₃.

Second, the MM model does not include aerosol chemistry. As presented above, the MM model as a box model with the detailed O₃–NO_x–VOC relationship allows us to quantify relative contributions of individual factors to increase in surface O₃. Overall, the model results are comparable to those using three-dimensional (3D) chemistry and transport models (CTMs) (e.g., Liu and Wang 2020a, b). For instance, the MM model result indicates that 45 % of increase in surface O₃ was attributed to reduction of anthropogenic emissions of NO_x in the NCP region during 2013–2019, which fell in the range of the results with 3D CTM. Among the 3D modeling studies, Li et al. (2019a) found that anthropogenic emissions contributed about 10 % of change in surface O₃ in summertime from 2013–2017, and Sun et al. (2019) showed the percentage contribution of anthropogenic emissions was 63 % over eastern China. However, there is some substantial difference between the MM model result and that of 3D CTMs in terms of percentage contribution of aerosol radiative effect to changes in surface O₃. The MM model showed that aerosol radiative effect was ranked as the second contributor to the change in surface O₃ in this region. The percentage contribution was larger than that presented by other studies (Li et al., 2019a; Xing et al., 2015). This is partly because this study is focused on the impact on MDA8 O₃ whereas their studies investigated the impact on diurnal variations in surface O₃. In addition, Li et al. (2019a, b) and Liu and Wang (2020a, b) pointed out that aerosol chemistry played the most important role in the enhancement of surface O₃ in this region through modification of HO₂ radicals that produce additional O₃ formation. However, the MM model does not include aqueous-phase chemistry that has been implemented in the 3D meteorology–chemistry models (e.g., Li et al., 2019a; Liu and Wang, 2020a, b), which could be another possible reason for the response to such a difference. Thus, inclusion of detailed aerosol chemistry and observation-based uptake coefficients in a box model like MM is necessary to provide more accurate assessment of the impact of aerosol radiative effect on surface O₃ change.

Third, compared to 3D meteorology–chemistry coupling model(s), the box model does not include complex physical processes such as regional transport, vertical transport, cloud formation, etc. The influence of changing meteorological factors on the change trend in surface O₃ may vary greatly with regions and time. In addition to air temperature and the boundary layer conditions, other meteorological factors such as cloud cover, precipitation, wind fields played an important role in driving the changes in surface O₃ observed in many places in China (Liu and Wang, 2020a). Computational resource and workload that a box model requires are much less than those that a 3D chemical transport model needs. This may allow us to complete a series of designed numerical experiments to quantify the roles of individual factors easily with limited computational resources. It is acceptable to use a box model if terrains are relatively flat in the box, horizontal gradients of emissions and air pollutant concentrations are not strong, and transport in and out reaches a relative equilibrium state. As shown in Figs. S1 and 2, the NCP region defined in this study represents the most polluted part in eastern China; anthropogenic emissions appear to distribute relatively uniformly across the region. To this extent, it is appropriate to examine O₃ formation and its response to changes of different factors such as emissions, meteorological conditions and aerosol radiative properties by using a box model in the NCP region. However, some other physical processes such as long-range transport may exert an important impact on change in surface O₃ (e.g., Han et al., 2018; Gaudel et al., 2018). It is reiterated that the box model results present an ensemble-mean behavior for the given box but need further evaluations by using a complex meteorology–chemistry coupling model such as the Weather Research and Forecasting model with Chemistry (WRF-Chem).

Furthermore, some other important factors may exert an important impact on surface O₃ concentrations, but they are not discussed in this study. Stratospheric intrusion and change in tropospheric O₃ could exert an important impact on O₃ in the atmospheric boundary layer (ABL) and near surface. For instance, Jiang et al. (2015) presented a factor analysis on an O₃ episode observed in the southeast coastal area of China and found that the downward transport of O₃ from the upper troposphere–lower stratosphere (UTLS) region driven by a typhoon was the key factor causing a large increase in surface O₃ by 21–42 ppb. Thus, the impact of tropospheric O₃ should be taken into account when the appropriate observational data are available in the NCP region. Another factor is synoptic patterns. As an example, high concentrations of surface O₃ or O₃ episodes that occurred in the western Mediterranean and central Europe were usually linked with an anticyclone synoptic pattern which led to large-scale subsidence, clear sky and high temperature (e.g., Kalabokas et al., 2013, 2017). In addition, Yin et al. (2019) found that synoptic patterns played a critical role in summer O₃ pollution events in eastern China. Under the control of the zonally enhanced East Asian deep trough, the local hot, dry air and intense solar radiation enhanced the photochemical reactions and produced more O₃. The inter-annual magnitude variations in the domain synoptic patterns may have an important impact on surface O₃, and their impact on the long-term change in surface O₃ needs further investigation.

5 Summary and conclusions

In this study, 7-year-long surface observational air quality data are presented together with satellite retrieval measurements of TCNO₂, AOD and SSA to investigate long-term change trends of surface O₃ over the NCP region in summer from 2013–2019. A comprehensive statistical analysis is completed to explore the relationship of MDA8 O₃ with PM_2.5 concentrations, tropospheric columns of NO₂, AOD and meteorological variables such as T_2 max, surface-reaching shortwave radiation, wind speed and PBLH. A box model representing the O₃–NO_x–VOC relationship is then utilized to quantify the relative contributions of different driving factors to the increase in surface O₃ in the NCP region over the period of 2013–2019.

The observational analysis indicates, while PM_2.5 concentrations continued to decrease with a rate of 9.5 µg m⁻³ a⁻¹, surface O₃ showed an accelerated increase trend at a rate of 4.6 ppb a⁻¹ over the NCP region during summertime from 2013–2019. Both decrease in PM_2.5 and reduction in NO₂ are the two key factors leading to such an increase in surface O₃. The former is closely associated with the attenuation of UV entering the PBL whereas the latter is related to the fact that O₃ photochemical production in the NCP region is dominated by a VOC-limited regime. The trend analysis of satellite retrieval measurements revealed an obvious increase in T_2 max at the rate of 0.34 ^∘C a⁻¹, a rapid decrease in AOD from 1.0 in 2013 to 0.75 in 2019, and a reduction in SSA from 0.95 to 0.93. The changes of both T_2 max and AOD were conducive to photochemical production of O₃ whereas the variability of aerosol scattering–absorption properties (i.e., decrease in SSA) may largely offset the impact of AOD reduction.

The sensitivity studies with the box model MM indicate that reduction of emissions (i.e., NO_x), meteorological conditions and aerosol radiative effect associated with decrease in aerosol concentrations were the three most important factors in driving such a large increase in surface O₃. They accounted for 45 %, 30 % and 23 % of the total increase in surface O₃, respectively, over the NCP region in summertime during 2013–2019. For the meteorological contribution, increases in the PBLH and air temperature (e.g., T_2 max) were responsible for 18 % and 12 % of the total change of surface O₃, respectively. The percentage contribution of aerosol radiative effect (23 %) represented the net changes caused by aerosol concentrations (i.e., AOD) and aerosol radiative properties (scattering–absorption, SSA) (70 % vs. −47 %). The model results further demonstrated that decrease in SSA (i.e., more absorptive) may lead to reduction in HO₂ radicals and NO₂ concentrations, and then less O₃ production, which may largely counteract the impact of the aerosol radiative effect associated with a decrease in AOD.

This study has a strong implication that development of more effective control strategies on surface O₃ reduction needs to consider the impact of aerosol radiative effect as well as the change of aerosol scattering–absorption properties (i.e., AOD and SSA).