The contribution of industrial emissions to ozone pollution: identified using ozone formation path tracing approach

Wintertime meteorological conditions are usually unfavorable for ozone (O3) formation due to weak solar irradiation and low temperature. Here, we observed a prominent wintertime O3 pollution event in Shijiazhuang (SJZ) during the Chinese New Year (CNY) in 2021. Meteorological results found that the sudden change in the air pressure field, leading to the wind changing from northwest before CNY to southwest during CNY, promotes the accumulation of air pollutants from southwest neighbor areas of SJZ and greatly inhibits the diffusion and dilution of local pollutants. The photochemical regime of O3 formation is limited by volatile organic compounds (VOCs), suggesting that VOCs play an important role in O3 formation. With the developed O3 formation path tracing (OFPT) approach for O3 source apportionment, it has been found that highly reactive species, such as ethene, propene, toluene, and xylene, are key contributors to O3 production, resulting in the mean O3 production rate (PO3) during CNY being 3.7 times higher than that before and after CNY. Industrial combustion has been identified as the largest source of the PO3 (2.6 ± 2.2 ppbv h−1), with the biggest increment (4.8 times) during CNY compared to the periods before and after CNY. Strict control measures in the industry should be implemented for O3 pollution control in SJZ. Our results also demonstrate that the OFPT approach, which accounts for the dynamic variations of atmospheric composition and meteorological conditions, is effective for O3 source apportionment and can also well capture the O3 production capacity of different sources compared with the maximum incremental reactivity (MIR) method.


INTRODUCTION
Air pollution has caused widespread concern due to human health risks and economic losses 1,2 . Over the past few years, the Chinese government has made great efforts to improve air quality, especially by reducing fine particulate matter (PM 2.5 ) concentrations. From 2013 to 2017, PM 2.5 concentrations decreased by more than 30% across the country 3 . However, ozone (O 3 ) concentrations show an increasing trend (around 2.9 ppbv yr −1 ) in many regions over the period 2013-2019 4

