Using Stable Sulfur Isotope to Trace Sulfur Oxidation Pathways during the Winter of 2017–2019 in Tianjin, North China

After the implementation of the Coal Replacing Project (CRP) in the northern parts of China in 2017, its effect on PM2.5 composition is still unclear. In the study, water-soluble ionic components (WSICs) and stable sulfur isotope ratios (δ34S) of SO42− in PM2.5 collected during the domestic heating period before and after the implementation of CRP in Tianjin were analyzed. Results showed that the average concentrations of both PM2.5 and WSICs have dropped dramatically after the CRP, especially for the SO42− (by approximately 57–60%). After the CRP, the range of δ34Ssulfate was significantly narrowed to 4.1–7.5‰ in January 2018 and 1.4–6.1‰ in January 2019, which suggested that the sulfur source was becoming simple. It was interesting that the δ34Ssulfate value in the pollution period before the CRP was higher than that in the clean period, whereas it showed the opposite tendency after the CRP, which implied that the contribution of sea salt was high during the pollution period before the CRP. The MIXSIAR model calculated that the contributions of the transition-metal ion (TMI) oxidation and NO2 oxidation pathways in the three sampling stages were higher than those of the OH radical oxidation and H2O2/O3 oxidation pathways, indicating that the formation pathway of sulfate was mainly dominated by heterogeneous oxidation. Before the CRP, the NO2 oxidation pathway was the dominant sulfate oxidation pathway during a haze episode, and the TMI oxidation pathway dominated the formation of sulfates after the CRP.


Introduction
In recent years, many urban areas in China have been severely affected by haze [1][2][3]. The primary emissions, secondary formation and transformation, and stable weather conditions have been considered as the main reasons for such widespread and continuous high loading of air pollutants, including haze, over China [4][5][6]. During certain special events, the Chinese government adopted some temporary policies to pursue a great improvement in air quality, such as shutting down heavily polluting industries, prohibiting construction activities, and reducing the number of vehicles in operation [7,8]. However, the impact of these measures was limited, and the improvement in air quality could not be maintained for a long time [7,9,10]. Among many control measures taken during the "Action Plan on Prevention and Control of Air Pollution" since 2013, the Coal Replacing Project (CRP), i.e., the consumption of electricity and natural gas instead of coal for energy, was an important solution for the Chinese government to achieve the goal of reducing fine particulate matter (PM 2.5 ) permanently.
It has been widely reported that the CRP has significant impacts on the energy, the economy, and the environment domains [11][12][13]. Lin and Jia (2020) concluded that the main significance of the coal-to-electricity project was not saving energy but improving air quality [12]. Ren et al. (2017) reported that the contributions of coal combustion to n-alkanes, PAHs, and OPAHs in Urumqi have been decreased from 21~75% before the implementation of CRP to 4.0~21% after the CRP [14]. Feng et al. (2020) found that the contribution of coal combustion to nitrate in the winter has been decreased by about 8%  (Table S1). PM 2.5 samples collection membranes were pre-combusted when collected (450 • C for 6 h). Quartz filters (23 × 18 cm) using a high-volume air sampler (TE-PM2.5HVP-BL, TISCH, Cleves, OH, USA) at a flow rate of 1.05 m 3 min −1 during 8:00-18:00 (local time) for daytime and from 18:00 to next morning 8:00 for nighttime in both the seasons, and stored at −20 • C. The membranes were weighed before and after sampling using a microelectronic balance (Mettler Toledo ML 204T, Switzerland) after a 48-h equilibration in a drying vessel. One blank filter was also collected in each campaign. The samples were classified into clean (PM 2.5 ≤ 75 µg m −3 ) and polluted (PM 2.5 > 75 µg m −3 ) periods based on the PM 2.5 mass concentration. Temperature (T), relative humidity (RH), and wind speed (WS) were measured at sampling points using a meteorological station (Gill MetPak, Gill Instruments, Lymington, Hampshire, UK).

