Modeling Ozone Source Apportionment and Performing Sensitivity Analysis in Summer on the North China Plain

: In recent years, air quality issues due to ﬁne particulate matter have been su ﬃ ciently treated. However, ozone (O 3 ) has now become the primary pollutant in summer on the North China Plain (NCP). In this study, a three-dimensional chemical transport model (the Nested Air Quality Prediction Model System, NAQPMS) coupled with an online source apportionment module was applied to investigate the sources of O 3 pollution over the NCP. Generally, the NAQPMS adequately captured the observed spatiotemporal features of O 3 during the period of July 1st to August 31st in 2017 on the NCP. The results of the source apportionment indicated that the contributions of local emissions and transport from the NCP accounted for the largest proportion of O 3 , with magnitudes of 25% and 39%, respectively. Compared with those in the average monthly results, the local contribution and regional transport during O 3 episodes on the NCP increased by 7% and 10%, respectively. Based on sensitivity tests, two thresholds of the sensitivity indicator P(H 2 O 2 ) / P(HNO 3 ) were detected, at 0.08 and 0.2. Ozone formation in the urban sites of Beijing, Tianjin, and the southern part of Hebei Province was controlled by VOCs, while the other sites were mainly controlled by NO X . Biogenic emissions contributed approximately 18% to O 3 formation in July in the southwestern part of Hebei Province.


Introduction
Ground-level ozone is a secondary pollutant that predominantly results from chain photochemical reactions involving nitrogen oxides (NO X = NO + NO 2 ), carbon monoxide (CO), and volatile organic compounds (VOCs) with the catalysis of sunlight in the troposphere [1,2]. Photochemical reactions contributed~90% to global surface ozone, 9 times the contribution of the exchange from the stratosphere to the troposphere [3,4]. Involving findings have been comprehensively reviewed that high levels of ozone have adverse effects on human health [5], ecology [6], agricultural productivity [7,8], and potentially affect climate change [9]. equation for terrain-following coordinates. The terrain-following coordinate [21] has advantages in implementing boundary conditions and consider the terrain influence. The carbon-bond mechanism Z (CBM-Z) [22], composed of 71 species and 134 chemical reactions, was implemented in NAQPMS to simulate gaseous chemistry. The RADM mechanism involving 22 species was employed to describe the aqueous chemistry and wet deposition process. The inorganic aerosol chemistry and secondary organic aerosol chemistry in NAQPMS were addressed with the ISSOROPIA1.7 model [23] and a bulk two-product yield parameterization [24], respectively. An accurate radiative-transfer model (TUV version 4.5) with an eight-stream discrete ordinate solver was coupled online with NAQPMS to evaluate the effects of aerosols on photolysis frequencies and tropospheric oxidants [25]. For heterogeneous chemical processes, 12 species and 28 reactions involving dust, sea salt, sulfate, and black carbon were included, and the updated model has recently been successfully applied to simulate mixing between dust and anthropogenic gaseous pollutants in East Asia [26]. The advection scheme used an accurate mass-conservative, peak-preserving algorithm [27]. The dry deposition velocity was parameterized using a scheme proposed by Wesely [28]. The meteorological input for the NAQPMS was provided by the Weather Research and Forecasting Model (WRF) [29], a new generation of prediction models and assimilation systems developed by a joint effort of the National Center for Atmospheric Research (NCAR), the National Centers for Environmental Prediction (NCEP) and other institutes and colleges. WRF has been widely applied in the fields of meteorological research and numerical weather prediction, and it can be used not only to simulate real cases but also to research fundamental physical processes with a single module.
As shown in Figure 1a, the modeling domain centered on Beijing in this study included three nested domains covering all of China, the North Plain, and the NCP region. The horizontal grid resolutions are 27 km for domain 1 (D1) with 148 × 121 grids, 9 km for domain 2 (D2) with 136 × 124 grids, and 3 km for domain 3 (D3) with 178 × 244 grids. Vertically, the model employed terrainfollowing coordinate, and the vertical grid spacing was divided into 30 layers increasing gradually from 20 m near the surface to 20 km at the top with 20 layers in the boundary layer. The simulation period was from July 2017 to August 2017, and the first 3 days were used as the spin-up period to reduce the influence of the initial conditions.

