Enhanced Receptor Modeling Using Expanded Equations with Parametric Variables for Secondary Components of PM2.5

Receptor modeling provides valuable information to help develop effective control strategies. Additionally, incorporating parametric variables into expanded receptor modeling improves the understanding of formation mechanisms and potential sources of secondary aerosol. This study was conducted in a rural township in central Taiwan, where the air pollution level was comparable with that in the urban area. Bihourly measurements were applied into an enhanced receptor modeling approach using positive matrix factorization (PMF). Eight potential sources, including oil combustion, coal combustion, secondary aerosol related, nitrate‐rich secondary aerosol, biomass burning, industry/vehicle, road dust, and SOM‐rich (dominated by secondary organic matter) secondary aerosol, were identified. SOM‐rich secondary aerosol (24%) contributed the most to PM2.5 mass, followed by biomass burning (19%) and nitrate‐rich secondary aerosol (18%). Contributions from three factors involving secondary formation features accounted for 55% of PM2.5 mass. Through the enhanced modeling approach, photo-oxidation formation, condensation and aqueous phase oxidation of volatile organic compounds, and transport of secondary nitrates from upwind urban area could be potential formation process and sources of secondary aerosol.


INTRODUCTION
Over the past decades, PM2.5 (particulate matter with aerodynamic diameter less than or equal to 2.5 µm) has raised great concern due to its adverse health effects and inevitable exposures (Bell et al., 2013;Kim et al., 2015;Lu et al., 2015). Identifying the PM2.5 sources and quantifying their contributions will be beneficial in improving emission control strategies, therefore reducing PM2.5 exposures. Quantification of source contribution estimates can be achieved through receptor-oriented source apportionment approaches, including effective variance-chemical mass balance (EV-CMB) and positive matrix factorization (PMF) (Hopke, 2016).
The multivariate PMF solution requires less a priori information regarding the emission profiles as compared to EV-CMB. Given a selected number of factors [p in Eq. (1)], the PMF solution iteratively calculates gik (factor contribution) and fkj (factor profile) with non-negativity constraint on them to fit xij (speciated concentration data). The iteration achieves convergence after obtaining small gradient decrease of objective function Q (Paatero and Tapper, 1994;Viana et al., 2008;Belis et al., 2013;Hopke, 2016): where D represents the WD matrix and V is the WS matrix. δi and νi are index parameters obtained from the classification of the wind data. Previous studies demonstrated the benefits of utilizing wind and time variables in the expanded PMF modeling, e.g., a diesel factor can be separated from vehicle exhaust due to reduced rotational ambiguity (Kim et al., 2003b;Begum et al., 2005). Rotational ambiguity, which results from rotation of the factor matrices and could lead to infinite solutions in PMF modeling, can be expressed as Norris et al., 2009;Paatero and Hopke, 2009): Conceptually, additional environmental characteristics, apart from the wind data, also can be considered in the expanded PMF modeling. For example, photochemical strength and relative humidity (RH) could play important roles in the formation process of secondary aerosol and serve as the parametric variables.
To address above issues and to better understand the sources of secondary aerosol, in this Technical Note we attempted to evaluate the effects of applying an expanded PMF model with bihourly speciation data. Parametric variables, including WS, photochemical strength and RH, were imposed to selected factors involving secondary formation features.