Chemical Composition Analysis
Two filter cuts with a diameter of 44 mm (or four filter cuts of 26 mm in diameter if the PM2.5 concentration was low) were used for the analysis of WSICs. The filters were extracted into 15 mL ultrapure water (18.2 MΩ cm, Millipore, Burlington, MA, USA) under ultrasonication for 30 min at room temperature. The extractants were filtered through a 0.22 μm PTFE membrane filter and then the ions were measured using ion chromatography (ICS-5000+, Thermal, Sunnyvale, CA, USA). Cations were measured using a Dionex IonPacTM CS12A analytical column and a Dionex CRDS 600 suppressor on a conduction detector, and 20 mM MSA was used as an eluent. At a flow rate of 1.0 mL min −1 for anion measurement, 30 mM KOH was used as an eluent with a Dionex IonPacTM AS11-HC RFIC TM analytical column. The detection limit for WSICs was lower than 2 ng mL −1 with an error < 5%.

Determination of Sulfur Isotopic (δ 34 S) Ratios
Half of each quartz fiber filter was cut into pieces for δ 34 S analysis. They were extracted into 50 mL of ultrapure water and sonicated for 15 min, and then repeated three times. The solution was then filtered with 0.45 μm syringe filters. The filtrate was acidified to a pH < 2 by adding HCl solution and then trapped as barium sulfate (BaSO4) by adding droplets of supersaturated BaCl2 solution. The BaSO4 precipitation was filtered using 0.22 μm acetate filters 24 h later. After heating in an oven with step-wise heating at 200 °C for 2 h, 400 °C for 2 h, and 850 °C for 4 h, 340-450 μg BaSO4 were weighted and wrapped in tin capsules. The tin cups were crushed into a small size and placed into the isotope ratio mass spectrometer (MAT253, Thermo, Waltham, MA, USA). The δ 34 S of BaSO4 was

Chemical Composition Analysis
Two filter cuts with a diameter of 44 mm (or four filter cuts of 26 mm in diameter if the PM 2.5 concentration was low) were used for the analysis of WSICs. The filters were extracted into 15 mL ultrapure water (18.2 MΩ cm, Millipore, Burlington, MA, USA) under ultrasonication for 30 min at room temperature. The extractants were filtered through a 0.22 µm PTFE membrane filter and then the ions were measured using ion chromatography (ICS-5000+, Thermal, Sunnyvale, CA, USA). Cations were measured using a Dionex IonPacTM CS12A analytical column and a Dionex CRDS 600 suppressor on a conduction detector, and 20 mM MSA was used as an eluent. At a flow rate of 1.0 mL min −1 for anion measurement, 30 mM KOH was used as an eluent with a Dionex IonPacTM AS11-HC RFIC TM analytical column. The detection limit for WSICs was lower than 2 ng mL −1 with an error < 5%.