Emission Inventory
The emissions data used in this study included an anthropogenic emissions inventory and biomass inventory data. The anthropogenic emission inventory was obtained from MEIC (Multiresolution emission inventory for China) developed by Tsinghua University (http://www.meicmodel.org/index.html, accessed on 19 November 2019) based on 2016 with a resolution of 0.1° × 0.1°. The MEIC was used in this study to investigate the emission responses to ozone concentration. MEIC is a dynamic technology-based inventory developed for China covering the years from 1990 to recent years by Tsinghua University [30] covering 10 major air pollutants and greenhouse gases (SO2, NOx, CO, NMVOC, NH3, CO2, PM2.5, PM10, BC, and OC) and more than 700 anthropogenic emission sources. Considering that China has been vigorously controlling pollution in recent years, the MEIC inventory includes recent control policies based on the available official reports to provide a reliable description of pollutant emissions. In general, the MEIC is a technologically advanced, reliable, and widely used inventory. Exactly, the latest MEIC has covered the year of 2016, but the data updated in the website is lagging, we, fortunately, get the data collection by personal contact and applied it to this study. Five source categories, i.e., residence, agriculture, industry, power plant and transport, and 25 species of pollutants were included.

Emission Inventory
The emissions data used in this study included an anthropogenic emissions inventory and biomass inventory data. The anthropogenic emission inventory was obtained from MEIC (Multi-resolution emission inventory for China) developed by Tsinghua University (http://www.meicmodel.org/index. html, accessed on 19 November 2019) based on 2016 with a resolution of 0.1 • × 0.1 • . The MEIC was used in this study to investigate the emission responses to ozone concentration. MEIC is a dynamic technology-based inventory developed for China covering the years from 1990 to recent years by Tsinghua University [30] covering 10 major air pollutants and greenhouse gases (SO 2 , NOx, CO, NMVOC, NH 3 , CO 2 , PM 2.5 , PM 10 , BC, and OC) and more than 700 anthropogenic emission sources. Considering that China has been vigorously controlling pollution in recent years, the MEIC inventory includes recent control policies based on the available official reports to provide a reliable description of pollutant emissions. In general, the MEIC is a technologically advanced, reliable, and widely used inventory. Exactly, the latest MEIC has covered the year of 2016, but the data updated in the website is lagging, we, fortunately, get the data collection by personal contact and applied it to this study. Five source categories, i.e., residence, agriculture, industry, power plant and transport, and 25 species of pollutants were included. The largest NOx hotspot is in the main urban region of Beijing, where vehicle exhaust emissions are high. Biogenic emissions were obtained from a biogenic emission model named MEGANv2 (Model of Emission of Gases and Aerosols from Nature) developed by NCAR (http://accent.aero.jussieu.fr/database_table_inventories.php, accessed on 21 December 2019). This model is driven by meteorological fields, plant functional types, leaf area index, and emission factors classified by VOC species. More details can be found in the study of Wang et al. [18]. For the areas of high vegetation coverage in the Yanshan Mountains and Taihang Mountain over the NCP, as shown in Figure 1a, the concentration of BVOC (Biogenic volatile organic compounds) was high in the northern part of the research area; this distribution pattern was the opposite of that for NOx. Anthropogenic volatile organic compounds known as AVOCs (the involving species considered in this study included OLEI (Internal olefin carbons (C = C)), OLET (Terminal olefin carbons (C-C)), PAR (Paraffin carbon), TOL (Toluene) and XYL (Xylene)) are mainly produced by four categories of activity: transport, solvent use, production and storage processes, and combustion processes [31]. Their distribution is similar to that of NOx emissions. In addition, the magnitude of AVOCs is low compared to that of BVOCs, suggesting the non-negligible effect of BVOCs on ozone formation.

Source Apportionment
To qualitatively and quantitatively evaluate the ozone source-receptor relationship on the NCP, an online source apportionment module that attributes pollutants to different geographical locations at each step of the simulation without influencing the standard calculations was implemented in NAQPMS. Compared with the classic sensitivity analysis that turns on and off emissions in targeted ozone production regions, the tagged tracer method provides a different and more efficient measurement of the relative importance of various ozone production regions and lacks the errors introduced by important non-linearities in the transport and fast photochemistry of ozone and its precursors. A detailed description of the methods of the trace tagging module can be found in previous studies [32,33]. In this study, the ten ozone production regions were tagged mainly according to administrative divisions, as were the initial conditions and lateral conditions. As shown in Figure 1b

Model Evaluation
The NAQPMS is driven by the mesoscale numerical weather prediction system (WRF), and the simulation of the meteorological field is vital to the simulation of pollutants, especially to the species sensitive to temperature. Therefore, we evaluate the performance of WRF at 6 stations consisting of 3 rural sites and 3 urban sites with the MICAPS data by calculating several important statistical parameters in Table 1. We can see that WRF rationally recaptured the magnitudes of the air temperature with R up to 0.9, and the MB and the RMSE of most sites are around 4 • C and 5 • C, providing sufficient confidence for the following discussion. To evaluate the performance of NAQPMS in both urban and rural regions of the NCP, both urban and rural sites were chosen in Beijing city, Tianjin city, and Hebei Province. Figures 2 and 3 presented the observed and simulated hourly surface ozone and O x (O 3 + NO 2 ) at the six selected sites over the NCP. The observed data from 1 June 2017 to 31 August 2017 were provided by the China National Environment Monitoring Centre (CNEMC). In general, the NAQPMS rationally recaptured the magnitudes and trends of the ozone pattern and was comparable to previous modeling studies in this region [25]. The correlation coefficients (R) between the observed and simulated ozone and O X were 0.56 and 0.58, respectively. A total of 66.36% of the simulated O 3 values were within a factor of  Figure 4). The statistical parameters of NMB, NME, and RMSE for O 3 ranged from 0.22-0.62, 0.4-0.71, 55.74-84.1, respectively; and the above indicators for O X ranged from 0.11-0.32, 0.27-0.4, and 48.55 to 65.13, respectively. Quantitatively, the simulated ozone concentration overestimated the observations on 22-24 July and 27-28 August. This overestimation is widespread in the current atmospheric transport model in East Asia [34,35]. Studies in MCIS-Asia III revealed that this underestimation is mostly attributed to uncertainties in the treatment of the heterogeneous chemistry and vertical transport in boundary layers in the current CTMs [35].
presented the observed and simulated hourly surface ozone and Ox(O3 + NO2) at the six selected sites over the NCP. The observed data from 1 June 2017 to 31 August 2017 were provided by the China National Environment Monitoring Centre (CNEMC). In general, the NAQPMS rationally recaptured the magnitudes and trends of the ozone pattern and was comparable to previous modeling studies in this region [25]. The correlation coefficients (R) between the observed and simulated ozone and OX were 0.56 and 0.58, respectively. A total of 66.36% of the simulated O3 values were within a factor of 0.5 to 2 of the observations (FAC2), while 88.15% of the simulated Ox values were within the FAC2 ( Figure 4). The statistical parameters of NMB, NME, and RMSE for O3 ranged from 0.22-0.62, 0.4-0.71, 55.74-84.1, respectively; and the above indicators for OX ranged from 0.11-0.32, 0.27-0.4, and 48.55 to 65.13, respectively. Quantitatively, the simulated ozone concentration overestimated the observations on 22-24 July and 27-28 August. This overestimation is widespread in the current atmospheric transport model in East Asia [34,35]. Studies in MCIS-Asia III revealed that this underestimation is mostly attributed to uncertainties in the treatment of the heterogeneous chemistry and vertical transport in boundary layers in the current CTMs [35].   presented the observed and simulated hourly surface ozone and Ox(O3 + NO2) at the six selected sites over the NCP. The observed data from 1 June 2017 to 31 August 2017 were provided by the China National Environment Monitoring Centre (CNEMC). In general, the NAQPMS rationally recaptured the magnitudes and trends of the ozone pattern and was comparable to previous modeling studies in this region [25]. The correlation coefficients (R) between the observed and simulated ozone and OX were 0.56 and 0.58, respectively. A total of 66.36% of the simulated O3 values were within a factor of 0.5 to 2 of the observations (FAC2), while 88.15% of the simulated Ox values were within the FAC2 ( Figure 4). The statistical parameters of NMB, NME, and RMSE for O3 ranged from 0.22-0.62, 0.4-0.71, 55.74-84.1, respectively; and the above indicators for OX ranged from 0.11-0.32, 0.27-0.4, and 48.55 to 65.13, respectively. Quantitatively, the simulated ozone concentration overestimated the observations on 22-24 July and 27-28 August. This overestimation is widespread in the current atmospheric transport model in East Asia [34,35]. Studies in MCIS-Asia III revealed that this underestimation is mostly attributed to uncertainties in the treatment of the heterogeneous chemistry and vertical transport in boundary layers in the current CTMs [35].    As the ozone formation is temperature-dependent, we also simply analyzed the impact of temperature uncertainty inert of the WRF model on the ozone simulation. As shown in Figure 5, the averaged observed ozone generally increased with the increasing air temperature, and the averaged simulated and observed air temperature are 31 ℃ and 26 °C, respectively. The discrepancy caused by uncertainty in the air temperature of ozone is about 30 µg/m 3 , which is within an acceptable range.

Indicator of Ozone-Sensitivity Regime
In this study, the P(H2O2)/P(HNO3) ratio was selected as the sensitivity indicator for studying the ozone formation regimes over the NCP. Ozone formation results from a chain of nonlinear chemistry reactions that are terminated by the cross-reactions of ROX and/or NOX depending on the abundance of NOX. Under high NOX conditions, the termination process is dominated by the reactions of NO2 with OH and RO2, as shown in reaction (R1)-(R2), and this process finally generates so-called NOZ species, the rate constant for reactions R1, R2, R3 can also be found in the reference As the ozone formation is temperature-dependent, we also simply analyzed the impact of temperature uncertainty inert of the WRF model on the ozone simulation. As shown in Figure 5, the averaged observed ozone generally increased with the increasing air temperature, and the averaged simulated and observed air temperature are 31°C and 26 • C, respectively. The discrepancy caused by uncertainty in the air temperature of ozone is about 30 µg/m 3 , which is within an acceptable range.  As the ozone formation is temperature-dependent, we also simply analyzed the impact of temperature uncertainty inert of the WRF model on the ozone simulation. As shown in Figure 5, the averaged observed ozone generally increased with the increasing air temperature, and the averaged simulated and observed air temperature are 31 ℃ and 26 °C, respectively. The discrepancy caused by uncertainty in the air temperature of ozone is about 30 µg/m 3 , which is within an acceptable range.

Indicator of Ozone-Sensitivity Regime
In this study, the P(H2O2)/P(HNO3) ratio was selected as the sensitivity indicator for studying the ozone formation regimes over the NCP. Ozone formation results from a chain of nonlinear chemistry reactions that are terminated by the cross-reactions of ROX and/or NOX depending on the abundance of NOX. Under high NOX conditions, the termination process is dominated by the reactions of NO2 with OH and RO2, as shown in reaction (R1)-(R2), and this process finally generates so-called NOZ species, the rate constant for reactions R1, R2, R3 can also be found in the reference

Indicator of Ozone-Sensitivity Regime
In this study, the P(H 2 O 2 )/P(HNO 3 ) ratio was selected as the sensitivity indicator for studying the ozone formation regimes over the NCP. Ozone formation results from a chain of nonlinear chemistry reactions that are terminated by the cross-reactions of RO X and/or NO X depending on the abundance of NO X . Under high NO X conditions, the termination process is dominated by the reactions of NO 2 with OH and RO 2 , as shown in reaction (R1)-(R2), and this process finally generates so-called NO Z species, the rate constant for reactions R1, R2, R3 can also be found in the reference [22]. NAQPMS is capable of calculating and outputting the values of OH and HO 2 online for each time step and each grid without assumptions. Under low NO X conditions, the dominant processes are self-reactions of HO 2 (R3) and cross-reactions of HO 2 and RO 2 (R4), resulting in the formation of hydrogen peroxides (H 2 O 2 ) and organic peroxides. Hence, the ratio of P(H 2 O 2 )/P(HNO 3 ) can reflect the ambient atmospheric conditions; a higher P(H 2 O 2 )/P(HNO 3 ) ratio indicates a NO X -limited regime, while a lower P(H 2 O 2 )/P(HNO 3 ) ratio corresponds to a VOC-limited regime. This indicator has been applied in other studies [16,36]. Figure 6 illustrates the simulated monthly average concentrations of O 3 over the NCP during the day (LST: 6:00-18:00) and night (LST: 18:00-6:00) in July and August. During the day, the highest concentrations, 190 µg/m 3 , were simulated with NAQPMS in the majority of Beijing in July, and the high ozone levels (>170 µg/m 3 ) extended hundreds of kilometers north in Hebei Province. The simulated O 3 concentration in August was generally lower than that in July, with a maximum value of 170 µg/m 3 . At night, high concentrations of O 3 occurred over the western part of Beijing and in the northern and western parts of Hebei Province in July. In August, the simulated O 3 distribution pattern was similar to that in July; however, the magnitudes of O 3 were lower. In both August and July, the daytime ozone concentration was simulated to be higher than that at night, which can be explained by the chemical reactions involving nitrogen oxide (NO) that occur at night and the physical process of dry deposition that results in the depletion of ozone [37]. There was little discrepancy in the simulated O 3 between day and night in the northern part of the NCP, where the O 3 concentration is high in July, which indicates a weak NO titration reaction in this region. In contrast, more serious pollution with a high concentration of NO occurred in the southern part of the NCP, resulting in a large gap between daytime and nighttime O 3 levels in the simulation.

O3 Source Apportionments for 13 Cities on the NCP
To highlight the sources for these O3 episodes, source apportionments for the monthly mean ozone and ozone-exceedance days (O3 > 160 µg/m 3 ) in 13 cities on the NCP are shown in Figure 7. In general, the contributions of local emissions and transport from the NCP accounted for the largest proportion of O3 during the episodes in most cities except for Tangshan and Qinhuangdao. When O3 episodes occurred, the average local contributions and the contributions from the NCP for the 13 cities are 32% and 49%, respectively. Compared to the results of monthly averaged source apportionments shown in Figure 7a, during O3 episodes, the contributions of local emissions and transport from the NCP increased by up to 7% and 10%, respectively. For further analysis, during O3 episodes, the self-contribution, the contribution from neighboring cities within the NCP region and the contribution from Henan, Shandong, and Shanxi to Beijing were 43%, 22%, and 10%, respectively, while on clean days, the values were 31%, 22%, and 10%, respectively, showing that the effect on Beijing decreases with increasing distance. Beijing's terrain is surrounded by mountains on three

O 3 Source Apportionments for 13 Cities on the NCP
To highlight the sources for these O 3 episodes, source apportionments for the monthly mean ozone and ozone-exceedance days (O 3 > 160 µg/m 3 ) in 13 cities on the NCP are shown in Figure 7. In general, the contributions of local emissions and transport from the NCP accounted for the largest proportion of O 3 during the episodes in most cities except for Tangshan and Qinhuangdao. When O 3 episodes occurred, the average local contributions and the contributions from the NCP for the 13 cities are 32% and 49%, respectively. Compared to the results of monthly averaged source apportionments shown in Figure 7a, during O 3 episodes, the contributions of local emissions and transport from the NCP increased by up to 7% and 10%, respectively. For further analysis, during O 3 episodes, the self-contribution, the contribution from neighboring cities within the NCP region and the contribution from Henan, Shandong, and Shanxi to Beijing were 43%, 22%, and 10%, respectively, while on clean days, the values were 31%, 22%, and 10%, respectively, showing that the effect on Beijing decreases with increasing distance. Beijing's terrain is surrounded by mountains on three sides, which favors the accumulation of local O 3 . Additionally, the emissions of precursors involving NO X and VOCs are high, as shown in Figure 1. The transport from neighboring Henan Province also plays an important role in the O 3 levels in Xingtai and Handan regardless of whether an O 3 episode is occurring, revealing the effect of prevailing southerly winds in these regions.
Atmosphere 2020, 11, x FOR PEER REVIEW 10 of 21 sides, which favors the accumulation of local O3. Additionally, the emissions of precursors involving NOX and VOCs are high, as shown in Figure 1. The transport from neighboring Henan Province also plays an important role in the O3 levels in Xingtai and Handan regardless of whether an O3 episode is occurring, revealing the effect of prevailing southerly winds in these regions. In addition to analyzing the spatial contributions, we also quantified the temporal source contributions. Figure 8 presents the temporal source apportionment contributions for the 13 cities within the NCP region during episodes and on clean days. Day 0 represents contributions from the "current day", while day 1 and day 2 denote contributions from one day and two days prior to day 0, respectively. Among the 13 cities, the contribution of the "current day" of most cities exceeded 50%, and the maximum value found in Beijing exceeded 60%. This indicates that contributions on the day of emission contribute the most to O3 formation, with average values of 44% and 46% on clean and episode days, respectively, and the contribution of the "current day" on episode days is slightly higher than that on clean days. The contributions on the day before the day of emission are the second-largest source, with average values of 27% and 28% on clean and episode days, respectively, revealing that the temporal contribution decreases over time. In summary, the local source emissions and transport from neighboring areas of the "current day" are the dominant factors affecting O3 formation on the NCP. These results indicate that policymakers should respond rapidly to O3 episodes, and both regional collaboration and stricter local measures for emissions control should be considered to meet air quality standards and reduce episodes of high O3 levels. In addition to analyzing the spatial contributions, we also quantified the temporal source contributions. Figure 8 presents the temporal source apportionment contributions for the 13 cities within the NCP region during episodes and on clean days. Day 0 represents contributions from the "current day", while day 1 and day 2 denote contributions from one day and two days prior to day 0, respectively. Among the 13 cities, the contribution of the "current day" of most cities exceeded 50%, and the maximum value found in Beijing exceeded 60%. This indicates that contributions on the day of emission contribute the most to O 3 formation, with average values of 44% and 46% on clean and episode days, respectively, and the contribution of the "current day" on episode days is slightly higher than that on clean days. The contributions on the day before the day of emission are the second-largest source, with average values of 27% and 28% on clean and episode days, respectively, revealing that the temporal contribution decreases over time. In summary, the local source emissions and transport from neighboring areas of the "current day" are the dominant factors affecting O 3 formation on the NCP.
These results indicate that policymakers should respond rapidly to O 3 episodes, and both regional collaboration and stricter local measures for emissions control should be considered to meet air quality standards and reduce episodes of high O 3 levels.  Figure 9 shows the diurnal variation in the simulated concentrations of O3 at 6 selected sites that were evaluated and their source apportionment averaged for July and August in 2017. In general, the diurnal variation in the O3 concentration presented a single peak, with a maximum observed in the afternoon. This peak was caused by the complex photochemical reactions occurring in the presence of sufficient solar radiation. The minimum was observed in the morning when the O3 concentration fell below 30 µg/m 3 . Among the 6 selected sites, the O3 concentration at the two sites in Beijing was the highest, and the O3 concentration simulated in urban sites was higher than that in the rural sites in Beijing and Tianjin.   Figure 9 shows the diurnal variation in the simulated concentrations of O 3 at 6 selected sites that were evaluated and their source apportionment averaged for July and August in 2017. In general, the diurnal variation in the O 3 concentration presented a single peak, with a maximum observed in the afternoon. This peak was caused by the complex photochemical reactions occurring in the presence of sufficient solar radiation. The minimum was observed in the morning when the O 3 concentration fell below 30 µg/m 3 . Among the 6 selected sites, the O 3 concentration at the two sites in Beijing was the highest, and the O 3 concentration simulated in urban sites was higher than that in the rural sites in Beijing and Tianjin.  Figure 9 shows the diurnal variation in the simulated concentrations of O3 at 6 selected sites that were evaluated and their source apportionment averaged for July and August in 2017. In general, the diurnal variation in the O3 concentration presented a single peak, with a maximum observed in the afternoon. This peak was caused by the complex photochemical reactions occurring in the presence of sufficient solar radiation. The minimum was observed in the morning when the O3 concentration fell below 30 µg/m 3 . Among the 6 selected sites, the O3 concentration at the two sites in Beijing was the highest, and the O3 concentration simulated in urban sites was higher than that in the rural sites in Beijing and Tianjin.  Generally, self-contributions accounted for approximately 40% and 60% of the ozone at rural and urban sites in Beijing in the afternoon. Variations in self-contribution were similar to those of the ozone concentration, revealing the significance of local emissions during the day. At night, when the O 3 concentration was low, transport from other regions both within the NCP and outside of the NCP played an important role. In Tianjin, regional transport from emissions outside the NCP occupied a larger proportion than in the sites in Beijing and Hebei Province, accounting for approximately 60% in the early morning, while a self-contribution below 35% was simulated in the rural sites even if the O 3 concentration reached a peak. Influenced by the prevailing southerly winds in the summer, the contribution of transport from regions upwind of Shijiazhuang, Baoding, Tangshan, Cangzhou, and Langfang played an important role in the source apportionment of the two sites in Hebei Province. In summary, at the two sites in Hebei Province, the variation in self-contribution maintained an analogous pace with that of O 3 concentration, with the maximum value exceeding 50%.

Ozone Sensitivity Analysis over the NCP
To understand the ozone precursor relationships over the NCP, we performed the following sensitivity tests to determine the thresholds of the selected ozone formation indicator P(H 2 O 2 )/P(HN O 3 ). P(H 2 O 2 )/P(HNO 3 ) was chosen because this ratio in particular is closely related to the chemical factors that create the division between NOx-sensitive and VOC-sensitive conditions [38].
The sensitivity experiment design scheme is shown in Table 2. In the Base condition, the baseline emission source was used, while in Case 1 and Case 2, a 50% reduction in the NO X emission source and a 50% reduction in the VOC emission source were used, respectively, in the NAQPMS to calculate the ozone formation rate (P(O 3 )), which is defined as follows: where C N and C V denote the simulated ozone concentrations in Case 1 and Case 2, respectively, and C B denotes the simulated ozone concentration when using the baseline emission source in the model. Figure 10 is the scatter plot of the ozone formation indicator (P(H 2 O 2 )/P(HNO 3 )) and the ozone formation rate (P(O 3 )). It should be mentioned that the extracted data were all urban grid points, and the extraction time is the period when the ozone concentration was high (11-13 July 2017 and 8-10 August 2017) to better reflect the relationships between the ozone formation rate and the ozone formation sensitivity indicator. We can see that in Case 1, when NOx emissions were reduced by 50%, the ozone formation rate increased with the increasing P(H 2 O 2 )/P(HNO 3 ) ratio, revealing that the ozone formation regime transformed from a NO X -limited regime to a VOC-limited regime; in contrast, the ozone formation regime transformed from a VOC-limited regime to a NO X -limited regime when the VOC emissions were reduced by 50%. The intersections of the ozone formation rate and the ozone formation indicator shown in Figure 10 are the thresholds for inferring the ozone formation regime, and two values, 0.08 and 0.2, were detected by correlation calculation. We can conclude that ozone was controlled by VOCs when the ratio of P(H 2 O 2 )/P(HNO 3 ) was below 0.08, ozone was controlled by VOCs when the ratio of P(H 2 O 2 )/P(HNO 3 ) was above 0.2, and ozone was controlled by both NO X and VOCs when the ratio was in the range of 0.08 to 0.2, and the results obtained from this study agree well with the results from the research of Tonnesen and Dennis [39], the thresholds of which are 0.06 and 0.2. Except for the ratio of P(H 2 O 2 )/P(HNO 3 ), several other indicators have been developed to diagnose the O 3 -precursor relationship, we have reviewed other photochemical indicators (NOy, HCHO/NO 2 ) and compared them to our results. The transition point between NOx-and VOC-sensitive locations for the NCP occurred at approximately 7.87 ppb NOy, which is close to that for the Barcelona area [40] and San Joaquin Valley [41]. However, the NOy transition value in our study is higher than the result conducted in Lake Michigan, which is considered a scenario with little or no biogenic VOCs [42]. For the ratio of HCHO/NO 2 , Tonnesen and Dennis [40] found that in situ measurements of the ratios of which between 0.8 and 1.8 indicating a "transition" environment where O 3 was both sensitive to radicals and NO X , and the two thresholds we calculated are 0.9 and 2.2. Considering the complexity of the real ambient air, we are supposed to note that differences and uncertainties still exist in applying the indicators to determine the ozone formation regime and thresholds often vary by time and region. considered a scenario with little or no biogenic VOCs [42]. For the ratio of HCHO/NO2, Tonnesen and Dennis [40] found that in situ measurements of the ratios of which between 0.8 and 1.8 indicating a "transition" environment where O3 was both sensitive to radicals and NOX, and the two thresholds we calculated are 0.9 and 2.2. Considering the complexity of the real ambient air, we are supposed to note that differences and uncertainties still exist in applying the indicators to determine the ozone formation regime and thresholds often vary by time and region.  Figure 10. Scatter plot for the ratio of P(H2O2)/P(HNO3) and the difference in P(O3); blue indicates the sensitivity test with a 50% reduction in NOX emissions, and red indicates the sensitivity test with a 50% reduction in VOC emissions. Figure 11 illustrates the monthly distribution of ozone-precursor relationships over the NCP with the indicator discussed above at different ozone levels. On clean days, when the ozone concentration was below 100 µg/m 3 , the majority of Beijing city and Tianjin city and several cities lying in the southwestern part of Hebei province experienced a VOC-limited regime, while the northern rural part of Beijing, Chengde city and Zhangjiakou city were controlled by NOx. The transition regime was small. In episodes when the ozone concentration was above 200 µg/m 3 , a transition in O3 production from a VOC-limited to a NOX-limited regime occurred near the urban area, which can be attributed to the poor meteorological diffusion conditions and high levels of pollutant emissions from urban areas. The precursors transported to the nearby regions corresponded to the change in the dominant factors that affect ozone formation. Furthermore, we also found that the opposite conditions occurred in other locations, which may be explained by the ; blue indicates the sensitivity test with a 50% reduction in NOX emissions, and red indicates the sensitivity test with a 50% reduction in VOC emissions. Figure 11 illustrates the monthly distribution of ozone-precursor relationships over the NCP with the indicator discussed above at different ozone levels. On clean days, when the ozone concentration was below 100 µg/m 3 , the majority of Beijing city and Tianjin city and several cities lying in the southwestern part of Hebei province experienced a VOC-limited regime, while the northern rural part of Beijing, Chengde city and Zhangjiakou city were controlled by NOx. The transition regime was small. In episodes when the ozone concentration was above 200 µg/m 3 , a transition in O 3 production from a VOC-limited to a NO X -limited regime occurred near the urban area, which can be attributed to the poor meteorological diffusion conditions and high levels of pollutant emissions from urban areas. The precursors transported to the nearby regions corresponded to the change in the dominant factors that affect ozone formation. Furthermore, we also found that the opposite conditions occurred in other locations, which may be explained by the control measures in the local region. Compared to that in July, the ozone formation in August was more controlled by VOCs. On clean days, the southern part of Beijing was almost entirely controlled by VOCs, and Baoding, Shijiazhuang, Handan, and Xingtai were transformed into a VOC-limited regime.
Atmosphere 2020, 11, x FOR PEER REVIEW 14 of 21 control measures in the local region. Compared to that in July, the ozone formation in August was more controlled by VOCs. On clean days, the southern part of Beijing was almost entirely controlled by VOCs, and Baoding, Shijiazhuang, Handan, and Xingtai were transformed into a VOC-limited regime.

Ozone Net Production Analysis over the NCP
We also examined the net chemical production of O3 (chemical production-chemical loss) within the boundary layer over the NCP. Here, the net chemical production, N(O3), is calculated by the equation, [olefin]} in NAQM. The NAQPMS models give the net chemical production as the difference of O3 mixing ratio between the calculation steps of the chemistry module with a process analysis package. The net chemical production was calculated in each grid and then the average was taken for all the selected grids. As shown in Figure 12, in general, the maximum net chemical production of ozone was simulated to be approximately 50 µg/m 3 over the urban sites of Beijing, the eastern part of Tianjin, and the southwestern part of Tangshan; meanwhile, negative ozone net [olefin]} in NAQM. The NAQPMS models give the net chemical production as the difference of O 3 mixing ratio between the calculation steps of the chemistry module with a process analysis package. The net chemical production was calculated in each grid and then the average was taken for all the selected grids. As shown in Figure 12, in general, the maximum net chemical production of ozone was simulated to be approximately 50 µg/m 3 over the urban sites of Beijing, the eastern part of Tianjin, and the southwestern part of Tangshan; meanwhile, negative ozone net chemical production was simulated over the sites adjacent to the urban and rural regions. Compared to the maximum 8-h ozone net chemical formation, the high-value zone of the maximum 1-h zone net chemical production was distributed over a larger scale. This phenomenon is particularly obvious in August, which can be attributed to the rapid photochemical reactions during periods with strong solar radiation. The net chemical production of O 3 in this paper was averaged throughout the day, and the net chemical production during the day may be higher than the net chemical production at night. Some studies on O 3 net chemical production have been conducted in Europe and North America. For example, Connor et al. [43] estimated that the daytime net ozone production in the boundary layer over central Europe during July/August 2000 was within a range of 27.87-83.62 µg/m 3 , and our results are consistent with this research. By the further analysis shown in Figure 1 (NOx, AVOC), we found that the pattern of ozone net production was similar to that of anthropogenic emissions, and pollution sources from anthropogenic emissions played a dominant role in ozone net chemical production, revealing the significant effect of anthropogenic source emissions on the model simulation.
Atmosphere 2020, 11, x FOR PEER REVIEW 15 of 21 chemical production was simulated over the sites adjacent to the urban and rural regions. Compared to the maximum 8-h ozone net chemical formation, the high-value zone of the maximum 1-h zone net chemical production was distributed over a larger scale. This phenomenon is particularly obvious in August, which can be attributed to the rapid photochemical reactions during periods with strong solar radiation. The net chemical production of O3 in this paper was averaged throughout the day, and the net chemical production during the day may be higher than the net chemical production at night. Some studies on O3 net chemical production have been conducted in Europe and North America. For example, Connor et al. [43] estimated that the daytime net ozone production in the boundary layer over central Europe during July/August 2000 was within a range of 27.87-83.62 µg/m 3 , and our results are consistent with this research. By the further analysis shown in Figure 1 (NOx, AVOC), we found that the pattern of ozone net production was similar to that of anthropogenic emissions, and pollution sources from anthropogenic emissions played a dominant role in ozone net chemical production, revealing the significant effect of anthropogenic source emissions on the model simulation.

