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 (PM2.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 (O3) 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 O3. In this study, 7-year surface observations and satellite retrieval data are analyzed to determine the long-term change in surface O3 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 O3. Numerical simulations with the National Center for Atmospheric Research (NCAR) Master Mechanism (MM) model suggest that reduction of O3 precursor emissions and aerosol radiative effect accounted for 45 % and 23 % of the total change in surface O3 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 O3, 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 O3 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 O3 as well as alleviation of national air quality standard exceedance events. Published by Copernicus Publications on behalf of the European Geosciences Union. 2 X. Ma et al.: Rapid increase in summer surface ozone

Abstract. 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 3 ) 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 3 . In this study, 7-year surface observations and satellite retrieval data are analyzed to determine the long-term change in surface O 3 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 3 . Numerical simulations with the National Center for Atmospheric Research (NCAR) Master Mechanism (MM) model suggest that reduction of O 3 precursor emissions and aerosol radiative effect accounted for 45 % and 23 % of the total change in surface O 3 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 3 , 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 3 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 3 as well as alleviation of national air quality standard exceedance events.

Introduction
Elevated ambient levels of ozone (O 3 ) 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 3 is produced by a series of photochemical reactions involving nitrogen oxides (NO x = NO + NO 2 ) and volatile organic compounds (VOCs) in the presence of solar radiation. Ambient levels of O 3 are highly dependent on emissions of O 3 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 3 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 3 -NO x -VOC relationship and complexity of processes governing ambient levels of O 3 , a large uncertainty remains in the determination of the impact of different driving factors on changes in surface O 3 concentrations under the polluted atmospheric conditions. Thus, accurate quantification of relative contributions of individual factors to the large increase in surface O 3 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 3 . With rapid development of industrialization and urbanization, anthropogenic emissions of NO x and VOCs, two major precursors of O 3 formation, have been increasing significantly in China over the past several decades . For instance, tropospheric columns of NO 2 (TCNO 2 ), 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 3 in the lower troposphere. Meanwhile, an increase in surface O 3 at a rate of 2 % a −1 was observed in Beijing from 1995-2005(Ding et al., 2008, and a similar increase with 1-2 ppb a −1 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 3 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 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 . 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 2 but not for VOCs, which showed an increase of 2 % instead over the period of 2013-2017 . As a result, monthly mean PM 2.5 concentrations de-creased 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 3 was observed in the NCP region during this period (Lu et al., 2018;Cooper et al., 2014). The fact that O 3 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 3 is attributed to anthropogenic emission control efforts.
Aerosol radiative effect is another factor imposing a large constraint on change in surface O 3 . Aerosols attenuate surface-reaching solar near-ultraviolet (UV) radiation effectively and reduce photolysis rate of NO 2 , a key parameter in determining O 3 formation. The impact of aerosol radiative effect on photolysis rate of NO 2 or O 3 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 3 but aerosols with a strong absorption property (e.g., mineral dust and soot) may inhibit O 3 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 3 increased at a rate of 3.1 ppb a −1 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 3 due to additional O 3 production associated with reduced sink of hydroperoxy radicals (HO 2 ) (Li et al., 2019a). They pointed out that increase in surface O 3 associated with decrease in PM 2.5 was more prominent than that with reduction of NO x emissions over the NCP region where O 3 formation was dominated by a VOC-limited regime. Liu and Wang (2020a, b) found the reduction of PM emissions increased the O 3 levels by enhancing the photolysis rates and reducing heterogeneous uptake of reactive gases (mainly HO 2 and O 3 ), of which the latter is more important than the former. A similar impact associated with aerosol radiative properties on O 3 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 3 production by up to 18 % compared to BC alone in the Amazon Basin (Mok et al., 2016). Thus, surface O 3 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 3 needs to be addressed. 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 . The NCP has become one of the most polluted regions in China. The highest PM 2.5 concentration reached 900 µg m −3 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 3 exceedance events became more frequent and more serious in the NCP region Lang et al., 2017;Zhai et al., 2019). Hourly surface O 3 concentrations went up to 150.0 ppb and the increase rate reached 3.1 ppb a −1 , 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 3 has become an emerging critical air quality issue in this region Shi et al., 2015). Understanding of the driving factors such as a rise in surface O 3 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 3 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 3 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 3 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 3 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 3 to the key driving factors. The specific objectives include (1) to identify the key factors driving the increase in surface O 3 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 3 in summertime during 2013-2019.

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 3 . 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 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 in- cluded 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 3 and other chemical species (e.g., Liu et al., 2019;Geng et al., 2007).
Photolysis rate j (NO 2 ) 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 2 and other chemical species (e.g., H 2 O 2 , O 3 , NO 3 , N 2 O 5 ) every 15 min dynamically. The TUV model is initialized with the monthly means of AOD, SSA and total columns of O 3 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 2 radicals are important to O 3 formation. HONO photolysis as the primary production of OH radicals and formaldehyde (HCHO) photolysis as the net radical source of HO 2 can lead to major changes in the HO x and NO x budget that may have an important effect on O 3 production and loss (e.g., Aumont et al., 2003;Brasseur and Solomon, 2006;Lin et al., 2012). The role of HO 2 radicals can be determined by the following reactions.
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 3 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 2 , O 2 , H 2 O, NO 2 , O 3 , 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 3 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 3 . 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.  Table 1.  for the BTH, YRD and PRD across eastern China in 2013. While the TCNO 2 was 2 times higher than that observed in North America (Stavrakou et al., 2008), exceedance events of the MDA8 O 3 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 2 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 −3 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 −1 to −11 % a −1 which was similar to that during 2013-2017 (Li et al., 2019a). A similar decrease trend was seen in both TNO 2 (Fig. 2g, h) and in situ NO 2 measurements (Fig. S2). On the other hand, a rapid increase in surface O 3 concentrations was observed in the NCP region over the past several years. The hot spots with MDA8 O 3 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 3 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 3 became more severe and more frequent. The frequency of NAAQS exceedance events for surface MDA8 O 3 (i.e., greater than 160 µg m −3 ) 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). 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 3 formation was dominated by a VOC-limited regime. To better understand the relationship of increase in surface O 3 with the decrease in NO 2 , the change in monthly mean O x (a sum of O 3 and NO 2 ) 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 2 radicals was reduced and more O 3 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 2 and eventually lead to more O 3 production. Importance of aerosol radiative effect in the increase in surface O 3 via acceleration of photolysis of NO 2 can be further evaluated through numerical experiments.
Meteorological conditions are another critical factor affecting O 3 production. Typically, higher air temperature is responsible for higher photochemical reaction rates and more O 3 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.