.
Tropospheric O 3 is produced via the photochemical cycle of nitrogen monoxide (NO)-O 3 -nitrogen dioxide (NO 2 ), which is sped up in the presence of volatile organic compounds (VOCs), resulting in an accelerated O 3 production rate in the environment 5 . Xiong et al. 6 investigated the VOCs pollution characteristics in Chengdu and found that combustion sources contributed more to VOCs in winter (23.9%) than that in summer (12.1%), and the main sources of VOCs were natural gas (NG)/liquefied petroleum gas (LPG) usage and industry. Seco et al. 7 conducted seasonal (winter and summer) measurements at a forest site and found that biogenic emissions of isoprene and monoterpenes occurred only during the warmest summer months. Recent studies have found that wintertime VOC concentrations in Heibei province are 1.5 to 2 times higher than those in summer 8,9 . Therefore, the potential source of O 3 in winter should be different from that in summer because O 3 pollution in most Chinese cities is located in a VOC-limited regime and the emission sources and components of VOCs vary greatly 10,11 . In addition, photochemical reactions greatly depend on solar radiation and temperature 12 . Highly reactive VOCs (k OH > 10 −11 cm 3 molecule −1 s −1 , photochemical age assumed to be 12 h) may lose~80% before being sampled due to photochemical loss in summer 13 . Previous studies investigated photochemical losses of VOCs in summer and found that total VOCs were underestimated by about 18.4-23.3% due to photochemical loss, while highly reactive species (isoprene, ethene, propene) were underestimated by about 30.0-61.9% 14,15 . Thus, the importance of these highly reactive VOCs to O 3 formation might be overlooked in summer due to the photochemical losses based on observation-based modeling (OBM) 16,17 . The specific atmospheric conditions such as weak solar radiation and low temperature in winter should provide us a good opportunity to characterize the highly reactive VOCs, in particular, their roles in O 3 production. However, most of the previous studies have mainly focused on summertime O 3 pollution, while wintertime O 3 pollution has attracted insufficient attention 12,18 .
Shijiazhuang (SJZ), the provincial capital of Hebei, is an important industrial base with 10 million residents. It is also an important site in the southwest transport channel of air pollutants in the North China Plain (NCP) 19 . Guan et al. 20 investigated the characteristics and sources of summertime VOCs in Shijiazhuang using an offline sampling method. They found that oxygenated VOCs (OVOCs) accounted for 37.9% of the total VOCs, followed by alkanes (33.9%), unlike the results in Beijing 21,22 , Chengdu 6 , etc., where alkanes had the highest contribution (42.6-44.5%) to total VOCs. Meanwhile, petrochemical (24.2%) and other industrial sources (15.2%) were the main emission sources of VOCs in SJZ 20 . However, less attention has been paid to O 3 pollution in SJZ at present although it is an urgent issue from the perspective of both regional transport of air pollutants and local photochemistry.
When we implement O 3 pollution control, it is crucial to determine the main regional and sectoral sources of O 3 23 . The trajectories analysis of air masses can identify the regional sources of surface O 3 , but not account for the chemical processes 24,25 , thus limiting its application in quantifying O 3 sources. Emission-based model (EBM), which predicts the transport and photochemical reactions of O 3 as well as its precursors utilizing air quality models with a known emissions inventory, is a powerful tool and has been widely used for O 3 source apportionment 26,27 . The OBM, which simulates O 3 formation processes constrained by measurements of the concentrations of O 3 , oxides of nitrogen (NO x ), sulfur dioxide (SO 2 ), carbon monoxide (CO), nitrous acid (HONO), VOCs, the photolysis rate constant (J), temperature (T), pressure (P), and relative humidity (RH), can avoid the uncertainties caused by emission inventories and the simulated dynamics of the boundary layer when compared with the EBM 28 . Thus, it has attracted much attention in China, in particular, due to the relatively insufficient research on emission inventories 25 . It has been well recognized that VOC species have different reactivities, resulting in different O 3 production rates (P O3 ). However, previous studies have rarely quantified the contribution of different VOC sources to O 3 production rates 29,30 . At present, a propylene-equivalent concentration method and a maximum incremental reactivity (MIR) method are usually used to calculate the reactivity of VOCs and the contributions of chemical species and sources to the ozone formation potential (OFP) 31,32 , because the production of O 3 is VOCs-limited in the most of urban areas of China 33,34 . These two methods estimate the amount of O 3 formed under optimum or ideal conditions, which may differ from the amount formed in the actual atmosphere. Ling et al. 35 proposed an O 3 source apportionment method combining positive matrix factorization (PMF) and OBM simulations, in which the OBM simulation was driven by the PMF extracted concentrations to calculate the relative O 3 reduction efficiency (RORE) or relative incremental reactivity (RIR) of a single given VOC source. However, using a given VOC source as OBM input may result in differences in both the concentrations of free radicals and P O3 in the model from that in the actual situation. Therefore, there is a need to develop a more reliable method to apportion O 3 sources and understand wintertime ozone pollution.
In this study, we propose an ozone formation path tracing (OFPT) approach for P O3 in the OBM under real observational conditions and combine it with the VOCs-PMF to quantify the sources of the O 3 production rate. Based on an observation campaign in SJZ during the Chinese New Year (CNY) in 2021, we investigated the O 3 pollution in SJZ from the perspective of both regional transport of air pollutants and local photochemistry. The contributions of anthropogenic VOCs, especially aromatics and alkenes with high reactivity, to O 3 production and formation sensitivity regime were discussed. The OFPT method was used to identify the O 3 sources from the perspective of O 3 production. Meanwhile, the OFPT method was compared with the traditional one based on the MIR method. To the best knowledge, this is the first time report about the wintertime O 3 pollution in SJZ and source apportionment in light of the P O3 based on the OFPT method. Our results highlight the role of highly reactive VOCs from industrial emissions in O 3 pollution in SJZ.