Contributions of Biological Sources
To investigate the impact of BVOCs on ozone formation over the NCP, two simulations were conducted using identical model configurations except that the first scenario considered both anthropogenic and BVOC emissions, while the second included only anthropogenic emissions. Figure 13 illustrates the spatial distribution of the differences in the maximum 1-h average, maximum 8-h average, and monthly average ozone concentrations simulated from the NAQPMS. The largest contributions of BVOCs to the maximum 1-h average, maximum 8-h average, and monthly average were over 30 µg/m 3 in July and August, suggesting a significant impact of BVOC emissions on ozone formation over the NCP; this result is consistent to a certain degree with the research conducted in the YRD region [44]. In July, the hotspot of BVOC contributions to the maximum 1-h average and the maximum 8-h average was mainly concentrated in the southwestern part of Hebei Province, in the range of 36 • N-39 • N, 113 • E-115 • E, and the hotspot of BVOC contributions to the monthly average was concentrated in a larger region. The spatial distribution of BVOC contributions to the maximum 1-h average and maximum 8-h average and monthly average in August is similar to that in July but with lower values due to the weaker solar radiation and lower air temperature in August. The hotspot for the monthly average BVOC contributions is only centered on the southwestern corner of Hebei Province. As shown in Figure 14, there is no major spatial distribution discrepancy between the absolute BVOC contributions and the relative BVOC contributions. The largest contributions of BVOCs to the maximum 1-h average, maximum 8-h average, and monthly average were approximately 18% in July and approximately 12% in August.
Atmosphere 2020, 11, x FOR PEER REVIEW 16 of 21 To investigate the impact of BVOCs on ozone formation over the NCP, two simulations were conducted using identical model configurations except that the first scenario considered both anthropogenic and BVOC emissions, while the second included only anthropogenic emissions. Figure 13 illustrates the spatial distribution of the differences in the maximum 1-h average, maximum 8-h average, and monthly average ozone concentrations simulated from the NAQPMS. The largest contributions of BVOCs to the maximum 1-h average, maximum 8-h average, and monthly average were over 30 µg/m 3 in July and August, suggesting a significant impact of BVOC emissions on ozone formation over the NCP; this result is consistent to a certain degree with the research conducted in the YRD region [44]. In July, the hotspot of BVOC contributions to the maximum 1-h average and the maximum 8-h average was mainly concentrated in the southwestern part of Hebei Province, in the range of 36° N-39° N, 113° E-115° E, and the hotspot of BVOC contributions to the monthly average was concentrated in a larger region. The spatial distribution of BVOC contributions to the maximum 1-h average and maximum 8-h average and monthly average in August is similar to that in July but with lower values due to the weaker solar radiation and lower air temperature in August. The hotspot for the monthly average BVOC contributions is only centered on the southwestern corner of Hebei Province. As shown in Figure 14, there is no major spatial distribution discrepancy between the absolute BVOC contributions and the relative BVOC contributions. The largest contributions of BVOCs to the maximum 1-h average, maximum 8-h average, and monthly average were approximately 18% in July and approximately 12% in August.   Interestingly, the spatial pattern of BVOC contributions did not align with that of the absolute ozone distribution or with biogenic emissions because ground-level ozone is formed by complex nonlinear reactions between NOx and VOCs. In summer, the amount of BVOC emissions in suburban districts with high vegetation coverage is larger than that in urban districts, and O3 formation in suburban districts and urban districts is commonly controlled by NOX and VOC. When BVOC emissions in urban districts increase, the increase in O3 levels in urban districts is more pronounced than that in suburban districts, and in some southern regions with large amounts of BVOC emissions, the O3 formation regime will be transformed when BVOC emissions exceed a certain threshold. The discussion above confirms the significance of the O3 formation regime, and we will discuss this issue in detail in Section 4. In addition, it is worth mentioning that the addition of the BVOC emission inventory effectively improves the performance of the simulation shown in Figures 2-4.

