Characterization of VOCs during Nonheating and Heating Periods in the Typical Suburban Area of Beijing, China: Sources and Health Assessment

: In recent years, the “coal to electricity” project (CTEP) using clean energy instead of coal for heating has been implemented by Beijing government to cope with air pollution. However, VOC pollution after CTEP was rarely studied in suburbs of Beijing. To ﬁll this exigency, 116 volatile organic compounds (VOCs) were observed during nonheating (P1) and heating (P2) periods in suburban Beijing. The results showed that the total of VOCs (TVOCs) was positively correlated with PM 2.5 , PM 10 , NO 2 , CO, and SO 2 but negatively correlated with O 3 and wind speed. The average TVOCs concentration was 19.43 ± 12.41 ppbv in P1 and 16.25 ± 8.01 ppbv in P2. Aromatics and oxygenated VOCs (OVOCs) were the main contributors to ozone formation potential (OFP). Seven sources of VOCs identiﬁed by the positive matrix factorization (PMF) model were industrial source, coal combustion, fuel evaporation, gasoline vehicle exhaust, diesel vehicle exhaust, background and biogenic sources, and solvent usage. The contribution of coal combustion to VOCs increased signiﬁcantly during P2, whereas industrial sources, fuel evaporation, and solvent usage exhibited opposite trends. The potential source contribution function (PSCF) and concentration weighted trajectory (CWT) were used to analyze the source distributions. The results showed that VOC pollution was caused mainly by air mass from southern Hebei during P1 but by local emissions during P2. Therefore, although the contribution of coal combustion after heating increased, TVOCs concentration during P2 was lower than that during P1. Chronic noncarcinogenic risks of all selected VOC species were below the safe level, while the carcinogenic risks of most selected VOC species were above the acceptable risk level, especially for tetrachloromethane and 1,2-dichloroethane. The cancer risks posed by gasoline vehicle emissions, industrial enterprises, and coal combustion should be paid more attention.


Introduction
Volatile organic compounds (VOCs) are defined as organic compounds with boiling points below 260 • C or saturated vapor pressures higher than 13.33 Pa at standard atmospheric pressure [1,2]. VOCs can form secondary aerosols [3,4] and enhance regional and even global atmospheric oxidizability as important precursors of tropospheric ozone (O 3 ) [5,6]. Some VOCs, such as polycyclic aromatic hydrocarbons (PAHs), oxygenated VOCs (OVOCs), and halohydrocarbons, have teratogenic and carcinogenic effects [7][8][9]. Additionally, tropospheric O 3 generated by photochemical reactions of VOCs leads to a decline in the CO 2 uptake of plants and crop yield [10,11].
Beijing, one of the world's megacities, with a population of more than 21 million and more than 5 million vehicles, has been suffering from severe photochemical pollution in recent decades [12]. Several studies have indicated that vehicle exhaust, solvent usage, fuel evaporation, and industrial sources were important sources of VOCs in Beijing [13,14] and that the species with high ozone formation potential (OFP) contribution were alkenes and aromatics [15]. Halohydrocarbons such as 1,2-dichloropropane, 1,2-dichloroethane, and chloroform might pose potential carcinogenic risks to the exposed populations in pollution episodes in urban Beijing [16]. Ambient total VOC (TVOC) concentration in Beijing shows significant seasonal variation, with the highest concentration during the heating period (15 November to 15 March of the next year) [17,18], which has been confirmed by satellitederived emission inventories [19]. This phenomenon is caused mainly by scattered coal combustion in suburbs of Beijing [20,21]. According to the announcement of the People's Government of Beijing Municipality in 2019, 850,000 rural households have completed the "coal to electricity" project (CTEP). According to the urban development plan of Beijing, suburbs will accelerate the task of urban population and industry transfer. However, most of the research on VOC pollution has been performed in urban Beijing. Few related studies have been conducted in suburban areas of Beijing, especially after CTEP.
The heating time of Beijing is generally from 15 November to 15 March of the next year. To illustrate the effect of heating after CTEP on VOC pollution in the suburbs of Beijing, the concentrations and compositions of ambient VOCs were measured during nonheating (P1, 5-14 November 2020) and heating (P2, 15-26 November 2020) periods in suburban Beijing. The maximum incremental reactivity (MIR) method was used to estimate the contribution of VOCs to O 3 formation. The major VOC sources and their temporal characteristics were identified by the positive matrix factorization (PMF) model. The potential source contribution function (PSCF) and concentration weighted trajectory (CWT) model was used to study the potential VOC source regions around Beijing. Finally, Monte Carlo simulation and the point estimate approach were used to assess residential inhalation noncancer and cancer risks.