RESULTS
Air quality and meteorological conditions during the observation Figure 1 shows the time series of meteorological parameters, including RH, T, solar radiation (SR), wind speed (WS), and wind direction (WD), and air pollutants (VOCs, NO x , O 3 , SO 2 , PM 2.5 ) during our observations. The RH was 45.0 ± 10.8% during the Spring Festival period, which was higher than that before and after CNY (22.8 ± 14.7% and 28.0 ± 18.2%, respectively). The concentration of PM 2.5 showed a similar trend to the RH, i.e., higher (187.3 ± 54.7 µg m −3 ) during CNY than that before and after CNY (46.7 ± 22.2 and 41.6 ± 25.3 µg m −3 , respectively). The temperature ranged from 0 to 26.9°C, with a mean value of 8.3 ± 6.3°C during the observation period. The maximum daily 8-hour average (MDA8) O 3 concentration increased significantly (P < 0.05) during CNY (48.8 ± 3.7 ppbv) when compared with that before (40.3 ± 3.2 ppbv) and after CNY (38.4 ± 2.7 ppbv), respectively. The maximal hourly O 3 concentration even reached 75.5 ppbv during CNY, which is close to the Chinese second-order air quality standard (160 µg m −3 , corresponding to 77.3 ppbv at 281 K and 1008 hPa). The WS was 1.5 ± 0.7 m s −1 during CNY, which was slightly lower than 1.7 ± 0.7 m s −1 before and 2.2 ± 0.9 m s −1 after CNY. Small wind speed is favorable for the accumulation of pollutants. Previous studies found that local emissions and chemistry significantly contribute to the accumulation of air pollutants when wind speeds were below 2 m s −1 36,37 . This implies that local emissions and chemistry may be the important contributors to O 3 pollution during CNY. In addition, O 3 pollution usually occurs in warm seasons under conditions of intense solar radiation and high ambient temperature 38 . The increased O 3 concentration during CNY implies vigorous photochemistry although the meteorological conditions (e.g., lower T and SR) are not in favor of the occurrence of O 3 pollution in winter.
As can be seen in Fig. 1f and Supplementary Table 1, the concentrations of alkanes, alkenes, alkyne, and aromatics increased significantly (P < 0.05) during CNY, with mean values of 23.3 ± 7.3, 7.3 ± 3.6, 5.2 ± 1.5, and 3.5 ± 1.7 ppbv, respectively, when compared with those before and after CNY period. The NO x concentrations showed a slightly decreasing trend with the mean values of 21.7 ± 11.3, 20.9 ± 9.9, and 17.5 ± 12.3 ppbv before, during, and after CNY, respectively. The NO concentration was lower during CNY than that before and after CNY. This may be attributed to the enhanced generation rates of alkyl peroxide radicals (RO 2 ) and hydrogen peroxide radicals (HO 2 ) from VOCs during CNY, which subsequently will promote the conversion of NO to NO 2 . This is generally consistent with the slightly higher NO 2 concentrations during CNY than that before and after CNY. It should be noted that the NO x concentration is slightly higher during CNY compared to that (17.2 ± 7.3 ppbv) in the previous 4 days (Fig. 1e). This highlights the important role of VOCs in O 3 pollution during CNY. As shown in Supplementary Fig. 1, CO concentrations increased significantly during CNY when compared with those before and after CNY, which implies the intensive combustion emissions from industry and/or residents in Shijiazhuang 39 .