Conclusions
The formation of O3 is sophisticated by nonlinear chemistry that involves multiple ambient pollutants. The measures for O3 pollution prevention should be refined, and the core issue is to determine the main sources and control factors. In this study, we focus on the NCP region, and the Nested Air Quality Prediction Modeling System (NAQPMS) coupled with an online source-tagged model was applied to quantify the O3 source-receptor relationship, to identify the ozone sensitivity threshold with the sensitivity indicator P(H2O2)/P(HNO3) and to investigate the contribution of biogenic emissions to ozone formation. The following conclusions can be drawn: In general, the simulated O3 was higher during the daytime and in July than at night and in August; the highest concentration, 190 µg/m 3 , was distributed across the majority of Beijing in July; and the maximum values extended to the region over Chengde and Baoding in Hebei Province. Monthly average relative contribution (%) of BVOCs to ozone formation in July and August 2017 over the NCP (a1-a6); the left column represents the maximum 1-h average, the middle column represents the maximum 8-h average, the right column represents the monthly average, and the rows from top to bottom represent July and August.
Interestingly, the spatial pattern of BVOC contributions did not align with that of the absolute ozone distribution or with biogenic emissions because ground-level ozone is formed by complex nonlinear reactions between NOx and VOCs. In summer, the amount of BVOC emissions in suburban districts with high vegetation coverage is larger than that in urban districts, and O 3 formation in suburban districts and urban districts is commonly controlled by NO X and VOC. When BVOC emissions in urban districts increase, the increase in O 3 levels in urban districts is more pronounced than that in suburban districts, and in some southern regions with large amounts of BVOC emissions, the O 3 formation regime will be transformed when BVOC emissions exceed a certain threshold. The discussion above confirms the significance of the O 3 formation regime, and we will discuss this issue in detail in Section 4. In addition, it is worth mentioning that the addition of the BVOC emission inventory effectively improves the performance of the simulation shown in Figures 2-4.