Sampling Site
As the main implementation area of CTEP and the transfer of urban population and industry, the suburbs of Beijing lack research on VOC pollution. Therefore, this campaign was carried out at the Yanqi Lake campus of the University of Chinese Academy of Sciences (UCAS, 40 • 40 83 N, 116 • 68 45 E, 20 m above ground level). As shown in Figure S1, the sampling site was approximately 60 km north of downtown of Beijing, without large-scale pollution sources within 3 km. Therefore, the condition of this site was representative of the suburban VOC pollution level in Beijing. The concentrations of NO 2 , O 3 , CO, SO 2 , PM 2.5 , and PM 10 were obtained from https://quotsoft.net/air/ (accessed on 1 July 2021), and meteorological parameters (temperature (T), wind speed (WS), and wind direction (WD)) were provided by Yanshan Earth's Critical Zone and Surface Flux Observation Station.
Leakage tests and cleaning were performed before the summa canisters were used. As a blank sample, 200 mL of N 2 was repeatedly introduced to check whether the equipment pipeline had been contaminated until the signal intensities of 116 VOC species were below their detection limits. Then, a standard curve was drawn by seven concentration gradients (0.5 ppbv, 1 ppbv, 2 ppbv, 4 ppbv, 6 ppbv, 8 ppbv, and 10 ppbv). The method detection limit (MDL), correlation coefficient (R 2 ), and relative standard deviation (RSD) of VOC species in this study are listed in Table S1. During the sampling period, 0.5 ppbv standard gas was injected into the instrument to calibrate the retention time performance at 00:00 each day. Since the deviation of VOCs did not exceed 30%, the standard curve was not recalibrated with multiple points. N 2 , standard gas, and internal standard gas were all prepared using a dynamic gas diluter (4700, Entech Inc., Los Angeles, CA, USA).

Ozone Formation Potential (OFP) Calculation
O 3 in the troposphere is mainly formed photochemically from the interactions of VOCs and NOx in the presence of sunlight [22,23]. In this study, the MIR method was used to calculate the OFP of VOCs [24], and the calculation formula was as follows: where OFP i represents the O 3 formation potential of species i, VOC i represents the mass concentration of species i (µg·m −3 ), and MIR i represents the generation coefficient of species i in the maximum incremental reaction of O 3 .