Regional transport vs local photochemistry
Although the mean wind speed on the surface was lower than 2 m s −1 during our observations, regional transportation of air pollutants (O 3 and VOCs) might still be important. Figure 2 shows the potential source contribution of VOCs and O 3 in different periods. VOCs and O 3 have similar patterns of geographical sources regardless of the observation period although a slight difference is observable. If O 3 mainly results from transport, it should show a different distribution pattern from that of VOCs, which are mainly from local emissions. Thus, the observed O 3 pollution event at our observation site should be highly related to the photochemistry of VOCs besides the transport of O 3 40 . Figure  2c, f and i show the fields of geopotential height and wind at a pressure of 925 hPa before, during, and after CNY, respectively. It can be seen that the pressure system in the NCP was influenced by the abnormal weakening of the West Pacific subtropical highpressure system in the southern regions ( Supplementary Fig. 2) and the temporary appearance of a low-pressure system in the northern regions of the observation site (dash box in Fig. 2f). Subsequently, the clean air masses from the northwest monsoon in the observation area are weakened, resulting in unfavorable dilution of local pollutants. Meanwhile, the area between the trough of low pressure and the ridge of high pressure is not conducive to local pollutants diffusion. This unfavorable meteorological condition can be further analyzed using the results of air mass trajectory. Supplementary Fig. 3a, d and g show the 24-h backward trajectories of air mass at a height of 0.5 km before, during, and after CNY, respectively. The air mass was clustered to 3 trajectories in these three periods, according to the variation of the total variation of spatial (TVS) to a possible number of clusters ( Supplementary Fig. 3b, e, h). Before and after CNY, long-distance trajectories from north-related directions (routes 1 and 1, 2, respectively) contributed 63.5-74.5% to the air masses (Supplementary Fig. 3a, g), while south and east directions (routes 2 and 3) with short-range trajectories dominated the air masses during CNY (69.1%, Supplementary Fig. 3d). As shown in Supplementary  Fig. 3f, clusters 2 and 3 were near the ground surface with shortrange trajectories during CNY. Meanwhile, cluster 1 and clusters 1&2 were transmitted at a high height during the periods before and after CNY, respectively. Long-range transport generally has higher wind speeds and a stronger dilution ability to air pollutants than short-range transport 41 . It should be noted that the southward ground surface winds only accounted for 24.0% during CNY (Fig. 1b), while the southward trajectories account for 45.8% during the same period ( Supplementary Fig. 3d). This should be ascribed to their different heights, i.e., the start point of the air mass trajectory is 500 m from the ground (much less susceptible to topographic influence), while the wind direction is measured at 25 m at the observation station.
Interestingly, both the concentrations of O 3 and VOCs well correlated with the relative contribution of short-range air masses ( Supplementary Fig. 4). The highest contributions of the shortrange air masses to VOCs and O 3 were observed during CNY. In addition, high O 3 concentrations generally come from the south direction, which coincided with the pollution distribution in SJZ (Supplementary Figs. 5 and 6). These results illustrate that the O 3 pollution event during CNY should be mainly affected by the lowspeed air mass intrusion from the southern regions of SJZ, leading to the accumulation of air pollutants. This can be ascribed to the fact that (1) the sudden changes in the meteorological system  43 . Large-scale changes in meteorological conditions can affect local air quality. Therefore, unfavorable meteorological and topographic conditions lead to the occurrence of wintertime O 3 pollution, which provides favorable conditions for the photochemistry of VOCs.
Ozone formation sensitivity and contribution of reactive VOCs to O 3 formation The empirical kinetic modeling approach (EKMA) curves are widely used to reveal the dependence of O 3 formation on its precursors after the local pollution levels and meteorological parameters have been accounted for. As shown in Supplementary Fig. 7, the model well-predicted O 3 concentrations with an R 2 of 0.97 during the observations. Meanwhile, the simulated OH and HO 2 concentrations are also at reasonable levels 44 . Figure 3a shows the O 3 formation sensitivity during our observations. O 3 formation was located in a VOC-limited regime as expected during our observations (Fig. 3a), while the ratio of VOCs/NO x during CNY was significantly higher than those before and after CNY. During CNY overlapping the COVID-19 lockdown in 2020, it was also found that an emission reduction of NO x led to an increase in O 3 concentrations with a constant emission of VOCs in the North China Plain 18 . This means that O 3 pollution during CNY in this study should be mainly ascribed to the enhanced VOC emissions. Furthermore, the emitted highly reactive species may also react with nitrate radicals (NO 3 ) in the afternoon when NO concentrations are low, which reduces the atmospheric NO x concentration and leads to a weak titration capability ( Supplementary Fig. 8) 45 .
Figure 3b-d shows the mean diurnal variations of the P O3 and the D O3 during the three periods. The mean O 3 production rates (P O3 ) were 2.2 ± 1.6, 6.4 ± 5.3, and 1.4 ± 1.2 ppbv h −1 before, during, and after CNY, respectively. The O 3 production rates before and after CNY (Fig. 3b, d) were comparable with that in Beijing in winter (~3 ppbv h −1 ) 44 , while they were lower than that in Beijing in summer (10.7 ppbv h −1 ) 40 . During the observation, alkenes were the dominant contributors (53.5 ± 14.5%) to the P O3 , followed by oxygenated VOCs (OVOCs, 20.1 ± 11.8%), aromatics (14.1 ± 4.7%) and alkanes (11.3 ± 3.7%). Similarly, alkenes accounted for 64.7 ± 8.7% of the total P O3 during CNY, followed by alkanes (12.5 ± 3.5%), aromatics (11.4 ± 2.6%), and OVOCs (10.1 ± 2.9%). The relative contribution of alkenes and aromatics to O 3 formation increased around 4.9 and 3.1 times, respectively, during CNY when compared with before and after CNY periods. As shown in Supplementary Fig. 9, alkenes had the highest RIR value (0.43) during CNY, followed by aromatics (0.10), and OVOCs (0.08). A previous study found that OVOCs play an important role in RO x (RO x = OH + HO 2 + RO 2 ) formation (22-44%), subsequently, affecting the O 3 production rate 46 . Thus, the RIR value of OVOCs was lower than that of alkenes and aromatics, but higher than that of alkanes (0.06) during CNY, as shown in Supplementary Fig. 9a. In terms of a single VOC compound, ethene, propene, toluene, and xylene showed high RIR values ( Supplementary Fig. 9b-d). These four compounds accounted for 46.6 ± 7.1% of P O3 during CNY and 38.7 ± 10.1% of P O3 during the whole observation period. Thus, the O 3 pollution event during CNY should be related to the increase in highly reactive VOCs. It should be noted that the peak values of P O3 usually appeared at around 11:00 a.m. This can be explained by the promotion effect of HONO on O 3 formation because the photolysis of HONO is an important source of OH radicals in winter 47 . This is also supported by the maximal photolysis rate of HONO that appeared at 10:00 a.m. as shown in Supplementary Fig. 10.
The reaction between NO 2 and OH was the major contributor to O 3 destruction, which contributed 66.2% of the total D O3 during the whole observation period, followed by the photolysis of O 3 (27.8%). These values are different from those summertime values