Yearly changes in surface O 3 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 3 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 3 has become a major air quality concern in summer. June was the month with the highest monthly mean MDA8 O 3 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 3 over the NCP region throughout the period of 2013-2019. Figure 4 shows the yearly changes in monthly means of MDA8 O 3 , PM 2.5 , AOD, SSA, TCNO 2 , T 2 max and PBLH over the NCP region in June from 2013-2019. The change in monthly mean of surface MDA8 O 3 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 3 (4.6 ppb a −1 ) during 2013-2019 was much higher than that observed in the same region during the period of 2005-2015 (1.1 ppb a −1 ) (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 −1 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 2 ( Fig. 4b-d) and in situ NO 2 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 no-ticed feature is that MDA8 O 3 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.
To understand the factors driving the change in surface O 3 , a series of scatter plots are presented to examine the relationships between the surface MAD8 O 3 and individual factors such as aerosol optical properties (i.e., AOD and SSA), TCNO 2 , 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 3 showed two different regimes with an opposite dependence of O 3 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 −3 whereas the second one showed no trend with increasing PM 2.5 when PM 2.5 concentrations were higher than 140.0 µg m −3 . 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 3 photochemical production seemed very minor or even negligible. Instead, O 3 production was suppressed significantly and MAD8 O 3 concentrations were less than 20.0 ppb. In this case, removal of surface O 3 through titration of NO was not effective and surface O 3 showed an increase rather than a decrease trend with increasing NO x concentrations under the strong NO x conditions as indicated by TCNO 2 higher than 40-45 × 10 15 (cm −2 ) in the troposphere. Here the threshold value of 140.0 µg m −3 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 3 . MDA8 O 3 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 3 in the NCP region. On the other hand, MDA8 O 3 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 3 pollution to some extent. The positive correlation between the PBLH and O 3 shown in Fig. 5g represents one case when radiation is stronger and temperature is higher, which is favorable for O 3 formation. Meanwhile, higher PBLH could enhance the transport down of O 3 -enriched air aloft, resulting in an increase in surface O 3 (Reddy et al., 2012). On the other hand, some studies found a negative correlation between the PBLH and O 3 . They claimed that a shallower PBL may suppress the dispersion of pollutants and lead to higher O 3 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 3 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 3 (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 −2 when surface PM 2.5 concentrations reach around 300 µg m −3 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 −3 , TCNO 2 of equal to or less than 5.0 × 10 15 (cm −2 ), T 2 max higher than 28.0 • C and surface-reaching shortwave radiation stronger than 250.0 W m −2 . 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 3 production.

Relative contributions of different driving factors to increase in surface O 3
In this section, the box model MM is utilized to quantify the relative contributions of individual driving factors to the increase in surface O 3 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 3 (Fig. 7), of which the O 3 observations averaged over all the stations in NCP are considered to be the standard observed concentrations. The simulated O 3 peak was about 1 h later than the observation, which was likely due to uncertainty of emission inventory and other meteorological fac-tors. 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 3 . 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 3 over   Table 2. The changes in emissions of O 3 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 3 . But increase in surface O 3 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 3 enhancement (Case E). Meanwhile, the increase in PBLHs also contributed about 18 % to the increment in surface O 3 (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 3 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 3 since O 3 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 3 . It is known that photolysis rate of NO 2 , j (NO 2 ), plays a critical role in O 3 formation. The parameter j (NO 2 ) 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 2 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 2 . Thus, decrease in SSA may counteract the impact associated with decrease in AOD, which may slow down the increase in surface O 3 to some extent. In addition, j (NO 2 ) showed the highest value at noontime (θ = 0 • or sec (θ ) = 1) and tended to decrease when θ became larger (i.e., early morning or late afternoon). Figure 8b further demonstrates that O 3 formation or MDA8 O 3 showed a near-linear increasing trend with j (NO 2 ). 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 3 (e.g., MDA8 O 3 ) still showed a large increase with the combined effect of AOD and SSA over the past several years. Now let us turn our attention to O 3 chemistry in the varying polluted regions. As illustrated in Fig. 9, HO 2 radicals were sensitive to aerosol properties (i.e., AOD and SSA) but the sensitivity was highly reliant on the solar zenith angle (θ). The HO 2 radical was more sensitive to AOD or SSA in the afternoon than in the morning while photolysis rate of HO 2 is more sensitive to AOD or SSA. It is noted that higher net O 3 production is highly associated with the faster decrease in j (O 3 ) than j (NO 2 ) in the afternoon (Gerasopoulos et al., 2006). HO 2 radical abundance reduced as aerosol optical property became more absorptive. This indicates that decrease in SSA may cause reduction of HO 2 , less NO 2 and then less O 3 production. The HO 2 peak hour was matched well with that of the O 3 peak (around 15:00 LT), further confirming its important role in O 3 formation. Decrease in AOD may accelerate production of HO 2 radicals or slow down their sink, which was conducive to production of NO 2 (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 2 . As a result, more O 3 is produced. This accounted for a substantial increase in surface O 3 while PM 2.5 decreased over the past several years (2013)(2014)(2015)(2016)(2017)(2018)(2019). The results are consistent with the findings by Li et al. (2019a).

Discussions
In this study, a box model NCAR MM with the detailed NO x -VOC-O 3 chemistry included is utilized to quantify percentage contributions of emissions, aerosol optical properties and meteorological variabilities to increase in surface O 3 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 3 as well as exceedance events. However, several points deserve further discussions.
First, the impact of aerosol radiative effect on surface O 3 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 3 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 3 formation (e.g., Dickerson et al., 1997;Wang et al., 2016a;Xing et al., 2015Xing et al., , 2017, studies on impact of change in SSA on surface O 3 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 3 .
Second, the MM model does not include aerosol chemistry. As presented above, the MM model as a box model with the detailed O 3 -NO x -VOC relationship allows us to quantify relative contributions of individual factors to increase in surface O 3 . 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 3 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 3 in summertime from 2013, 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 3 . The MM model showed that aerosol radiative effect was ranked as the second contributor to the change in surface O 3 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 3 whereas their studies investigated the impact on diurnal variations in surface O 3 . 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 3 in this region  through modification of HO 2 radicals that produce additional O 3 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 3 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 3 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 3 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 3 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 3 (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 3 concentrations, but they are not discussed in this study. Stratospheric intrusion and change in tropospheric O 3 could exert an important impact on O 3 in the atmospheric boundary layer (ABL) and near surface. For instance, Jiang et al. (2015) presented a factor analysis on an O 3 episode observed in the southeast coastal area of China and found that the downward transport of O 3 from the upper troposphere-lower stratosphere (UTLS) re-gion driven by a typhoon was the key factor causing a large increase in surface O 3 by 21-42 ppb. Thus, the impact of tropospheric O 3 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 3 or O 3 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., 2013Kalabokas et al., , 2017. In addition, Yin et al. (2019) found that synoptic patterns played a critical role in summer O 3 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 3 . The inter-annual magnitude variations in the domain synoptic patterns may have an important impact on surface O 3 , and their impact on the long-term change in surface O 3 needs further investigation.

Summary and conclusions
In this study, 7-year-long surface observational air quality data are presented together with satellite retrieval measurements of TCNO 2 , AOD and SSA to investigate long-term change trends of surface O 3 over the NCP region in summer from 2013-2019. A comprehensive statistical analysis is completed to explore the relationship of MDA8 O 3 with PM 2.5 concentrations, tropospheric columns of NO 2 , AOD and meteorological variables such as T 2 max , surface-reaching shortwave radiation, wind speed and PBLH. A box model representing the O 3 -NO x -VOC relationship is then utilized to quantify the relative contributions of different driving factors to the increase in surface O 3 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 −3 a −1 , surface O 3 showed an accelerated increase trend at a rate of 4.6 ppb a −1 over the NCP region during summertime from 2013-2019. Both decrease in PM 2.5 and reduction in NO 2 are the two key factors leading to such an increase in surface O 3 . The former is closely associated with the attenuation of UV entering the PBL whereas the latter is related to the fact that O 3 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 −1 , 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 3 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 con-ditions 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 3 . They accounted for 45 %, 30 % and 23 % of the total increase in surface O 3 , 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 3 , 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 2 radicals and NO 2 concentrations, and then less O 3 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 3 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).
Data availability. Data used in this paper can be provided by Xiaodan Ma (xaiodanma_nuist@163.com) upon request.
Author contributions. JH came up with the original idea of this study. XM and JH designed the numerical simulations. XM conducted the data analysis and the first draft of the manuscript and JH edited it. TZ, CL, KZ, JX and WX were involved in the scientific interpretation and discussions. All the authors commented on the paper.
Review statement. This paper was edited by Maria Kanakidou and reviewed by Guy Brasseur and one anonymous referee.