Positive Matrix Factorization (PMF)
Factors contributing to atmospheric VOC concentration were determined using the U.S. Environmental Protection Agency (EPA) source apportionment PMF model, which has been widely used in previous studies [25,26]. This model assumes that the concentrations at receptor sites are impacted by linear combinations of source emissions, which are then derived as factors in the model. A PMF source apportionment analysis does not require source profiles as model inputs but does require an assumption about the number of factors. The source contribution and source profile can be solved by the least square method [27]. The principle of the PMF model is as follows: where x ij is the concentration of species j in sample i (ppbv), p is the number of pollution sources, g ik is the contribution of the emission source k to sample i (%), f kj is the concentration of species j from the emission source k (ppbv), and e ij is the residual.
Factor contribution and distribution were obtained by the minimizing objective function Q with g ik ≥ 0 and f kj ≥ 0 constraints, which was calculated as follows: where n is the number of samples, m is the number of species, and u ij is the uncertainty of species j in sample i. Equations (4) and (5) were utilized to calculate the uncertainty (Unc): where MDL is the method detection limit of species j and EF ij is the error fraction of species j in sample i. In this study, MDL values were obtained by repeatedly injecting 0.2 ppbv standard samples seven times, calculating the mean standard deviation, and multiplying by the statistical coefficient t (3.14) at the 99% confidence interval. The EF of active species such as alkenes, trimethylbenzenes, aldehydes, and ketones were set to 30%, and the EF of other species were set to 10%. The selection of VOC species in this study was based on the following conditions: (1) retaining important VOC species as the marker of sources; (2) removing VOC species with missing values; (3) removing VOC species with more than 25% of data below MDL; (4) excluding species with highly photochemical reactivity [28]; (5) further filtering the remaining VOC according to signal-to-noise ratio (S/N) and excluding species with S/N less than 0.5 [29]. Finally, 39 VOC species (listed in Table S2) were input in the PMF model. The PSCF was used to analyze the global meteorological data provided by the National Oceanic and Atmospheric Administration (NOAA) and calculate a 48 h backward trajectory from a 500 m altitude with the TrajStat software (MeteoInfo v2.3.1) The trajectory times were set at 1:00, 4:00, 7:00, 10:00, 13:00, 16:00, 19:00 and 22:00 in Beijing time (UTC + 8), based on the summa canisters' sampling time. The area covered by the backward trajectory was divided into 0.5 • × 0.5 • grid cells, and potential sources were identified and located in space [30,31]. The PSCF is defined as follows: where i represents latitude, j represents longitude, n ij represents the total number of trajectories passing through the ijth cell grid, and m ij represents the number of trajectories passing through the ijth cell grid when measured value exceeds the mean concentration of VOCs (17.60 ppbv) during the sampling period. The PSCF value is usually multiplied by a weight function W ij to reduce uncertainty [32]. W ij was defined per the following equation:

Concentration Weighted Trajectory (CWT)
Compared with PSCF, CWT combines pollutant concentration with residence time to calculate the contribution of different grids to the study area, which can more effectively distinguish medium from strong sources [33]. Therefore, it is necessary to refer to the analysis result of CWT, which was calculated as follows: where C l is the pollutant concentration corresponding to trajectory l in the ijth cell, τ ijl is the residence time of trajectory l in the ijth cell, and M is the total number of backward trajectories. The CWT value was corrected by the same weight function W ij as used with the PSCF method to reduce uncertainty [34,35].

Health Risk Assessment
Among three pathways of exposure to pollutants (dermal absorption, ingestion, and inhalation), this study assessed the health risks of toxic VOCs through the inhalation exposure pathway because of VOCs' high vapor pressure [36]. Among all measured species in this study, only 35 VOC species with dose-effect relationships were specified by the California Environmental Protection Agency (CalEPA), including 29 noncarcinogenic species and 20 carcinogenic species. The noncancer hazard quotient (HQ) and the hazard index (HI) for target organ systems were calculated as follows [37]: where HQ i is the noncancer inhalation hazard quotient of species i, C i is the average concentration of species i (µg/m 3 ), REL i is the reference exposure level for chronic inhalation of species i (Table S3), and HI is the sum of the HQs of multiple VOC species for the target organ system (Table S4). When HQ i increases to above 1.0, the probability of adverse health effects on the target organ system increases by an undefined amount [37].
The following algorithms were used to calculate residential inhalation cancer risk (Risk inh ): where Dose air is the residential inhalation dose (mg/kg-day) for carcinogenic VOCs, BR/BW is the daily breathing rate for normalized body weight (L/kg body weight-day); A is the inhalation absorption factor (unitless); EF is the exposure frequency (unitless); Risk inh is the residential lifetime inhalation cancer risk (unitless); CPF is the inhalation cancer potency factor (kg-day/mg); ASF and ED are the age sensitivity factor (unitless) and exposure duration (year) for a specific age group, respectively; AT is the average time (year) for lifetime cancer risk; and FAH is the fraction of exposure time (unitless). All the recommended default values for Equations (8) and (9) can be found in Tables S3 and S5. In this study, stochastic exposure assessment (Monte Carlo simulation) was applied to evaluate the noncarcinogenic and carcinogenic risks of specific VOC species. In this approach, according to the BatchFit tool (i.e., Anderson-Darling, chi-square, and Kolmogorov-Smirnov tests) of the CrystalBall ® Oracle software (v11.1.24), the optimal concentration distribution model of VOCs was fitted to obtain cancer risk caused by specific VOC species. This study performed 10,000 iterations in each independent run to ensure the stability and convergence of the simulation results. However, for residential carcinogenic risks caused by VOC sources, only deterministic assessment (point estimate approach) was conducted, because VOCs have different concentration fraction models and CPF values.