Source apportionment of VOCs and O 3 production
To identify the sources of anthropogenic VOCs connecting with O 3 pollution during CNY, source apportionment of the observed VOCs was performed using the PMF model. The source profiles are shown in Supplementary Fig. 11a, and the details about source identification are described in Supplementary Discussion 1. Briefly, six types of VOC sources were identified, including biomass burning mixed with fireworks, industrial combustion, cooking, vehicle exhaust, solvent use, and pharmaceutical process exhaust. Figure 4a shows the time series of each source's contribution to VOCs. As shown in Supplementary  Fig. 12). Besides these six sources of VOCs in winter, a biogenic source characterized by high loading of isoprene was also identified in summer. The relative contributions of these sources in summer are in the following order, vehicle exhaust (26.2%) > industrial combustion (24.0%) > cooking (21.7%) > solvent use (9.2%) > biomass burning (9.0%) > biogenic emissions (8.1%) > pharmaceutical exhaust (1.2%). Thus, the contribution of industrial combustion to the observed VOCs (37.8%, 8.9 ± 5.6 ppbv) was higher in winter than that in the summer, i.e., 24.0% ( Supplementary Fig. 12c, 2.3 ± 0.9 ppbv). In addition, it can be seen from Supplementary Table 1 that the contributions of OVOCs and aromatics to the total VOCs in summer (25.6% and 11.1%, respectively) are higher than those in winter (17.9% and 7.3%, respectively), while the contributions of alkenes and alkyne to the total VOCs are lower in summer (6.3% and 4.4%, respectively) than those in winter (13.9% and 9.1%, respectively). These results indicate that the role of highly reactive species such as alkenes and alkyne in O 3 formation is visible due to the weak irradiation in winter, while it might be underestimated in summer.
To assess the effect of different sources on O 3 production capacity, we calculated P O3 and OFP using the OFPT method and MIR method, respectively. The latest MIR values are used to calculate the OFP 49 . Figure 4e, f shows the contribution of each VOC source to the O 3 production rate according to the OFPT method. The largest source of the O 3 production rate was industrial combustion (2.6 ± 2.2 ppbv h −1 , 39.9 ± 3.7%) during CNY, followed by vehicle exhaust (1.3 ± 1.3 ppbv h −1 , 20.6 ± 3.6%), biomass burning/fireworks (0.9 ± 0.8 ppbv h −1 , 14.7 ± 0.6%), solvent use (0.9 ± 0.8 ppbv h −1 , 14.3 ± 1.0%), cooking (0.6 ± 0.5 ppbv h −1 , 8.6 ± 1.6%), and pharmaceutical exhaust (0.1 ± 0.1 ppbv h −1 , 1.9 ± 1.0%). For the MIR method, as shown in Fig. 4e and Supplementary Fig. 11c, industrial combustion contributed to 33.0 ± 14.3 ppbv (40.5 ± 3.7%) of the OFP during CNY, and the second important source contributing to the OFP was vehicle exhaust (14.5 ± 8.8 ppbv, 17.3 ± 2.9%), followed by solvent use (13.6 ± 5.5 ppbv, 17.1 ± 1.1%), biomass burning/fireworks (11.2 ± 5.1 ppbv, 13.7 ± 0.5%), cooking (7.0 ± 3.0 ppbv, 8.8 ± 1.3%), and pharmaceutical exhaust (2.1 ± 1.2 ppbv, 2.6 ± 1.0%). As shown in Fig. 4b, f, VOCs and OFPs had similar daily variation patterns, e.g., low values presenting at noon, because the OFPs were obtained by multiplying the VOC concentration and the MIR coefficient, which is a constant value under the specific condition being most sensitive to VOCs 49-51 . Thus, the OFP calculated with the MIR coefficient cannot reflect the actual atmospheric O 3 production in terms of diurnal variations, although the order of the source contribution to O 3 formation is similar to that of P O3 (Supplementary Figs. 11c, d or  Supplementary Table 2). Meanwhile, the P O3 showed a similar diurnal pattern with O 3 concentration in the daytime, but different from that of VOCs or OFP, as shown in Fig. 4d. The source apportionment based on the O 3 production rate should be more accurate than that based on the MIR method due to the following reasons. (1) The O 3 production rates are calculated based on the observed chemicals and meteorological parameters with a hightime resolution rather than that at fixed values of MIR because the formation of O 3 is highly correlated with free radicals in the daytime 44,52,53 . (2) The mean change rate of measured O 3 (dc O3, Measured /dt) and net P O3 (have considered the O 3 loss) are comparable during the observation period with the peak values of 4.3 and 5.1 ppbv h −1 at 10 a.m., respectively, and the OFP is obviously different from the actual situation ( Supplementary Fig.  13). (3) The P O3 is better to indicate the O 3 pollution than the OFP during the whole observation period. For example, similar to the O 3 concentrations, high P O3 values only occurred during CNY, while high OFP values appeared in all three periods (Figs. 1 and 4). (4) The P O3 generally shows similar diurnal variation patterns with O 3 during the observation period, while the OFP shows opposite variation patterns with O 3 (Fig. 4 and Supplementary Fig. 8). Therefore, the P O3 is able to correctly reflect the actual diurnal variations of each source to O 3 formation.
It should be noted that the total P O3 during CNY was 3.7 times higher than before and after CNY. Besides being the largest source of the P O3 (2.6 ± 2.2 ppbv h −1 ) during CNY, industrial combustion showed the biggest increment (4.8 times) regarding the P O3 during CNY when compared with the mean value before and after CNY. This is consistent with the fact that the air masses during CNY were dominantly from the southern and eastern regions in the local and neighboring areas (Fig. 2), featured by intensive industrial emissions. In addition, the RIR values of high reactive species (ethene, propene, toluene, benzene, and xylene) from industrial combustion accounted for 36.5% of the total RIR value during CNY (Supplementary Fig. 14). These results highlight the importance of industrial emissions in O 3 formation in Shijiazhuang, and highly reactive species play an important role under the conditions of weak irradiation and low temperature in winter. The different emission reduction scenarios of industrial sources were also further studied as shown in Supplementary Fig. 15. It can be seen that an 80-100% reduction of industrial emissions during CNY is needed to return to the O 3 levels before CNY, which highlights the importance of industry in the observed O 3 pollution event. Meanwhile, the P O3 attributed to vehicle exhaust and biomass burning/fireworks with large absolute values (1.3 and 0.9 ppbv h −1 ) also showed obvious increases (3.7-3.8 times) during CNY. Thus, attention should also be paid to these local sources of O 3 pollution. Although the VOC concentrations in summer ( Supplementary Fig.  12d) were lower than that in winter, the O 3 production rates in summer (13.1 ± 10.0 ppbv h −1 ) were larger than that in winter (2.9 ± 3.8 ppbv h −1 ) due to the stronger solar radiation. Meanwhile, the emission sources of VOCs varied obviously between winter and summer. Thus, the dominant source of O 3 also changed obviously. For example, the summertime sources of O 3 follows the order of: vehicle exhaust (3.3 ± 2.7 ppbv h −1 , 25.0 ± 7.7%) > biogenic (2.7 ± 2.7 ppbv h −1 , 20.4 ± 18.6%) > industrial combustion (2.4 ± 1.9 ppbv h −1 , 18.4 ± 3.7%) > cooking (1.6 ± 1.3 ppbv h −1 , 11.9 ± 3.6%) > solvent use (1.5 ± 1.4 ppbv h −1 , 11.8 ± 8.2%) > biomass burning (1.3 ± 1.1 ppbv h −1 , 10.3 ± 3.6%) > pharmaceutical exhaust (0.3 ± 0.3 ppbv h −1 , 2.3 ± 1.3%). In particular, biogenic VOCs only contributed 8.1% to the total VOCs in summer, while it was the second largest source of the P O3 . Industrial combustion, as the observed largest O 3 source in winter, became the third source of O 3 in summer. In addition, to understand the influence of the components of VOCs on O 3 production, a cross-simulation of meteorological conditions was performed. As shown in Supplementary Table 3, the P O3 of industrial combustion increased from 2.6 to 4.1 ppbv h −1 when using summer meteorological conditions during CNY. The virtual P O3 (4.1 ± 3.6 ppbv h −1 ) is even higher than the P O3 in summer (2.4 ± 1.9 ppbv h −1 ). The relative contribution of alkenes to P O3 in winter (52.8%) was also higher than that in summer (43.1%) as shown in Supplementary Fig. 16. These results highlight the highly reactive VOCs from industrial combustion in wintertime O 3 photochemistry, while they might be overlooked in summer due to chemical loss.