Conclusions
The formation of O 3 is sophisticated by nonlinear chemistry that involves multiple ambient pollutants. The measures for O 3 pollution prevention should be refined, and the core issue is to determine the main sources and control factors. In this study, we focus on the NCP region, and the Nested Air Quality Prediction Modeling System (NAQPMS) coupled with an online source-tagged model was applied to quantify the O 3 source-receptor relationship, to identify the ozone sensitivity threshold with the sensitivity indicator P(H 2 O 2 )/P(HNO 3 ) and to investigate the contribution of biogenic emissions to ozone formation. The following conclusions can be drawn: In general, the simulated O 3 was higher during the daytime and in July than at night and in August; the highest concentration, 190 µg/m 3 , was distributed across the majority of Beijing in July; and the maximum values extended to the region over Chengde and Baoding in Hebei Province.
The results of the source apportionment indicated that the contributions of local emissions and transport from the NCP accounted for the largest proportion of O 3 , with magnitudes of 25% and 39%, respectively. Compared with the monthly averaged results, the local contribution and regional transport within the NCP during O 3 episodes increased by 7% and 10%, respectively. Considering the diurnal variation in the source apportionment, we found that the local contribution was the highest in the daytime, while at nighttime, the contribution of regional transport increased.
Based on the sensitivity tests, two thresholds, 0.08 and 0.2, were detected. The urban sites of Beijing, Tianjin, and the southern part of Hebei Province were controlled by VOCs, with a P(H 2 O 2 )/P(HNO 3 ) ratio below 0.08, while the other sites were mainly controlled by NO X , with a P(H 2 O 2 )/P(HNO3) ratio above 0.2. Similar to the distribution of anthropogenic emissions, the maximum ozone net chemical production was simulated to be approximately 50 µg/m 3 over the urban sites in Beijing, the eastern part of Tianjin, and the southwestern part of Tangshan, revealing the significant effect of anthropogenic emissions on the net chemical production of ozone.
Biological sources exert an important impact on O 3 formation. The maximum absolute contributions for the maximum 1-h average, maximum 8-h average, and monthly average were over 30 µg/m 3 , and a hotspot was centered on southwestern Hebei Province.
In this study, we applied the NAQPMS to analyze the main sources and control factors for O 3 , but further improvement is needed to compensate for the deficiencies in this work. The uncertainty of the emission inventory affects the simulation performance of the NAQPMS and further affects the sensitivity test results. To obtain a more scientific sensitivity threshold, it is necessary to compare the model results with the observed data combined with a box model and to conduct more sensitivity tests. Third, the ambient aerosol concentration has changed considerably recently, and the effects of radiation and heterogeneous chemistry related to aerosols on O 3 formation should be considered in future work.