Overview of Meteorological Parameters and Air Quality
The hourly variation of meteorological parameters and air pollutants from 5-26 November 2020 is shown in Figure 1. During the sampling period, the mean values of T and WS were 10.66 • C and 3.03 m/s, respectively. As shown in the wind rose diagram in Figure S2, the dominant WDs during P1 were westerly and southeasterly, and the WD during P2 was mainly southeasterly. Additionally, the WS in P2 (3.39 ± 2.65 m/s) was significantly higher than that in P1 (2.58 ± 1.97 m/s). The average concentrations of NO 2   During the whole campaign, the Pearson correlation coefficient between TVOCs and PM 2.5 was the highest (r = 0.83, p < 0.01), indicating that they may have the same sources. In addition, TVOCs were positively correlated with PM 10 , NO 2 , CO, and SO 2 , and their correlation coefficients were 0.64, 0.68, 0.66 and 0.33 (p < 0.01), respectively. Furthermore, this study used trimethylbenzene (1,2,3-trimethylbenzene, 1,2,4-trimethylbenzene, 1,3, 5-trimethylbenzene) and benzene, which have a 30-fold difference in OH reaction rate and similar sources, to evaluate the air aging level. As shown in Figure 1e, the ratio of benzene and trimethylbenzene (B/T ratio) was positively correlated with TVOCs (r = 0.44, p < 0.01) and higher in heavily polluted weather, indicating that the consumption of trimethylbenzene was higher than that of benzene under regional transmission. TVOCs were negatively correlated with wind speed, indicating that stagnant weather with low wind speed was more conducive to accumulation of VOCs. The RO 2 generated by VOCs can oxidize NO to NO 2 , thereby reducing the consumption of NO to O 3 and finally increasing the O 3 concentration [38,39]. However, O 3 was negatively correlated with TVOCs and NO 2 in this study, and the occurrence of the peak concentration of O 3 was about 12 h later than that of the peak concentrations of TVOCs and NO 2 . This phenomenon might be due to heavy-pollution weather increasing NO concentration and maintaining a NO-O 3 -NO 2 dynamic equilibrium.
The diurnal variations in VOC groups, PM 2.5 , NO 2 , O 3 , and meteorological parameters are shown in Figure S3. The diurnal variations of alkanes, alkenes, acetylene, and aromatics were similar, with the lowest concentrations at 11:00-14:00 and higher concentrations at 7:00 and 19:00, which symbolized the bimodal features of vehicle emissions. Their diurnal variation was similar to that of NO 2 and PM 2.5 , but completely opposite to that of O 3 , because stronger solar radiation and higher temperatures at noon promoted photochemical reaction. OVOCs had higher concentrations at noon as well and had another peak at 4:00 due to heavy pollution at night. Moreover, atmospheric pollutants could be easier to diffuse under the conditions of higher wind speed and atmospheric boundary layer at 7:00-16:00.

Concentrations and Compositions of VOCs
The descriptive statistical results for 116 VOC species during different periods are summarized in Table S6. The averaged TVOCs concentration in P1 was 19.43 ± 12.41 ppbv, higher than that in P2 (16.25 ± 8.01 ppbv). Figure 2 shows average concentrations of VOC groups, their compositions, and their relative differences between P1 and P2. Most alkanes, aromatics, and OVOC species were more abundant during P1 than P2, and their relative differences were 17.21%, 30.26%, and 17.40%, respectively. The concentration of acetylene in P2 was 4.52% higher than that in P1, which may be related to scattered coal burning in the heating period. The concentrations of alkenes with high reactivity and halohydrocarbons representing background sources were similar in P1 and P2. These results indicated that the relative differences between P1 and P2 were impacted mainly by aged and regional air mass transport.