DISCUSSION
Although the solar radiation was much weaker in wintertime compared to that in summer, wintertime O 3 pollution events still occurred in SJZ during our observations. The weak solar radiation facilitates the measurement of highly reactive VOCs. The high RIRs and concentrations of alkenes highlight their importance to O 3 formation in winter. In summer, the contribution of these highly reactive species might be underestimated due to photochemical loss. The potential source contribution function (PSCF) results showed the vital role of photochemistry in the observed O 3 pollution event. What's more, O 3 pollution during CNY was strongly affected by the changes in weather systems, e.g., the abnormal weakening of the West Pacific subtropical high-pressure system in the southern regions and the temporary appearance of the cut-off low-pressure system in the northern regions of the observation site, which provided favorable conditions for the accumulation and reaction of pollution during CNY. Further analysis in combination with the trajectories of air masses reveals that the levels of VOCs and O 3 were well correlated with the ratio of short-ranged trajectories from the southern and eastern regions, implying that local and neighboring emissions contributed to the occurrence of O 3 pollution events during CNY. Highly reactive species such as ethylene, propene, toluene, and xylene with high RIR values positively correlated with O 3 production during CNY.
To accurately quantify the contribution of different VOC species to ozone production, the OPFT method was developed by tracing the O 3 production rate in different generations of VOC oxidation in the OBM, which paves the way for source apportionment of the O 3 production rate. The OPFT method reflects the dynamic O 3 production under atmospheric conditions at an observational site unlike the O 3 source apportionment method based on the MIR values at a fixed atmospheric scenario. In addition, the OPFT method does not interrupt the ratio of VOCs to NO x , which is a problem for the O 3 production rate source apportionment based on the RORE or RIR method 35 , during box model simulations. Thus, the OPFT provides more refined and accurate information on O 3 sources when compared with these previous studies. Our results also confirmed that the sequence of O 3 source contribution based on the OFPT method is credible when compared with that based on the MIR method if O 3 production is in a VOC-limited regime. Moreover, the source apportionment of the P O3 based on the OFPT method could capture not only the diurnal variations of O 3 formation but also the long-term time series of O 3 pollution in the ambient environment.
With the aid of the developed OFPT approach, the wintertime O 3 pollution event was analyzed in detail. The OFPT results found that the increased O 3 concentrations during CNY could be ascribed to enhanced O 3 production related to industrial combustion, vehicle exhaust, and biomass burning/fireworks, from which highly reactive VOCs were emitted during the CNY. A comparison with the summer results further confirms that highly reactive VOCs were visible in winter due to the lower radiation. It should be pointed out that the OFPT method in this study does not account for vertical or horizontal transport because it is based on a box model. A combination with a 3D air quality model should provide a more comprehensive understanding of O 3 pollution in the future.   40,56 . In addition, the reaction between NO 2 and OH leads to a net loss of O 3 in the daytime, thus it is also considered in Eq. 3 56 .