Determination of Sulfur Isotopic (δ 34 S) Ratios
Half of each quartz fiber filter was cut into pieces for δ 34 S analysis. They were extracted into 50 mL of ultrapure water and sonicated for 15 min, and then repeated three times. The solution was then filtered with 0.45 µm syringe filters. The filtrate was acidified to a pH < 2 by adding HCl solution and then trapped as barium sulfate (BaSO 4 ) by adding droplets of supersaturated BaCl 2 solution. The BaSO 4 precipitation was filtered using 0.22 µm acetate filters 24 h later. After heating in an oven with step-wise heating at 200 • C for 2 h, 400 • C for 2 h, and 850 • C for 4 h, 340-450 µg BaSO 4 were weighted and wrapped in tin capsules. The tin cups were crushed into a small size and placed into the isotope ratio mass spectrometer (MAT253, Thermo, Waltham, MA, USA). The δ 34 S of BaSO 4 was standardized with laboratory barite standards (referenced to NBS-127 and IAEA , and reported to the international sulfur isotope standard V-CDT [22]. The δ 34 S of SO 4 2− was calculated by Equation (1): where the isotopic ratio of ( 34 S/ 32 S) standard = 0.044. Overall uncertainty for δ 34 S sulfate values reported here were ±0.2‰ [22].

MIXSIAR Model to Calculate the Formation Pathways of Sulfate Aerosols
In this estimation, isotope fractionation effects from SO 2 to SO 4 2− have to be considered. Many previous studies have shown that the fractionation factors of δ 34 S from SO 2 to SO 4 2− depend on the extent of the reaction [24,29]. Moreover, the isotope fractionation factor between SO 2 and SO 4 2− could be calculated by the Rayleigh distillation model, which was based on an instantaneous phase equilibrium process under an open system (Rayleigh 1896). In the model, the δ 34 S value of SO 2 is a function of δ 34 S emission , the fraction: where δ 34 S SO2 is the δ 34 S of remaining SO 2 in the atmosphere; ε obs is the fractionation factor between SO 2 and SO 4 2− ; SOR is the SO 2 oxidation ratio; δ 34 S emission is the δ 34 S of SO 2 in emission sources, which can be estimated by the following formula based on the isotope mass balance [21]: where δ 34 S sulfate is the δ 34 S value of observed SO 4 2− in the atmosphere. When SOR is close to 1, δ 34 S emission equals to δ 34 S sulfate . Assuming the emitted SO 2 is completely converted to sulfate (SOR = 1), δ 34 S emission is calculated according to the linear relationship between δ 34 S sulfate and SOR in each sampling stage. According to the Equations (2) and (3), the fractionation factor (ε obs ) can be calculated: Using the observed δ 34 S sulfate , SOR, and the estimated δ 34 S emission , we found the mean values of ε obs during the sampling period were 6.6 ± 2.2‰, 0.6 ± 1.1‰, −3.1 ± 0.4‰ and 4.2 ± 1.1 ‰ during January 2017, July 2017, January 2018, and January 2019, respectively. The difference in ε obs values reflected the different formation mechanisms of sulfate aerosols. The oxidation pathways of sulfate mainly involve gas-phase oxidation and aqueous-phase reactions [23][24][25]. Several laboratory experiments measured the fractionation factors (ε) for SO 2 oxidation through H 2 O 2 and O 3 , and these two pathways showed similar ε values ranging from +15.1‰ to +17.4‰ [28]. Therefore, we used ε H2O2/O3 to represent the combined isotopic effect of O 3 and H 2 O 2 pathways. The fractionation factors (ε) for SO 2 oxidation through OH radical, H 2 O 2 /O 3 , NO 2 , and TMI pathways can be calculated [28]: lnε NO2 (‰) = (0.2437 ± 0.0457)/T − (0.0008 ± 0.019) where T is the ambient temperature ( • C). Based on the observed temperature in Tianjin during the sampling period, the fractionation factors of different oxidation pathways were obtained. Instead, the ε obs should be a result of the mixing of multiple oxidation pathways [21]: where f i is the contribution of oxidation pathway i to sulfate formation, and f OH + f H2O2/O3 + f TMI + f NO2 = 1. Based on Equation (8), a Bayesian isotope mixing model (ran in the Stable Isotope Analysis in the R, MIXSIAR) was used to assess the contributions of oxidation pathways to sulfate. A logical prior distribution was established in this model first, and the contribution of each oxidation pathway to the mixture was then determined.

HYSPLIT Backward Trajectory Analysis
This study mainly used the Hybrid Single Particle Lagrangian Integrated Trajectory Model (HYSPLIT4) of the MeteoInfo software to analyze the trajectory of backward gas and backward air mass. In this study, 100 m was regarded as the average flow field height of atmospheric boundary layer, and the backward trajectory of air mass movement is simulated every 48 h. Combined with the Pearson correlation analysis of sulfur and RH, T, SOR, tried to clarify the main factors affecting the change of S concentration in PM 2.5 .

Boundary Layer Height (BLH) Simulation
In order to explore the influence of BLH change on PM 2.5 pollution in Tianjin during the pollution period, the study was conducted in the region of 116.75 • to 118 • east longitude and 38.5 • to 39.5 • north latitude, including most of Tianjin and a small part of Bohai Sea. The study area consisted of 30 grids (one grid corresponds to one height value), the resolution is 0.25 • × 0.25 • (about 27 km × 27 km), and the BLH per hour was obtained by connecting 30 grids (height value) together. The Data from "The Climate Data Store" (https://cds.climate.copernicus.eu/cdsapp#!/home, accessed on 5 May 2022).

Effect of CRP on PM 2.5 Loading and Its Composition
As summarized in Table 1, the average concentration of PM 2.5 in January 2018 and January 2019 dropped by approximately 60% and 57%, respectively, compared to that in January 2017. In addition, the mean value of boundary layer height (BLH) in 2017 was 313 m, which was close to 300 m in 2019 ( Figure S1), indicating a great improvement in air quality after the implementation of CRP in Tianjin. In the winter of three years, three typical haze episodes (Ep1, Ep2 and Ep3) during sampling were discussed in Section 3.3.
However, the PM 2.5 concentration two years after the implementation of CRP (January 2019) was slightly higher than that one year after the implementation of the project (January 2018), and the SO 2 and NO 2 concentrations in the atmosphere also showed similar annual changes (Table 1). In addition, the RH and WS were similar in January 2018 and January 2019 (Table 1). We speculated that, first, BLH in 2019 (300 m) was much smaller than that in 2018 (524 m), which makes it difficult to diffuse pollutants ( Figure S1). Secondly, we assumed that in the first year of the implementation of CRP, the "coal-to-gas" and "coal-to-electricity" transformation projects were implemented strictly, and the amount of coal consumption was greatly reduced, resulting in a very low PM 2.5 concentration in January 2018. Whereas in the second year of the CRP, some areas have the situation that the natural gas supply was not enough to support the heating demand and the consumption of coal was partly resumed to supplement the natural gas, and the coal consumption in January 2019 was slightly higher than that in January 2018 ( Figure S2), which should be one of the main reasons for causing the increased PM 2.5 loading in January 2019 compared to that in January 2018.
As shown in Figure 2, sulfate, nitrate, and ammonium (SNA) were the most important three components of WSICs, and the variation of SNA was consistent with the PM 2.5 concentrations. The concentration of SNA in January 2017 was higher than that in January 2018 and January 2019. Otherwise, as shown in Figure S3, the proportion of sulfate in SNA decreased, and the proportion of nitrate and ammonium in SNA increased after the CRP was implemented. However, the PM2.5 concentration two years after the implementation of CRP (January 2019) was slightly higher than that one year after the implementation of the project (January 2018), and the SO2 and NO2 concentrations in the atmosphere also showed similar annual changes (Table 1). In addition, the RH and WS were similar in January 2018 and January 2019 (Table 1). We speculated that, first, BLH in 2019 (300 m) was much smaller than that in 2018 (524 m), which makes it difficult to diffuse pollutants ( Figure S1). Secondly, we assumed that in the first year of the implementation of CRP, the "coal-togas" and "coal-to-electricity" transformation projects were implemented strictly, and the amount of coal consumption was greatly reduced, resulting in a very low PM2.5 concentration in January 2018. Whereas in the second year of the CRP, some areas have the situation that the natural gas supply was not enough to support the heating demand and the consumption of coal was partly resumed to supplement the natural gas, and the coal consumption in January 2019 was slightly higher than that in January 2018 ( Figure S2), which should be one of the main reasons for causing the increased PM2.5 loading in January 2019 compared to that in January 2018.
As shown in Figure 2, sulfate, nitrate, and ammonium (SNA) were the most important three components of WSICs, and the variation of SNA was consistent with the PM2.5 concentrations. The concentration of SNA in January 2017 was higher than that in January 2018 and January 2019. Otherwise, as shown in Figure S3, the proportion of sulfate in SNA decreased, and the proportion of nitrate and ammonium in SNA increased after the CRP was implemented.

Sulfur Isotope Composition of Sulfate Aerosols
As shown in Figure 3, before the implementation of the CRP (January-2017), the δ 34 S sulfate varied widely from 3.5 to 9.5‰ with a mean value of 5.5 ± 1.6‰, and the range of δ 34 S sulfate was narrowed after the CRP was implemented, which were 4.1~7.5‰ in January 2018 with a mean value of 5.6 ± 0.9‰, and 1.4~6.1‰ in January 2019 with a mean value of 4.4 ± 1.1‰. After the CRP, the range of δ 34 S sulfate value gradually became smaller and the fluctuation was not high (Figure 3), which showed that with the replacement of bulk coal, the sulfur source in the atmosphere might be more stable and simpler. Moreover, the δ 34 S sulfate in the pollution period (PM 2.5 ≥ 75 µg m −3 ) before the implementation of CRP was higher than that in the clean period (PM 2.5 < 75 µg m −3 ), while the δ 34 S sulfate in the pollution period after the implementation of CRP was lower than that in the clean period. Such a decreasing trend of δ 34 S sulfate from the clean period to the pollution period after the CRP might be interpreted by the changes in the contributions of various sources and formation pathways to sulfate under different atmospheric conditions.

Sulfur Isotope Composition of Sulfate Aerosols
As shown in Figure 3, before the implementation of the CRP (January-2017), the δ 34 Ssulfate varied widely from 3.5 to 9.5‰ with a mean value of 5.5 ± 1.6‰, and the range of δ 34 Ssulfate was narrowed after the CRP was implemented, which were 4.1~7.5‰ in January 2018 with a mean value of 5.6 ± 0.9‰, and 1.4~6.1‰ in January 2019 with a mean value of 4.4 ± 1.1‰. After the CRP, the range of δ 34 Ssulfate value gradually became smaller and the fluctuation was not high (Figure 3), which showed that with the replacement of bulk coal, the sulfur source in the atmosphere might be more stable and simpler. Moreover, the δ 34 Ssulfate in the pollution period (PM2.5 ≥ 75 μg m −3 ) before the implementation of CRP was higher than that in the clean period (PM2.5 < 75 μg m −3 ), while the δ 34 Ssulfate in the pollution period after the implementation of CRP was lower than that in the clean period. Such a decreasing trend of δ 34 Ssulfate from the clean period to the pollution period after the CRP might be interpreted by the changes in the contributions of various sources and formation pathways to sulfate under different atmospheric conditions.

Sulfur Isotopic Compositions of Sulfate during the Haze Episodes
To better explore the reasons for the change of δ 34 Ssulfate values in polluted days after the implementation of CRP, three typical haze episodes during sampling were discussed in this study. The first haze episode (Ep1) started on 1 February 2017 and lasted until 1 August 2017, with average PM2.5 concentrations of 217.5 μg m −3 and the highest value is 339.7 μg m −3 . The average δ 34 Ssulfate value is 5.9‰ during Ep1, which was 25.5% higher than that in the clean period. Furthermore, the measured δ 34 Ssulfate displayed significant positive correlations with temperature (T) (R = 0.86; p < 0.01). However, it was negatively correlated with relative humidity (RH) (R = −0.45; p < 0.01) and SO2 oxidation ratio (SOR = [SO4 2− ]/[SO4 2− + SO2]) (R = −0.60; p < 0.01) (Figure 4). The same negative correlation between the SOR and δ 34 Ssulfate value during the haze events observed in Nanjing, Beijing, and other places [21,29], indicating that the δ 34 Ssulfate values would be affected by different SO2 oxidation processes. In addition, the significant negative correlation between δ 34 Ssulfate and RH indicated that high RH promoted the isotope fractionation process in the water phase oxidation reaction of SO2 [29]. Combined with the backward trajectory analysis results ( Figure 5), it was found that the air mass during Ep1 mainly came from cluster 2 and cluster 3, which might be affected by the sea salt sources with high δ 34 Ssulfate values (21‰) ( Table S2) [30]. Therefore, before the implementation of CRP, the contribution of sea salt with high δ 34 Ssulfate might be one of the main reasons for the higher δ 34 Ssulfate in the pollution period than that in the clean period.

Sulfur Isotopic Compositions of Sulfate during the Haze Episodes
To better explore the reasons for the change of δ 34 S sulfate values in polluted days after the implementation of CRP, three typical haze episodes during sampling were discussed in this study. The first haze episode (Ep1) started on 1 February 2017 and lasted until 1 August 2017, with average PM 2.5 concentrations of 217.5 µg m −3 and the highest value is 339.7 µg m −3 . The average δ 34 S sulfate value is 5.9‰ during Ep1, which was 25.5% higher than that in the clean period. Furthermore, the measured δ 34 S sulfate displayed significant positive correlations with temperature (T) (R = 0.86; p < 0.01). However, it was negatively correlated with relative humidity (RH) (R = −0.45; p < 0.01) and SO 2 oxidation ratio (SOR = [SO 4 2− ]/[SO 4 2− + SO 2 ]) (R = −0.60; p < 0.01) (Figure 4). The same negative correlation between the SOR and δ 34 S sulfate value during the haze events observed in Nanjing, Beijing, and other places [21,29], indicating that the δ 34 S sulfate values would be affected by different SO 2 oxidation processes. In addition, the significant negative correlation between δ 34 S sulfate and RH indicated that high RH promoted the isotope fractionation process in the water phase oxidation reaction of SO 2 [29]. Combined with the backward trajectory analysis results (Figure 5), it was found that the air mass during Ep1 mainly came from cluster 2 and cluster 3, which might be affected by the sea salt sources with high δ 34 S sulfate values (21‰) (Table S2) [30]. Therefore, before the implementation of CRP, the contribution of sea salt with high δ 34 S sulfate might be one of the main reasons for the higher δ 34 S sulfate in the pollution period than that in the clean period.  The haze episode 2 (Ep2) from 12 January 2018 to 14 January 2018, with average PM 2.5 concentrations of 109.1 µg m −3 and the highest value was 155.5 µg m −3 . Compared with Ep1, the mass concentration of PM 2.5 in Ep2 decreased by 49.8%, which implies that the air quality of the Tianjin urban area has been significantly improved after the CRP. The average δ 34 S sulfate value of Ep2 is 5.1‰, which was 15.7% lower than Ep1. Moreover, it was noteworthy that the δ 34 S sulfate value of Ep2 was lower than that of the clean period during sampling, which was significantly different from the sampling period before the implementation of CRP (Figure 3b). In addition, the measured δ 34 S sulfate of Ep2 was positively correlated with T (R = 0.70, p < 0.01). Since the BLH of Ep2 was significantly higher than that of Ep1 and Ep3 ( Figure S1), the air mass had better diffusion ability, so the source and oxidation pathway of sulfate might change. Combined with the backward trajectory analysis (Figure 5), we found that the air mass during Ep2 mainly came from cluster 2 and cluster 4, which were mainly affected by the local sources and regional transportation in the southwest.
The haze episode 3 (Ep3) started on 10 January 2019 and lasted until 14 January 2019, with average PM 2.5 concentrations of 127.9 µg m −3 and the highest value was 145.7 µg m −3 . Compared with Ep2, the PM 2.5 concentration of Ep3 did not decrease significantly but increased slightly, which implies that the environmental effect of the replacement of bulk coal was limited, and sulfate was still mainly derived from coal during the pollution period. Therefore, additional environmental protection policies (i.e., desulfurization and denitrification of coal) should be implemented to improve the atmospheric environment in the future. After the implementation of CRP, the air masses during Ep2 mainly came from cluster2 and cluster4, while the air masses during Ep3 mainly came from cluster1 and cluster3 ( Figure 5). It was shown that the sources of air masses in the two episodes were similar, and they were mainly from local sources and regional transportation in the southwest ( Figure 5). However, it was noteworthy that the δ 34 S sulfate values in Ep3 decreased by 47.5% and 39.2% compared to Ep1 and Ep2, respectively. Since the regional transport of Ep2 and Ep3 is similar, it is speculated that the source has little influence on δ 34 S sulfate , and the variations of this value may be mainly due to isotopic fractionation during the SO 2 conversion to SO 4 2− after the CRP. The haze episode 2 (Ep2) from 12 January 2018 to 14 January 2018, with average PM2.5 concentrations of 109.1 μg m −3 and the highest value was 155.5 μg m −3 . Compared with Ep1, the mass concentration of PM2.5 in Ep2 decreased by 49.8%, which implies that the air quality of the Tianjin urban area has been significantly improved after the CRP. The average δ 34 Ssulfate value of Ep2 is 5.1‰, which was 15.7% lower than Ep1. Moreover, it was noteworthy that the δ 34 Ssulfate value of Ep2 was lower than that of the clean period during sampling, which was significantly different from the sampling period before the implementation of CRP (Figure 3b). In addition, the measured δ 34 Ssulfate of Ep2 was positively correlated with T (R = 0.70, p < 0.01). Since the BLH of Ep2 was significantly higher than that of Ep1 and Ep3 ( Figure S1), the air mass had better diffusion ability, so the source and oxidation pathway of sulfate might change. Combined with the backward trajectory analysis (Figure 5), we found that the air mass during Ep2 mainly came from cluster 2 and cluster 4, which were mainly affected by the local sources and regional transportation in the southwest.

Contributions of Oxidation Pathways to Sulfate Formation
Assuming the relative contribution of each source to sulfate aerosols was constant under the different air pollution conditions, the fractions of different sulfate formation pathways were calculated by the MIXSIAR model (Figure 6), and the contributions and uncertainties of different oxidation pathways were shown in Table S3. It could be found that the proportion of sulfate oxidation pathway in the three stages was different. As shown in Figure 6, the contribution of the TMI oxidation pathway and NO 2 oxidation pathway in the three sampling stages was higher than that of the OH radical oxidation pathway and H 2 O 2 /O 3 oxidation pathway, indicating that the formation pathway of sulfate in the Tianjin urban area was mainly dominated by a heterogeneous oxidation pathway. Except for Ep2, the contribution of TMI to sulfate during sampling was consistent with other cities in China [31,32], and the OH radical pathway was also consistent with values from the previous study (19-34%) in Beijing during the polluted season (Table S3) [32]. Thus, the source and formation pathways of sulfate have changed significantly during Ep2 when compared with the other two episodes.
that the proportion of sulfate oxidation pathway in the three stages was different. As shown in Figure 6, the contribution of the TMI oxidation pathway and NO2 oxidation pathway in the three sampling stages was higher than that of the OH radical oxidation pathway and H2O2/O3 oxidation pathway, indicating that the formation pathway of sulfate in the Tianjin urban area was mainly dominated by a heterogeneous oxidation pathway. Except for Ep2, the contribution of TMI to sulfate during sampling was consistent with other cities in China [31,32], and the OH radical pathway was also consistent with values from the previous study (19-34%) in Beijing during the polluted season (Table S3) [32]. Thus, the source and formation pathways of sulfate have changed significantly during Ep2 when compared with the other two episodes. Due to similar regional transport and local source contributions, we assumed that sulfate sources are similar in the winters of the 3 years. The NO2 oxidation pathway was the dominant sulfate oxidation pathway during Ep1 ( Figure 6). Before the implementation of CRP, due to the NO2 pathway being more important in sulfate production under high NO2 and high RH conditions (Table 1) [33]. In addition, the high PM2.5 concentration during Ep1 (Table 1) might have amplified the formation of sulfate aerosol via the NO2 oxidation pathway through heterogeneous reaction, which was consistent with the research results in Beijing [17,27]. The TMIs oxidation pathway dominated the formation of sulfates during Ep2, and the contribution of Ep2 was about two times higher than that of Ep1. After the implementation of CRP, the increased contribution of the TMI pathway likely resulted from the higher efficiency of TMIs (e.g., Fe and Mn) in the atmosphere, which were emitted from local industrial sources in the haze events that enhanced the TMI oxidation rate [29].
Compared with Ep1, the contribution of sulfate formation by the OH oxidation pathway in Ep2 was significantly reduced after the CRP ( Figure 6). This low contribution (5%) of gas-phase sulfate production was probably attributed to weak photochemistry in winter, high PM2.5 concentration, and extremely high heterogeneous and aqueous oxidations in haze episodes [21,34]. The present study highlighted that the oxidations of SO2 with Due to similar regional transport and local source contributions, we assumed that sulfate sources are similar in the winters of the 3 years. The NO 2 oxidation pathway was the dominant sulfate oxidation pathway during Ep1 ( Figure 6). Before the implementation of CRP, due to the NO 2 pathway being more important in sulfate production under high NO 2 and high RH conditions (Table 1) [33]. In addition, the high PM 2.5 concentration during Ep1 (Table 1) might have amplified the formation of sulfate aerosol via the NO 2 oxidation pathway through heterogeneous reaction, which was consistent with the research results in Beijing [17,27]. The TMIs oxidation pathway dominated the formation of sulfates during Ep2, and the contribution of Ep2 was about two times higher than that of Ep1. After the implementation of CRP, the increased contribution of the TMI pathway likely resulted from the higher efficiency of TMIs (e.g., Fe and Mn) in the atmosphere, which were emitted from local industrial sources in the haze events that enhanced the TMI oxidation rate [29].
Compared with Ep1, the contribution of sulfate formation by the OH oxidation pathway in Ep2 was significantly reduced after the CRP ( Figure 6). This low contribution (5%) of gas-phase sulfate production was probably attributed to weak photochemistry in winter, high PM 2.5 concentration, and extremely high heterogeneous and aqueous oxidations in haze episodes [21,34]. The present study highlighted that the oxidations of SO 2 with NO 2 and TMIs were the important formation pathways of sulfate aerosols in Tianjin, especially during the haze event.

Conclusions
After the implementation of CRP, the average concentration of PM 2.5 in January 2018 and January 2019 dropped by approximately 60%. The trend of WSICs was consistent with PM 2.5 . The δ 34 S sulfate in the clean period was lower than that in the pollution period before the CRP, which was mainly affected by sea salt sources with high δ 34 S. However, the δ 34 S sulfate in the clean period was higher than that in the pollution period after the CRP, which was mainly due to the increase of heterogeneous oxidation contribution during the pollution period. The contribution of the TMI oxidation pathway and NO 2 oxidation pathway in the three sampling stages is higher than that of the OH radical oxidation pathway and H 2 O 2 /O 3 oxidation pathway, which shows that the formation pathway of sulfate in the Tianjin urban area is mainly dominated by heterogeneous oxidation pathway. Before the implementation of CRP, the NO 2 oxidation pathway was the dominant sulfate oxidation pathway during Ep1, and the TMI oxidation pathway dominated the formation of sulfates after the implementation of CRP. Therefore, it should be noted that multiple energy-saving and emission reduction policies should be employed to further optimize the energy structure and therefore maximize the environmental benefits of energy conservation and emission reduction policies.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ijerph191710966/s1, Figure S1: Time series of BLH during January 2017, 2018 and 2019; Figure S2: Consumption of electricity and industrial energy above scale, including coal, natural gas, crude oil, and others (coke, oil product, and thermal) in Tianjin from January 2016 to January 2019. Data from Tianjin Municipal People's Government. Values for January and February each year are average consumption estimated by the sum of them. The SO 2 and NO 2 data during the sampling period come from the Air Quality Monitoring Station; Figure S3: Proportion of water-soluble ions in PM2.5 during sampling period; Table S1: Features of sampling phases; Table S2: The concentration of PM2.5, SO 2 , SO 4 2− (µg/m 3 ), and δ 34 S sulfate values of air mass from different clusters; Table S3: The contribution and uncertainties of different oxidation pathway. Institutional Review Board Statement: Not applicable.