OFP of VOCs
The average OFP value during P1 was 119.49 µg/m 3 , which was significantly higher than that during P2 (92.63 µg/m 3 ). The top fifteen OFPs of VOC species during P1 and P2 are shown in Figure 3a [40,41]. Controlling emissions of related sources is an effective measure to reduce near-surface O 3 pollution in the suburbs of Beijing. Although alkanes were dominant in terms of concentration (35.32%), their contribution to OFP was only nearly 12.01%, indicating that alkanes were not significant species for O 3 formation. The contributions of acetylene and halohydrocarbons to OFP were almost negligible.

Source Apportionment
In the running of the PMF model, numbers of factors ranging from 4 to 10 were evaluated to determine the optimum solution. The variation of Q value with number of factors is shown in Figure S4. When the factor number was 7 to 10, Q(True)/Q(Robust) converged to 1.09. Theoretically, when the ratio of Q(true)/Q(robust) tends to 1, the impact of data points with high-scaled residuals is least. Based on the physical plausibility of the resulting profiles and residual distributions, the use of seven factors was determined to be the optimal solution. Values of Fpeak ranging from −1 to 1 at intervals of 0.2 were used in the model. An Fpeak of −0.2 was determined to be the optimum solution, and source profiles of the resolved factors are shown in Figure 4. Factor 1 was identified as industrial sources, characterized by high percentages (73.46%) of CS 2 , which is mainly used in the production of viscose, rubber, and rayon [42,43]. This factor also had high contributions to 1,1-dichloroethane (25.52%), 1,2-dichloroethane (20.94%), and benzene (30.38%). Such high loadings of halohydrocarbons and aromatics were caused by industrial emissions [44,45]. 2-butanone is used as solvent in the manufacturing industry, and aldehyde is the main pollutant generated by the paint and textile furniture manufacturing industry [46,47]. This result was in good agreement with the contribution of factor 1 to butanal (22.67%) and 2-butanone (19.35%).
Factor 7 possessed the property of high content of aromatics such as ethylbenzene (41.25%), m/p-xylene (50.67%), o-xylene (45.62%), and toluene (16.60%). The high ratio of aromatics was related to solvent usage in the adhesive, furniture, and shoemaking industries [40]. Additionally, 1,1-dichloroethane (19.78%) and 1,2-dichloroethane (22.86%) are commonly used as industrial adhesives in furniture, varnish, and paint stripper [64]. Hence, factor 7 was attributed to solvent usage.The time series of the seven factors, as well as their contribution percentages, are shown in Figure 5. Temporal variation in the concentration contribution of each source is shown in Figure 5a. All sources except for background and biogenic sources and diesel vehicle exhaust had peaks on 11 November, which were caused by persistent southern winds with low WS and temperature inversions. The contribution peaks of industrial sources and organic solvents appeared on several severely polluted days, and their Pearson correlation coefficient was 0.58, indicating that the two sources could be attributed to regional transport. In addition, the contribution peaks of fuel evaporation and gasoline vehicle exhaust occurred simultaneously, indicating that both factors were affected by traffic conditions. The contribution peak of diesel vehicle exhaust usually appeared at night, because trucks were allowed to enter downtown Beijing only at 23:00-6:00. The coal source had a high concentration after 14 November, in accordance with heating time. The background source contribution was relatively stable, with the smallest standard deviation, at 0.54, among all the sources (the average contribution of each source was set to 1).  Figure 5b shows that the average concentrations of sources during P1 and P2 were significantly different. The contribution of coal combustion during P2 increased significantly, accounting for 21.23% (3.45 ppbv) compared with only 11.46% (2.18 ppbv) during P1. This indicated that scattered coal might still be used in suburbs around Beijing. Contrarily, the relative contributions of industrial sources, fuel evaporation, and solvent usage all decreased more than 35% from P1 to P2. The contributions of background and biogenic sources and gasoline and diesel vehicle exhaust between P1 and P2 had no significant differences, indicating that heating activities had little impact on plant photosynthesis and traffic conditions.