Field measurements
where k i means the corresponding reaction rate; [i] is the concentration of the species i. Figure 5 shows the workflow of the OFPT method to obtain the P O3 of each VOC. Briefly, for a VOC j , RO 2 , and HO 2 are searched in all the first-generation products after the reactions are initialized by oxidants (OH, NO 3 , and O 3 ). If RO 2 or HO 2 is found, the corresponding reaction rates of HO 2 /RO 2 with NO are recorded. Subsequently, all the first-generation products will further react with oxidants. RO 2 and HO 2 will be searched in all second-generation products again. If RO 2 or HO 2 is present, a constraint factor needs to be determined by performing a source-sink analysis of its precursors. The cause is that the traced RO 2 or HO 2 can be generated from both the first-generation products of the VOC j and other VOCs as well as their products. Then, the factor is used to constrain the reaction rate of the searched RO 2 or HO 2 with NO. The third-generation reactions are treated in the same way. The reaction rates for the fourthgeneration products are close to zero, thus we only account for three generations of RO 2 or HO 2 . In addition, considering that not all RO 2 and HO 2 radicals eventually generate O 3 due to homogeneous or non-homogeneous losses 57 , the P O3, HO2 and P O3, RO2 obtained from the source-sink equilibrium of O 3 are used to constrain the P O3, HO2 and P O3, RO2 from the VOC j , respectively. Moreover, besides oxidants, photolysis is also considered in the OFTP method.
VOC source apportionments were calculated by the PMF model (US EPA 5.0). In this study, forty-two species were selected for PMF analysis according to the principles as shown in Supplementary Note 2. The selected VOCs accounted for 94.3 ± 6.9% of the total observed VOC concentrations with an R 2 of 0.7, and 82.7 ± 11.9% of the total P O3 during the observation period. More information about PMF model setup and validation can be seen in Supplementary Note 2. A six-factor solution for wintertime VOCs was identified in this study. The contribution of each VOC source to the P O3 at a given time was further calculated by multiplying the fraction of the VOC attributed to a source and its corresponding P O3 traced in the OBM, i.e.,