Field Data
An intensive field campaign was conducted in Zhushan (23.76°N, 120.68°E), a hillside township located in central Taiwan, in late autumn (27 October-6 November) of 2018. Zhushan is in a rural area with few industrial sources and low road traffic intensity. Nonetheless, PM2.5 mass concentrations in Zhushan were comparable with those in the urban area (Tsai et al., 2016). This phenomenon is partially explained by the fact that particles can be transported from upwind urban area by sea breezes during daytime (Tsai et al., 2008;Tsai et al., 2016). The sampling site was set up in the downtown area of Zhushan township. The northwest and north-northeast winds were most predominant during the study period. Hourly measurements (mean ± SD) of wind speed, temperature, and relative humidity were 0.98 ± 0.59 m s −1 , 23.0 ± 3.3°C, and 77.8 ± 14.0%, respectively.
Bihourly PM2.5 samples were collected on 150 mm Whatman cellulose filters using high-volume samplers (HVS, DHA-80, DIGITEL) at a flow rate of 500 LPM for the measurement of elements and water-soluble ions. Elemental compositions were analyzed using an inductively coupled plasma mass spectrophotometer (ICP-MS, NexION 300X, Perkin Elmer). Water-soluble ions were quantified by ion chromatography (IC, Dionex ICS 1000, Thermo Scientific). Cellulose and Teflon filters are preferred for analysis of inorganic species, especially trace elements, because of their low backgrounds (Upadhyay et al., 2009). The advantage of using cellulose filters is that they can be completely acid-digested and have high recoveries of trace metals in the ICP-MS analysis (Pekney and Davidson, 2005;Upadhyay et al., 2009).
Additionally, three continuous/semi-continuous instruments were incorporated at the same site. Organic matter (OM) concentrations were obtained from Aerosol Mass Spectrometer (AMS, mini-AMS, Aerodyne) datasets. The AMS datasets were not directly used due to the difficulty of harmonizing mass spectra with inorganic speciation data, but deserve further investigation. Elemental carbon (EC) was monitored using a Multi-Angle Absorption Photometer (MAAP, Model 5012, Thermo-Fisher Scientific). PM2.5 mass concentrations were measured using a Tapered Element Oscillating Microbalance (TEOM, 1405DF, Thermo Electron) monitor. Above measurements were further averaged to bihourly means to match the temporal resolutions of filter samples.
Typically, secondary inorganic aerosol (SIA) can be estimated as the summation of ammonium, nitrate and sulfate ions, which represents the reactions of ammonia (NH3) with nitric acid (HNO3 from NO2) and sulfuric acid (H2SO4 from SO2) (Behera and Sharma, 2011;Huang et al., 2016). Secondary organic aerosol (SOA) is thought to be formed through photochemical reactions of volatile organic compounds (VOCs), followed by nucleation and/or condensation, and cannot yet be directly measured (Ervens et al., 2011;Al-Naiema et al., 2018;Belis et al., 2019). In order to determine secondary organic contributions to PM2.5, primary and secondary OM (POM and SOM) were calculated using Eq. (6) (Turpin and Huntzicker, 1995;Chou et al., 2010;Wu and Yu, 2016): where (OM/EC)pri was determined using the minimum R squared (MRS) method (Wu and Yu, 2016). In brief, the MRS method approximate (OM/EC)pri that generates the minimum R 2 between SOM and EC. As shown in Fig. S1, the (OM/EC)pri was estimated to be 3.20 during the study period.