PSCF and CWT Results
As shown in Figure 6a, the northwest long-distance transport from Russia, Mongolia, and the Inner Mongolia Autonomous Region contributed 60.22% to Beijing's air mass during P1. Contributions from short-distance transport from the southwest and mediumdistance transport from the northeast were relatively small, amounting to 27.27% and 12.50%, respectively. From the results of PSCF, the regions that featured high PSCF values were concentrated mainly in southern Hebei Province and Tianjin. Based on the CWT results in Figure 6b, the strong sources of VOCs during P1 were mainly Baoding, Hebei Province, and local sources in Beijing. As shown in Figure 6c, during P2, the main air mass came from the medium-distance transport from the western (30.21%) and northern Inner Mongolia Autonomous Region (27.08%), followed by long-distance transport from Mongolia (19.79%), Heilongjiang Province (14.58%), and medium-distance transport from Shandong Province (8.33%). The area with high PSCF value in P2 was smaller than that in P1, mainly in Beijing, northern Hebei Province, and Tianjin. As presented in Figure 6d, areas with CWT > 20 were located mainly in urban Beijing.
To the west of Beijing are the Taihang Mountains. The north and northeast are surrounded by the Yanshan Mountains, and the southern provinces of Hebei, Henan, and Shandong have many industrial enterprises [65]. It is easy for Beijing to accumulate air pollutants in the south wind and stagnant weather conditions, which was the main reason for the serious VOC pollution during P1. Nevertheless, air masses were mainly from the northwest. Therefore, the TVOCs concentration during P2 was lower than that during P1, although the contribution of coal combustion increased during the heating period.
3.6. Health Risk Assessment 3.6.1. Noncarcinogenic and Carcinogenic Risks of VOC Species The mean values of HQ of the 29 selected VOC species were all <1 during both periods. As shown in Figure S5, acrolein, benzene, and tetrachloromethane had relatively high HQ values during P1 and P2. The HI was calculated by summing the HQ values of VOC species affecting the same target organ (Table S7). The respiratory system had the highest HI values during P1 and P2, which were 0.83 and 0.65, respectively. Nevertheless, the HI values of all target organ systems were lower than 1, indicating that the chronic health risks caused by inhalation of selected VOCs in the suburbs of Beijing could be negligible.
For residential lifetime cancer risk assessment, the U.S. EPA divides the population into third trimester, infants (0-2), children (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16), and adults . The normalized respiration rate of children was the highest after calculation, which was caused by the combined effect of respiration rate, body weight, and daily exposure time. Hence, the Monte Carlo simulation in this study used the correlation coefficient of children. As shown in Figures 7 and 8 VOCs including benzene, tetrachloromethane, and 1,2-dichloroethane were higher than the U.S. EPA acceptable standard (1.0 × 10 −6 ) during P1 and P2, and benzene was even higher than the stricter acceptable standard of the Canadian EPA (1.0 × 10 −5 ) [66]. Furthermore, benzene and tetrachloromethane had high HQ values that might lead to noncarcinogenic hazards as well. Therefore, it is particularly important for suburban Beijing to control the VOC sources related to benzene and tetrachloromethane.

Carcinogenic Risks of VOC Sources
In order to propose more effective source control measures, we further calculated the carcinogenic risks of selected VOC species of PMF-derived sources by point estimation (Figure 8). The mean cancer risks of six sources ranged from 1.3 × 10 −6 (diesel vehicle exhaust during P2) to 1.40 × 10 −5 (gasoline during P1). The integrated carcinogenic risks during P1 and P2 were 4.30 × 10 −5 and 2.80 × 10 −5 , which were below the U.S. EPA tolerable risk level (1.0 × 10 −4 ). Out of six PMF-derived sources, gasoline vehicle exhaust, industrial source, and coal combustion exhibited higher risk values during both periods due to their high contributions to benzene or 1,2-dichloroethane. This was inconsistent with the finding that background sources had the highest cancer risks in some cities in Canada, which may be related to the local industrial structure [57,67]. In short, for public health, gasoline vehicle emissions, industrial enterprises, and coal combustion need more effective raw material pretreatment and emission control. Considering the regional transmission of VOCs, not only suburban Beijing but the surrounding provinces need to have corresponding control measures.