OFPT output
The novelty of the OFPT approach where S PO3,i is the source contribution of the ith VOC source including n VOC species to O 3 production rate at a given time (ppbv h −1 ); F VOCj,i is the relative contribution of a VOC j to the ith VOC source (%), and P O3, VOCj is the O 3 production rate of a VOC j at a given time (ppbv h −1 ) according to the corresponding product rate of RO 2 and HO 2 produced from the VOC j .
Empirical kinetic modeling approach and RIR calculation The EKMA is used to describe the nonlinear relationship between O 3 and its precursors NO x and VOCs 44 3 formation is promoted by HONO, which is an important source of OH radicals 58 . Thus, HONO was accounted for in the model simulations in our study.
To identify highly reactive species, the RIR method was used to evaluate the sensitivity of different precursors to O 3 production 59 . Briefly, we changed the model input for a target precursor to a certain extent and obtained the correspondingly relative changes of P O3 compared to that of the target species, i.e., RIRðXÞ ¼ ΔPO 3 ðXÞ PO 3 ðXÞ ΔcðXÞ cðXÞ (5) where c(X) and Δc(X) are the measured concentration and the concentration changes of a precursor X (ppbv), respectively; P O3 and ΔP O3 represent the simulated P O3 in the base scenario and the changes of P O3 resulting from the concentration changes of the precursor. In this study, a relative change of 20% was used as the same as that in a previous study 52 . The O 3 concentrations were calculated in the daytime (8:00-18:00 local time).

Potential source contribution function analysis
The PSCF analysis was performed to understand the possible spatial distribution of the emission sources, by which the air mass transport trajectories at the observed site are calculated 60 . The 24 h backward trajectories of air mass were calculated with a 1-h temporal resolution from a global database using MeteoInfo software 61 . The area was divided into 0.5°× 0.5°grid to study VOCs and O 3 potential sources. Furthermore, a cluster analysis method was also performed to study the variation of the air mass trajectories during the observation period.

DATA AVAILABILITY
The data of air mass backward trajectories are available at ftp://arlftp.arlhq.noaa.gov/. The data of geopotential height and the wind is available at http:// www.cpc.ncep.noaa.gov/. Solar radiation (SR) was obtained from Copernicus Services (www.copernicus.eu/en). All additional data to perform the analyses are available upon reasonable request from the corresponding author (liuyc@buct.edu.cn).