Receptor Modeling
The script-based Multilinear Engine (ME-2) program that allowed modifying and implementing additional equations was used to execute the expanded PMF modeling. In this study, Eq. (4) is modified to incorporate hourly parametric variables, including WS, RH and Ox, as follows: where V, R and O represent the WS, RH and Ox matrices. The classifications of the parametric data were defined to give approximately equivalent sample sizes across the intervals of each index parameter (ν, ρ and ο). During the study period, more than half of the WS measurements were smaller than 1 m s -1 , which could result in unreliable WD data (Kim et al., 2003b;Buset et al., 2006). Therefore, WD was not incorporated in this study. The expanded PMF modeling was conducted through a 2-step run. The base run retrieved source contribution estimates (gik) and source factors (fkj) using the ordinary bilinear equations as Eq.
(2). The expanded run applied additional equations as Eq. (7) to selected factors involving secondary formation features. The objective function Q was re-defined accordingly through adding an expanded term: In addition to the observed data, measurement uncertainties were required as the model input. The uncertainty for each observation was calculated as (Liao et al., 2017b): where MDLj denotes the method detection limit (MDL) of species j. Each concentration value below the MDL was replaced by MDL/2 with its uncertainty set to (5/6) × MDL. Additionally, only species having more than 30% of detectable values were included in the PMF modeling. Although the matrices, V, R and O, in Eq. (7) contained factor elements to be estimated, the chemical mass balance equation [Eq.
(2)] should be more dominant in the fitting process. Therefore, the expanded equations were assigned greater uncertainties [u ' ij in Eq. (8)] than the corresponding measurement uncertainties . The maximum individual column mean (IM) and the maximum individual column standard deviation (IS) were calculated from the residual matrix to approximate the optimal number of factors in the base model run (Lee et al., 1999): An obvious decrease in IM and IS was expected around the optimal number of factors (Lee et al., 1999;Liao et al., 2020). The retrieved profiles were interpreted by examining the mass fractions and explained variations (EVs) of specific species (Lee et al., 1999;Liao et al., 2017b):

RESULTS AND DISCUSSION
The final model input is a matrix with 116 rows (bihourly concentration data) and 28 columns (PM2.5 and its species) after screening invalid data. Table 1 shows the distribution of concentrations of each species included in the PMF modeling. Ammonium ion, nitrate ion, sulfate ion, POM and SOM were the major constituents (> 12% for each) of PM2.5. These components, except for POM, are regarded as secondary aerosol and accounted for over 60% of PM2.5 mass. After a sequence of model runs with different number of factors, an 8-factor solution ( Fig. 1) seemed most reasonable, having a high correlation (0.87) between observed/predicted PM2.5 mass and well model fitting (> 0.68) for all species (Kuo et al., 2014;Liao et al., 2017a).
Factor01, enriched in SOM, sulfate ion, ammonium ion and POM, is characterized by high EV of V that is associated with residual oil combustion (Swietlicki and Krejci, 1996;Lee et al., 1999;Pakkanen et al., 2003;Sun et al., 2004;Pandolfi et al., 2011). In Taiwan, oil-fired boilers, ship emissions and port activities are considered main sources of oil combustion (Hsu et al., 2005;Hsu et al., 2017). Factor02, characterized by high EV of As, Se and Tl along with high loadings of sulfate ion and ammonium ion, is associated with coal combustion (Han et al., 2006;Antonia López Antón et al., 2013). Factor03 is associated with a secondary aerosol related source by high loadings of sulfate ion, nitrate ion, ammonium ion and SOM without specific markers. Factor04, characterized by both high loading and high EV of nitrate ion, is interpreted as nitrate-rich secondary aerosol. Factor05 is associated with biomass burning, characterized by potassium ion, POM and EC (Streets et al., 2001;Cheng et al., 2009;Gugamsetty et al., 2012). Factor06, loaded on nitrate ion, ammonium ion, POM, SOM and EC along with high EVs of Mn, Fe, Cu, Zn, Mo and chloride ion, is interpreted as a mixed source of traffic and industrial emissions (Larsen and Baker, 2003;Han et al., 2006;Gugamsetty et al., 2012;Jain et al., 2017). Factor07 has high EVs of crustal (Al and Ti) and rare earth (Y, La and Ce) elements, and is thus associated with road dust (Kim et al., 2003a;Wimolwattanapun et al., 2011;Gugamsetty et al., 2012;Kara et al., 2015). Factor08 is interpreted as SOM-rich secondary aerosol, characterized by both high loading and high EV of SOM.
The constituents of secondary aerosol (i.e., sulfate ion, nitrate ion, ammonium ion and SOM) were partially attributed to primary sources, e.g., the profile of biomass burning contained a small fraction of these components. Nevertheless, a large proportion of them was differentiated into three factors (Factor03, Factor04 and Factor08) that showed non-specific source characteristics. To better understand the sources of secondary aerosol, the above-mentioned three factors were chosen as the target factors in expanded PMF modeling.
As shown in Fig. 2, both Factor03 and Factor04 showed higher contribution estimates during daytime but with different diurnal patterns. Factor03 exhibited a peak maximum in the early afternoon (12-14 local time), whereas Factor04 had highest contribution in the late morning (10-12 local time). On the other hand, Factor08 showed higher contribution estimates during nighttime. Difference in diurnal patterns indicated different formation mechanisms or source origins among these three factors. Expanded model results for these secondary factors are shown in Fig. 3. Factor03 was not correlated with wind speed but showed obvious trends with increasing photochemical strength and reducing humidity. It might indicate that Factor03 was involved with locally accumulated secondary aerosol formation from photochemical reaction. In contrast, Factor08 was not well correlated with Ox or RH but demonstrated an apparent drop as WS increased, indicating that Factor08 could be associated with condensation of VOCs under low wind speed conditions during nighttime. Saffari et al. (2016) revealed that the nighttime SOC/OC ratio was increased at RH above 70%, indicating enhanced aqueous phase formation of SOC. It may partially explain the decrease of contribution of Factor08 at RH < 69%. Factor04 exhibited slightly increasing trends with increasing WS, escalating Ox and reducing RH, suggesting that Factor04 might involve with photochemical reaction and transportation of secondary aerosol originated from more distant areas.
To examine the transport of air pollutants from upwind source regions, back trajectory analyses were applied to the top ten contribution periods of Factor04 using the HYbrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT-4) model (Draxler and Rolph, 2015;Rolph, 2015;Stein et al., 2015). Twenty-four-hour back trajectories were computed using meteorological data from the Weather Research and Forecasting (WRF) model simulation with a 1-km grid spacing (Skamarock and Klemp, 2008;Lin et al., 2017). The initial meteorological condition was from ds083.3 NCEP GDAS/FNL 0.25 Degree Global Tropospheric Analyses and Forecast Grids (DOI: 10.5065/D65Q4T4Z, https://rda.ucar.edu/datasets/ds083.3/). In order to get a better meteorological field, the WRF modeling applied four-dimensional data assimilation with observation nudging. The back trajectory analyses results were demonstrated on a map created using the free and open source QGIS (Fig. S2). NOx emission inventory data obtained from Taiwan Emission Data System version 10.0   (TEDS 10.0) was also illustrated on the map. Results indicated that high contribution periods of Factor04 had the pathway passing through urban areas with high NOx emissions, specifically the Taichung metropolis. Since low NOx emissions were reported at the sampling site, upwind urban area was likely the origin of high contributions of secondary nitrates. Average source contribution estimates to PM2.5 mass during the campaign period are shown in Fig. 4. The greatest contributor was SOM-rich secondary aerosol (Factor08, 24%), followed by biomass burning (Factor05, 19%) and nitrate-rich secondary aerosol (Factor04, 18%). Contributions from three factors involving secondary formation features (Factor03, Factor04 and Factor08) accounted for 55% of PM2.5 mass, which could have a considerable impact on air quality at the Zhushan area. In addition to photo-oxidation formation, condensation and aqueous phase oxidation of VOCs and transport of secondary nitrates could be important PM2.5 sources during the campaign period. On the other hand, previous studies also indicated biomass burning as an important source in central Taiwan partially due to the intense burning activities in the harvest season (Yang et al., 2006;Cheng et al., 2009;Liao et al., 2015;Tsai et al., 2016;Lee et al., 2019).
It should be noted that although several previous studies have revealed better results in the expanded PMF modeling (Kim et al., 2003b;Begum et al., 2005;Buset et al., 2006), no improvement was found than the normal PMF analysis in other researches, e.g., Zhou et al. (2009). Thus, multiple models should be explored to ensure that the most reliable solution has been obtained. Another limitation of this study is that the intensive 10-day field campaign was unable to assess long-term exposure to PM2.5 in the study area. In addition, time variables such as season and weekend effects were not considered in this study due to the same issue. Nonetheless, this study exhibited the benefits of application of expanded PMF modeling with parametric variables in improving characterization of the formation mechanisms and potential sources of secondary aerosol. Further investigation of SOA tracers would be beneficial in identifying the source origins of VOCs that are associated with SOA formation.

CONCLUSIONS
In this study an expanded PMF modeling approach was applied to bihourly PM2.5 speciation data. Results indicated that SOM-rich secondary aerosol contributed the most (24%) among the eight apportioned factors. Through the enhanced modeling approach, this factor could be associated with condensation and aqueous phase oxidation of VOCs, suggesting the need for further investigation of SOA tracers. The second largest contributor (biomass burning, 19%) might be abated through regulating open field burning activities. The third largest factor (nitrate-rich secondary aerosol, 18%) was characterized to be more likely the transport of secondary nitrates from upwind metropolitan area during the campaign period. Cooperation between neighboring cities would be more cost-effective in reducing its contribution to PM air pollution.

DISCLAIMER
The authors declare that there are no conflicts of interest.

SUPPLEMENTARY MATERIAL
Supplementary data associated with this article can be found in the online version at https://doi.org/10. 4209/aaqr.200549