Conclusions
To assess the impact of heating on pollution characteristics and sources and health assessment in regard to VOCs, 116 VOC species were continuously monitored at UCAS, Beijing during 5-26 November 2020. The averaged TVOCs concentration during P1 was 19.43 ± 12.41 ppbv, higher than that (16.25 ± 8.01 ppbv) during P2. Most alkanes, aromatics, and OVOC species gave higher contributions during P1 than P2, and the concentrations of alkenes and halohydrocarbons were close during P1 and P2 because of the aged and regional air masses. The concentration of acetylene in P2 was 4.52% higher than that in P1, which may be related to scattered coal burning during the heating period.
During the whole sampling campaign, the concentration of TVOCs was positively correlated with the concentrations of PM 2.5 , PM 10 , NO 2 , CO, and SO 2 . However, O 3 was negatively correlated with TVOCs and NO 2 in this study, and the peak concentration of O 3 was about 12 h later than those of TVOCs and NO 2 . Alkanes contributed most to the VOC mixing ratio during both P1 and P2, while aromatics and OVOCs contributed most to the OFP, with contributions of 36.44% and 34.43% during P1 and P2, respectively. These results indicate that vehicle exhaust and solvent usage were the main sources of the precursors of ozone, while the plant emissions and secondary generation contributed more OFP during P2. Controlling VOC and NOx emissions of gasoline vehicle exhaust, solvent usage, and industrial enterprises would constitute an effective measure to reduce near-surface O 3 pollution in the suburbs of Beijing.
The 48 h backward trajectories, PSCF, and CWT showed that the air quality of Beijing during P1 was affected mainly by short-distance transportation from the southern Hebei Province. VOC pollution in Beijing was basically caused by local source emissions during P2, so the TVOCs concentration during P2 was lower than that during P1. After CTEP, although heating still contributed to VOCs in the suburbs of Beijing, the main cause of VOC heavy pollution events was regional transmission. Therefore, regional joint prevention and control of the 2 + 26 cities in North China is vital to Beijing's air quality.
The seven sources identified from the PMF model were industrial source, coal combustion, fuel evaporation, gasoline vehicle exhaust, diesel vehicle exhaust, background and biogenic sources, and solvent usage. The contribution of coal combustion increased significantly during P2, indicating that scattered coal burning should be strictly controlled in suburban Beijing. Industrial sources, fuel evaporation, and solvent usage exhibited the opposite trend because of the regional transmission from southern Hebei Province during P1. Background and biogenic sources and gasoline and diesel vehicle exhaust showed no significant variation. Regarding the health risk assessment, the chronic noncarcinogenic risks of all selected VOC species were within the safe range, while the carcinogenic risks of most selected VOCs were above the acceptable risk level. From the perspective of public health, the cancer risks posed by gasoline vehicle emissions, industrial enterprises, and coal combustion should be paid more attention.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/atmos13040560/s1, Table S1. List of MDL, R 2 and RSD of VOCs species; Table S2. Summary of PMF input data. Concentration units in ppbv ; Table S3. The reference exposure level (REL) and cancer potency factor (CPF) of selected VOC species; Table S4. The target organ systems affected by non-carcinogenic adverse effects of selected VOC species; Table S5. Exposure parameters for calculating residential inhalation cancer risk for each age group; Table S6. Statistics of VOCs species during P1 and P2; Table S7. Substance-specific chronic inhalation hazard quotients (HQs) and the hazard index (HI) by target organ system during P1 and P2; Figure S1.

Conflicts of Interest:
The authors declare no conflict of interest.