Development of a CMAQ–PMF-based composite index for prescribing an effective ozone abatement strategy: a case study of sensitivity of surface ozone to precursor volatile organic compound species in southern Taiwan

Photochemical ozone pollution is a serious air quality problem under weak synoptic conditions in many areas worldwide. Volatile organic compounds (VOCs) are largely responsible for ozone production in urban areas where nitrogen oxide (NO_x) mixing ratios are high while usually not a limiting precursor to ozone (O₃). In this study, the Community Multiscale Air Quality model higher-order decoupled direct method (CMAQ-HDDM) at an urban-scale resolution (1.0 km×1.0 km) in conjunction with positive matrix factorization (PMF) was used to identify the dominant sources of highly sensitive VOC species to ozone formation in southern Taiwan, a complex region of coastal urban and industrial parks and inland mountainous areas. First-order, second-order, and cross sensitivities of ozone concentrations to domain-wide (i.e., urban, suburban, and rural) NO_x and VOC emissions were determined for the study area. Negative (positive) first-order sensitivities to NO_x emissions are dominant over urban (inland) areas, confirming ozone production sensitivity favors the VOC-limited regime (NO_x-limited regime) in southern Taiwan. Furthermore, most of the urban areas also exhibited negative second-order sensitivity to NO_x emissions, indicating a negative O₃ convex response where the linear increase of O₃ from decreasing NO_x emissions was largely attenuated by the nonlinear effects. Due to the solidly VOC-limited regime and the relative insensitivity of O₃ production to increases or decreases of NO_x emissions, this study pursued the VOC species that contributed the most to ozone formation. PMF analysis driven by VOCs resolved eight factors including mixed industry (21 %), vehicle emissions (22 %), solvent usage (17 %), biogenic sources (12 %), plastic industry (10 %), aged air mass (7 %), motorcycle exhausts (7 %), and manufacturing industry (5 %). Furthermore, a composite index that quantitatively combined the CMAQ-HDDM sensitivity coefficient and PMF-resolved factor contribution was developed to identify the key VOC species that should be targeted for effective ozone abatement. Our results indicate that VOC control measures should target (1) solvent usage for painting, coating and the printing industry, which emits abundant toluene and xylene; (2) gasoline fuel vehicle emissions of n-butane, isopentane, isobutane, and n-pentane; and (3) ethylene and propylene emissions from the petrochemical industry.

1 Introduction

Photochemical production of tropospheric ozone (O₃) depends in a nonlinear manner on the availability of nitrogen oxides (NO_x) and volatile organic compounds (VOCs). Understanding ozone formation sensitivity to NO_x and VOC emissions is key to developing effective abatement strategies on O₃ pollution in heavily polluted cities. The brute-force method (BFM) has often been used to address the relationship between O₃ and its precursors (Hakami et al., 2004b; Li et al., 2013a; Zhang et al., 2009). In a BFM approach within a 3D-modeling framework (e.g., Community Multiscale Air Quality modeling system – CMAQ), individual emissions are perturbed by a small amount, and the model response is recorded against the baseline, representing the sensitivity coefficient. However, the linear response of BFM is insufficient to account for changes in secondary pollutants (i.e., ozone) generated by nonlinear interactions of various substances. Therefore, more sophisticated numerical techniques have been introduced such as the adjoint method (Hakami et al., 2006; Wang et al., 2021); the decoupled direct method in three dimensions, DDM 3D (Dunker et al., 2002b; Luecken et al., 2018); and the higher-order decoupled direct method, HDDM 3D (Cohan et al., 2005; Koplitz et al., 2021), where the 3D aspect is specific to implementation in chemical transport models (CTMs) with a 3D-modeling framework (e.g., CMAQ). These methods offer an alternative to BFM by directly solving the auxiliary equations that represent the change in concentration over change in emission ($\mathrm{d}C/\mathrm{d}E$) derived from the governing equations of the model.

DDM 3D, using a semi-implicit finite-difference scheme with different time steps (Hakami et al., 2004a; Yang et al., 1997), has been widely implemented in CTMs (Dunker et al., 2002b, a) to calculate the first-order sensitivities (i.e., tangential, $\partial C/\partial E$) of ozone with respect to initial concentration, boundary concentrations, and precursor emission rates. First-order sensitivity describes the linear response of the model to a perturbed input parameter, while higher-order sensitivity describes the nonlinear quadratic, cubic, and higher-power responses. In the presence of strong nonlinearity (e.g., in the transition between the VOC-limited and NO_x-limited regime), first-order sensitivity alone may be insufficient to characterize the response if the magnitude of the emission changes are large (i.e., a manifestation of the nonlinear chemistry involved in ozone formation), and second-order sensitivity is necessary to provide the additional nonlinear responses of the system (Xing et al., 2011; Xiao et al., 2010). Thus, in HDDM, the sensitivity is calculated by solving the first-order ($\partial C/\partial E$) and second-order (i.e., local slope, ${\partial }^{\mathrm{2}}C/\partial E^{\mathrm{2}}$) derivatives in the model when a relative perturbation (e.g., ±10 %) is applied to a targeted parameter while keeping all other factors constant.

Apportioning pollutant concentrations to their sources is an ideal strategy to guide emission control policy (Dunker et al., 2015; Liu et al., 2008). However, source apportionment of a secondary pollutant such as ozone is complex, where a relatively crude method is to simply conduct simulations with and without NO_x and VOC emissions from a given source so that the difference in ozone concentration is a measure of the source contribution (Bergin et al., 2008). Other source apportionment approaches such as the ozone source apportionment technology (OSAT) and integrated source apportionment method (ISAM) are embedded into air quality models (Comprehensive Air Quality Model with Extensions (CAMx) and CMAQ, respectively) and rely on tagging NO_x and VOCs as tracers from emission to ozone production to estimate the contributions of different sources to the eventual ozone concentration (Dunker et al., 2002a; Kwok et al., 2015). Unless a high-resolution sector profile emission inventory is available, the application of OSAT and ISAM is often limited to four main sector groups: on-road, non-point (area), point, and biogenic sources (Wang et al., 2009; Li et al., 2013b). In addition, the spatiotemporal and sector group distribution uncertainties associated with the emission inventory greatly influence the accuracy of these source apportionment methods (Zheng et al., 2009).

Positive matrix factorization (PMF), a receptor-based approach, offers an alternative to those embedded methods by generating a set of ozone precursor source composition profiles, each identifying a mix of compounds associated with a particular category of emissions. Driven by measured concentrations at receptor sites, PMF can be used as an ozone source apportionment method directly and is independent of the uncertainties associated with the emission inventory. In recent years, PMF has been widely used to estimate VOC apportionment because it requires only identifying characteristics of the source profiles to interpret the PMF factors (Ji et al., 2022; Fan et al., 2021; Huang and Hsieh, 2020; Chen et al., 2019; Wu et al., 2016). For instance, Pallavi et al. (2019) found that traffic contributed 47 % of the total benzene in India followed by residential biofuel use and waste disposal (25 %) and industrial emissions and solvent use (20 %). Zhao et al. (2020) used PMF in conjunction with an ozone formation potential (OFP) calculation to infer that VOCs from industrial and vehicular emissions were the dominant ozone precursors in Nanjing, East China, particularly xylenes, toluene, and propene. Chen et al. (2019) used PMF to resolve the dominant VOC sources at an industrial complex in central Taiwan and found that the monitored VOC concentrations of vehicle exhaust, solvent use, and diesel attributed to high OFP were associated with easterly to southeasterly winds. Huang and Hsieh (2020) analyzed VOC data by PMF from western coastal Taiwan and suggested that on a mass concentration basis industrial emissions are the greatest contributors to OFP. However, previous PMF–OFP studies have only identified key VOC sources to OFP without investigating the sensitivity of OFP to these VOC emissions and thus failed to provide comprehensive advice on which sources to prioritize for effective ozone abatement.

In this work, we used CMAQ-HDDM-3D in conjunction with PMF, enabling us to interpret the sensitivity and source apportionment analyses together for a more comprehensive investigation of the key VOC sources to OFP. We combined the two approaches as a novel methodology to provide additional insights on the dominant sources of highly sensitive VOC species to ozone formation in a VOC-limited urban area of Taiwan, but this should be widely applicable across urban areas that experience similar O₃ episodes such as Hong Kong (Ling and Guo, 2014); Beijing and Hebei (Chi et al., 2018); Mexico City (Lei et al., 2007); and Houston, United States (Mazzuca et al., 2016). Specifically, HDDM describes the O₃ sensitivity response to a speciated emission and PMF identifies the sources of these species for effective ozone abatement strategy. While most of the HDDM sensitivity studies are performed at coarse resolution (>4.0 km), this study was conducted at an urban-scale resolution (1.0 km×1.0 km). We investigated the first-order, second-order, and cross sensitivities of ozone concentrations to domain-wide NO_x and VOC emissions and provide an overview of the O₃-precursor sensitivity across the study area. We then identified the VOC species to which ozone formation was most sensitive. Finally, we mapped these highly sensitive VOC species to the PMF model source apportionment and identified their dominant sources for effective emission control strategies using a composite index of CMAQ sensitivities and PMF factor contributions. The results of our study provide important information as to which VOC species are key to ozone formation and where to reduce sources of these VOC species for effective ozone abatement.

2 Methods

2.1 Study period and area

In this study, we selected the period from 7–20 October 2018 for conducting simulations of photochemical ozone production and transport in southern Taiwan. The selected case in October 2018 is the seasonal transition period when the summer season is in transition to the winter. The case can reasonably represent the typical ozone pollution conditions during the seasonal transition period in Taiwan because the synoptic weather pattern of the event features a weak intrusion of the Asian continental anticyclone system which slowly propagated eastward, causing the prevailing wind in Taiwan dominated by weak northeasterly (NE) flows due to continental high-pressure peripheral circulation (see Fig. S1 in the Supplement). Hsu and Cheng (2019) identified six synoptic weather patterns common in Taiwan and studied the characteristics of corresponding air pollutants in each pattern using 6 years (2013–2018) of daily averaged wind fields and sea-level pressure observed at surface weather stations in Taiwan. Among the six patterns (C1–C6), C4 has the highest mean O₃ concentrations and occurs predominantly in October. It features a weak anticyclone over the Asian continent, and the Pacific subtropical high-pressure system does not have an apparent influence in Taiwan, which is similar to our selected case in October 2018. Although the photochemistry is strong in the summer season, the seasonal O₃ variation in Taiwan shows that the monthly O₃ concentration is relatively higher during the seasonal transition months (i.e., October) compared to other seasons (Hsu and Cheng, 2019; Chen et al., 2021; Cheng et al., 2015). This is because during the seasonal transition months, when the photochemical reaction is still strong compared to that of the winter months together with the reduced ventilation capability, the O₃ concentration can accumulate (Yen and Chen, 2000; Tsai et al., 2008) (see Fig. S2 in the Supplement for ozone seasonality in Taiwan). During the seasonal transition period in autumn, southern Taiwan often suffers from high-O₃ episodes (Hsu and Cheng, 2019; Chen et al., 2004, 2021). Regional synoptic weather in autumn usually features a weak anticyclone over the Asian continent, allowing for local accumulation of pollutants, while the Pacific subtropical high-pressure system has shifted eastward with no apparent influence on Taiwan (Hsu and Cheng, 2019). This synoptic weather pattern is also categorized as high-pressure pushing (HPP), which occurred when the leading edge of Asian continental high-pressure systems moves over China coastal provinces and carries pollutants southward towards downwind area in Taiwan (Chuang et al., 2008). In addition, precipitation is lower and vertical dispersion is weaker than in summer (Yen and Chen, 2000; Tsai et al., 2008), which contribute to relatively more frequent ozone episodes in autumn. These conditions are also common occurrences for high-O₃ episodes in other places such as Hong Kong and East China cities in autumn (Lee et al., 2009; Liu et al., 2021; R. Yu et al., 2021). A 5 d spin-up period (2–6 October 2018) was discarded from the analysis to eliminate the effects of initial conditions on O₃ simulations.

Since the 1990s, ozone concentrations have been increasing in many areas of southern Taiwan (Chang et al., 2005; Chen et al., 2021). Kaohsiung city, the biggest city in southern Taiwan and second biggest in all of Taiwan, hosts many heavy industries, including petrochemical plants, refineries, steel-making plants, and power generation plants. Xiaogang (XG) district is located in the southeastern portion of the city and hosts a petrochemical industrial park with an overall area of 403 ha. Three power plants, each producing 1200–4300  MW d⁻¹, are located within a distance of 35 km of XG and lead to some of the poorest air quality in Taiwan.

2.2 WRF–CMAQ model configuration

The Weather Research and Forecasting (WRF) model version 3.9.1 and CMAQ model version 5.2.1 (WRF–CMAQ) is used in the current study with a domain configuration of a four-nested grid system centered at 28^∘ N and 125^∘ W (Fig. 1a). The outermost domain (D01) covers most of mainland China with a horizontal resolution 27 km×27 km and 166×169 grids, which is then nested to the second domain (D02) of 9 km×9 km and 223×223 grids over East China. The third domain (D03) covers the whole island of Taiwan with a resolution of 3 km×3 km and 223×223 grids, and the innermost domain (D04) focuses on southern Taiwan with a 1 km×1 km resolution and 136×136 grids. There are 41 vertical sigma layers spaced unequally from the ground to the model top 50 hPa, with the bottom 8 layers resolved below 1.5 km (Fig. S3 in the Supplement).

Figure 1 (a) Domain configuration of a four-nested grid system; (b) land use of the innermost domain; urban and inland areas are represented by class 31, 32, 33, and 6, respectively, (c, d) monthly averaged NO_x and VOC emissions in the innermost domain obtained from the 2016 TEDS-10 emission inventory. The locations of each TEPA air quality station (red stars) and PAMS station (red dots with label) used in the current study are displayed in the innermost domain. Refer to Fig. S3 and Table S3 in the Supplement for details of each TEPA and PAMS station.

The initial and boundary meteorological conditions are adopted from the National Centers for Environmental Prediction (NCEP) global final analysis (FNL) data at $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ resolution updated every 6 h. The projected 2017-year Multi-resolution Emission Inventory for China (available at http://meicmodel.org, last access: 30 January 2022) with $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ resolution is used for D01 and D02, while the 2016 Taiwan Emission Data System (TEDS) version 10 (Fig. 1c and d) is used for D03 and D04. Biogenic emissions are calculated offline using the Model of Emissions of Gases and Aerosols from Nature (MEGAN) version 2.10 (Guenther et al., 2006).

Other than local anthropogenic emissions, the contribution of long-range transport (LRT) from East Asia (e.g., Chinese emissions) is also substantial to Taiwan's air quality, especially under strong northeasterly wind conditions (Wu and Huang, 2021; C.-T. Chang et al., 2022). Lin et al. (2004) identified three types of common LRT events in Taiwan: (1) dust storm (DS), (2) LRT with pollutants (frontal pollution; FP), and (3) LRT of background air masses (BG). When there is no frontal system, local pollution (LP) dominates the air quality in Taiwan. During wintertime and springtime, the occurrence of LP cases was 70 %, and about 30 % were LRT cases (Lin et al., 2004). Lin et al. (2005) estimated that the long-range transport of pollutants contributes to about 30 µg m⁻³, 230 ppb, and 0.5 ppb to the PM₁₀, CO, and SO₂ concentrations, respectively, in northern and eastern Taiwan. Meanwhile a smaller contribution is estimated in southern Taiwan due to the geography (Lin et al., 2005).

In our work, the oceanic chlorine emission is calculated online by CMAQ as a function of meteorology using an OCEAN file which specifies the fraction of each grid cell that is open ocean (OPEN) or surf zone (SURF). Figure 2a–c present the spatial distribution of CMAQ-calculated sea-salt aerosol cations (ASEACAT – Na⁺, K⁺, Ca²⁺, and Mg²⁺), fine-mode chlorine, and coarse-mode chlorine emission rates averaged during the entire simulation period. The sea-spray emissions were higher in the surf zone area, and the highest emission rates were found along the eastern offshore area of southern Taiwan. This is because of the enhanced formation of sea-spray aerosols associated with higher relative humidity and greater offshore winds along the eastern offshore area of southern Taiwan that is open to the western Pacific Ocean. Besides, the anthropogenic chlorine (PCL) emissions are obtained from TEDS v10 emissions, and they are concentrated over the heavily industrialized urban areas of southern Taiwan (see Fig. 2d).

Figure 2 CMAQ-calculated (a) sea-salt aerosol cation (ASEACAT) emissions (Na⁺, K⁺, Ca²⁺, and Mg²⁺), (b) fine-mode chlorine sea-spray aerosol (SSA) emissions, (c) coarse-mode chlorine SSA emissions, and (d) TEDS v10 anthropogenic chlorine (PCL) emissions averaged during the entire simulation period.

To obtain more accurate dynamical downscaling, grid nudging was applied in the coarse domains D01 and D02; observation nudging was applied in the fine domains D03 and D04. Grid nudging is applied to the horizontal wind components, potential temperature, and water vapor mixing ratio; it is only applied above the planetary boundary layer (PBL). The observational data for observation nudging include hourly surface observations such as atmospheric pressure, air temperature, relative humidity, wind speed, and wind direction from 36 surface meteorological stations (https://www.epa.gov.tw/, last access: 30 January 2022), as well as the twice-daily at 00:00 and 12:00 UTC sounding data such as potential height, temperature, dew point temperature, RH, wind direction, and wind speed at each specified isobaric level from two radiosonde observation stations in Taiwan. The nudging coefficients, which determine the strength of the assimilation tendency term, were set to be $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{4}}$ for observation nudging and $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ for grid nudging. These values of coefficients were recommended by the WRF user guide and tested to be appropriate in previous studies (Li et al., 2022; Borge et al., 2008).

The simulation adopted the Carbon Bond Mechanism CB6 (Yarwood et al., 2010), which was developed as an update to CB05 to provide a condensed chemical mechanism for use in atmospheric models. CB6 includes five additional organic compounds that are long-lived and relatively abundant (i.e., propane, acetone, benzene, ethyne, and higher ketones) and a more detailed representation of organic nitrate reactions. CB6 has been tested against measurements over a wide variety of spatial, temporal, chemical, and meteorological conditions and has shown good agreement with measurements of ozone and nitrogen oxides in both urban and rural areas (Luecken et al., 2019). The halogen chemistry in CB6 considers chlorine-related reactions such as ClNO₂, HCl, and HNO₃ production from heterogeneous uptake of N₂O₅ on the aerosol surface, which are important to ozone pollution over coastal cities in southern Taiwan.

The Noah land surface model (LSM) is used to describe the land–atmosphere interaction (Chen and Dudhia, 2001). The urban effect in Noah LSM is invoked by implementing a single-layer urban canopy model (UCM), which assumes infinitely long 2D street canyon urban geometry to improve the processes associated with the exchange of momentum, heat, and moisture in the urban environment (Kusaka and Kimura, 2004; Kusaka et al., 2001). To take full advantage of the UCM scheme, three additional urban classes (i.e., low-density residential, high-density residential, and industrial/commercial) are further classified for better representation of the urban features (Fig. 1b).

The asymmetric convective model version 2 (ACM2) boundary layer scheme is selected to represent the boundary layer process (Pleim, 2007). It includes the nonlocal scheme of the original ACM combined with an eddy diffusion scheme. Thus, ACM2 is able to represent both the supergrid- and subgrid-scale components of turbulent transport in the convective boundary layer. Also used were the Goddard Cumulus Ensemble (GCE) microphysics scheme (Tao et al., 2014), Rapid Radiative Transfer Model (RRTM) longwave radiation scheme (Gallus and Bresch, 2006), Goddard shortwave radiation scheme (Chou and Suarez, 1999), Monin–Obukhov similarity scheme, and Kain–Fritsch cumulus parameterization scheme (only at D01 and D02).

The modeled results of both meteorology (i.e., 2 m temperature, wind speed and direction, and relative humidity) and air quality (i.e., ozone, nitrogen oxides, and volatile organic compounds) are validated against 15 Taiwan Environmental Protection Administration (TEPA) air quality stations. Overall, the modeling system reproduced the observed meteorological and air quality conditions within the benchmark values (see Supplement – “Model Evaluations”).

The urban and inland grid cells are defined according to the USGS-24 land use category. The urban area is represented by class 1 – urban and built-up land, which we further classified into class 31, 32, and 33 (see Fig. 1b) for the WRF single-layer urban canopy model (SLUCM) simulation. We refer the readers to our previous work for detailed discussion on the land use classification and SLUCM implementation in J. H. W. Chang et al. (2022). The inland area is represented by class 6 – cropland/woodland mosaic (see Fig. 1b). The general meteorological pattern of the event features a weak intrusion of the Asian continental anticyclone system which slowly propagated eastward, causing the prevailing wind at a synoptic scale in Taiwan dominated by weak northeasterly (NE) flows due to continental high-pressure peripheral circulation (see Fig. S1). At a local scale in southern Taiwan, the steering of weak NE flows by the orographic effect of the Central Mountain Range (CMR) enhanced the local circulations (i.e., land–sea breeze) and eventually pushed the locally produced urban O₃ as well as its precursors NO_x and non-methane hydrocarbon (NMHC) towards the inland areas (see Fig. 3).

Figure 3 (a) Topography of the innermost domain. (b, c) Spatial distribution of O₃ concentration averaged during the entire simulation at 10:00 and 20:00 LST, respectively. (d–f) Vertical profile of NMHC concentration cross sectioned at AB (see Fig. 3a) averaged during the entire simulation at 10:00, 15:00, and 20:00 LST, respectively. (h–i) Same as (d–f) but for NO₂ concentration.

2.3 Higher-order decoupled direct method (HDDM)

The response of a chemical concentration to perturbations in model parameters (i.e., emissions, initial condition, boundary condition, and reaction rate constants) can be evaluated through sensitivity analysis. A perturbed sensitivity parameter, p_i, is related to the unperturbed sensitivity parameter, P_i, in the baseline simulation as

$$\begin{array}{}\text{(1)}& p_i={\mathit{\epsilon }}_i{P}_i=(\mathrm{1}+\mathrm{\Delta }{\mathit{\epsilon }}_i)P_i,\end{array}$$

where ε_i is a scaling factor with a nominal value of 1.0 when there is no perturbation, and Δε_i is a perturbed scaling factor. The response of a chemical concentration, C, against the perturbations in a sensitivity parameter, p_i, is defined as the sensitivity coefficient, S_i. The first- and second-order sensitivity coefficients, $S_i^{\left(\mathrm{1}\right)}$ and $S_{i,j}^{\left(\mathrm{2}\right)}$, are defined as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }{S}_i^{\left(\mathrm{1}\right)}=P_i{\displaystyle \frac{\partial C}{\partial p_i}}=P_i{\displaystyle \frac{\partial C}{\partial \left({\mathit{\epsilon }}_i{P}_i\right)}}={\displaystyle \frac{\partial C}{\partial {\mathit{\epsilon }}_i}},\text{(3)}& {\displaystyle }{S}_{i,j}^{\left(\mathrm{2}\right)}=P_i{\displaystyle \frac{\partial C}{\partial p_i}}{P}_j{\displaystyle \frac{\partial C}{\partial p_j}}=P_i{\displaystyle \frac{\partial C}{\partial \left({\mathit{\epsilon }}_i{P}_i\right)}}{P}_j{\displaystyle \frac{\partial C}{\partial \left({\mathit{\epsilon }}_j{P}_j\right)}}={\displaystyle \frac{{\partial }^{\mathrm{2}}C}{\partial {\mathit{\epsilon }}_i\partial {\mathit{\epsilon }}_j}}.\end{array}$$

Both S_i(1) and S_i,j(2) have the same units as the chemical concentration, C. S_i(1) measures the impact of one variable p_i on a concentration, C; S_i,j(2) measures the impact of another variable p_j on a first-order sensitivity, S_i(1), which can be used to investigate the second-order cross sensitivity of the system. When i=j, S_i,i(2) represents the local curvature of the relationship between a concentration and a single parameter, thus indicating the responsiveness of C to a broad range of p_i and describing a nonlinearity of the system. A second-order sensitivity of zero indicates a perfectly linear response, while greater values represent proportionally greater nonlinearity effects. Matching signs of the sensitivities (e.g., positive first order and positive second order) represent a convex response while opposite signs represent a concave response (Fig. S4 in the Supplement). In the context of our study, reduced emission is the focus so a convex response represents a nonlinear ozone concentration decrease with decreasing emissions, and a concave response represents a nonlinear ozone concentration increase with decreasing emissions.

Taylor-series expansion

After determining the sensitivity coefficients, we can apply these coefficients for estimating the concentration changes in emission reduction scenarios based on Taylor-series expansion. The concentration changes from any fractional perturbations in a sensitivity parameter are approximated by (Cohan et al., 2005)

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{C}_i& {\mathrm{|}}_{p_i=P_i+\mathrm{\Delta }{\mathit{\epsilon }}_i{P}_i}\\ & \approx C_o{\mathrm{|}}_{p_i=P_i}+\mathrm{\Delta }{\mathit{\epsilon }}_i{S}_i^{\left(\mathrm{1}\right)}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }{\mathit{\epsilon }}_i^{\mathrm{2}}{S}_{i,i}^{\left(\mathrm{2}\right)}+\mathrm{\dots }+{\displaystyle \frac{\mathrm{1}}{n}}\mathrm{\Delta }{\mathit{\epsilon }}_i^n{S}^{\left(n\right)},\end{array}\end{array}$$

where C_i is the concentration when p_i has been perturbed by an amount Δε_iP_i. Subscript i represents the targeted emission species (i.e., a NO_x or VOC species). When multiple sensitivity parameters are perturbed simultaneously, a second-order Taylor approximation includes the interaction between the two parameters (Cohan et al., 2005):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{C}_{i+j}\cong & C{\mathrm{|}}_{p_{i,j}=P_{i,j}+\mathrm{\Delta }{\mathit{\epsilon }}_{i,j}{P}_{i,j}}\\ \approx & C_o+\mathrm{\Delta }{\mathit{\epsilon }}_i{S}_i^{\left(\mathrm{1}\right)}+\mathrm{\Delta }{\mathit{\epsilon }}_j{S}_j^{\left(\mathrm{1}\right)}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }{\mathit{\epsilon }}_i^{\mathrm{2}}{S}_{i,i}^{\left(\mathrm{2}\right)}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }{\mathit{\epsilon }}_j^{\mathrm{2}}{S}_{j,j}^{\left(\mathrm{2}\right)}+\mathrm{\Delta }{\mathit{\epsilon }}_i\mathrm{\Delta }{\mathit{\epsilon }}_j{S}_{i,j}^{\left(\mathrm{2}\right)}+\mathrm{\dots },\end{array}\end{array}$$

where S_i,j is the cross sensitivity between the two parameters, which differs from the sum of the two sensitivities, S_i and S_j.

In this study, CMAQ-HDDM-3D v5.2.1 is used to calculate the sensitivity coefficients to predict the O₃ response to perturbations in NO_x and VOC emissions. Perturbations at ±10 % were made on the domain-wide emissions of NO_x and VOCs (D01–D04) from both anthropogenic and biogenic sources, but only sensitivity coefficients from the innermost domain (D04) are used for data analysis. A total of 25 sensitivity tests were performed to quantify the change in O₃ concentration due to the perturbations made in each test (Table 1); the first 5 sensitivity tests, S1–S5, account for the first-order, second-order, and cross sensitivity due to the perturbed NO_x and total VOC emissions, and the other 20 sensitivity tests, S6–S25, account for sensitivities due to perturbed individual VOC model species. Note that the HDDM approach was only used in experiment S3–S5, which involves calculation of the higher-order sensitivity, while the DDM approach was used in all other experiments and in conjunction with PMF analysis.

Table 1 Perturbed emissions considered in the 25 sensitivity tests. The first 5 sensitivity tests S1–S5 account for the first-order, second-order, and cross-order sensitivity due to the domain-wide NO_x and VOC emissions, and the other 20 sensitivity tests S6–S25 account for the individual VOC model species.

2.4 Positive matrix factorization (PMF) model

The positive matrix factorization (PMF) model is a source-receptor statistical factor analysis method, which is widely used for apportioning sources of air pollution and resolves the dominant factor profile without prior knowledge of sources; the measured data uncertainty is used to optimize the model (Daellenbach et al., 2017; Fountoukis et al., 2014; Fan et al., 2021). In this study, the US EPA PMF 5.0 receptor model was used for the source apportionment of measured VOC species:

$$\begin{array}{}\text{(6)}& x_{ij}=\sum _{k=\mathrm{1}}^p{g}_{ik}{f}_{kj}+e_{ij},\end{array}$$

where x_ij is the jth species concentration measured in the ith sample, g_ik is the airborne mass contribution from the kth source in the ith sample, f_kj is the jth species fraction to the kth source, and e_ij is the residual associated with the jth species concentration measured in the ith sample. p is the total number of independent sources. Note that the subscript “i” in Eq. (6) refers to the ith sample, which is different from Eqs. (1–5) that refer to the emission of the ith species. In the source parsing process, the objective function Q is solved with an iterative minimization algorithm. The objective function is defined as

$$\begin{array}{}\text{(7)}& Q=\sum _{i=\mathrm{1}}^n\sum _{j=\mathrm{1}}^m{\displaystyle \frac{x_{ij}-{\sum }_{k=\mathrm{1}}^p{g}_{ik}{f}_{kj}}{u_{ij}}},\end{array}$$

where u_ij is the uncertainty of the jth species in the ith sample.

Figure 4 shows the overall framework of the PMF source apportionment methodology. The data used to drive the PMF model source apportionment were obtained from Taiwan's Photochemical Assessment Monitoring Stations (PAMS); details of the PAMS location and sampling protocol are provided in the Supplement. Species with >55 % of samples missing data or below the minimum detection limit (MDL) were discarded (Chen et al., 2019; Wu et al., 2016). Among 54 PAMS-VOC species, a total of 16 species were discarded due to abundant missing data but had a low OFP (Fig. S5 in the Supplement). The sample data uncertainty is calculated using

$$\begin{array}{}\text{(8)}& \text{unc}=\left\{\begin{array}{ll}\mathrm{5}/\mathrm{6}\times \text{MDL},& x_{ij}<\text{MDL}\\ \sqrt{(\mathit{\sigma }\times x_{ij}{)}^{\mathrm{2}}+(\mathrm{0.5}\times \text{MDL}{)}^{\mathrm{2}}},& x_{ij}\ge \text{MDL},\end{array}\right.\end{array}$$

where MDL represents the minimum detection limit, σ is the error fraction (i.e., 10 % was used in this study), and x_ij is the concentration of the jth species in the ith sample. In the PMF, we tested a range of three to eight VOC sources to determine the optimal number. Details of this protocol will be discussed further below.

Figure 4 Overall framework of the positive matrix factorization (PMF) model methodology. Processes in the dashed-line box are repeated for three to eight factor number combinations. MDL: minimum detection limit.

The Taiwan PAMS data provided a total of 54 VOC species, but not all species were used for the PMF model due to an abundance of values below MDLs. Among the selected 38 VOC species (n=744 samples), they were further categorized into three categories: strong, weak, and poor according to two steps. The first step calculated the signal-to-noise ratio ($S/N$) where species with a $S/N$ value less than 0.5 were classified as poor and removed from the PMF model. This threshold value 0.5 was determined according to the EPA PMF v5.0 user guide and also recommended by other PMF studies (Rajput et al., 2016; Reff et al., 2007). Next, we performed the model base run and calculated the correlation between the modeled and measured concentration for each species. In this step, species with a correlation value ≥0.6 were classified as “strong” and otherwise as “weak”, which were then down-weighted by tripling the analytical uncertainty (Pallavi et al., 2019; Chen et al., 2019). Details of the species categorization are summarized in Table 2. Finally, a total of 21 VOC species were identified and used in the PMF source apportionment analysis, which accounted for 75.0 %, 76.4 %, and 76.1 % of the total VOC concentrations in CZ, QT, and XG, respectively.

Table 2 Categorization of PAMS-VOC species for PMF model source apportionment analysis. Bold species represent unused species with abundant missing data >55 % below the minimum detection limit (MDL). Poor category species are identified for low $S/N<\mathrm{0.5}$. Weak (Strong) category species are identified for $S/N\ge \mathrm{0.5}$ and R²<0.6 (R²≥0.6) between the measured and modeled concentration predicted by the PMF model. Both unused and bad species are removed from the PMF model analysis.

3 Results and discussions

3.1 Decomposition of ozone response

First-order sensitivity coefficients indicate the linear response of ozone concentrations to small changes in emissions. The ozone response to small changes of daytime (09:00–15:00 LST) NO_x emissions is dominated by negative sensitivities (VOC-limited) in the urban area and positive sensitivities (NO_x-limited) in the inland area near the mountainous region (Fig. 5a). Small areas of extreme negative sensitivity were concentrated near the most intense NO_x emissions in XG (magenta borderline in Fig. 5), a heavy industrial park district in the study domain. Although high VOC emissions are also present in the industrial park, the magnitude of first-order sensitivity to NO_x emissions (−31 ppb h⁻¹) is relatively larger than that of VOC emissions (+18 ppb h⁻¹), indicating that the O₃ linear response to NO_x emissions is proportionally greater than that of VOC emissions from the industrial park (see Fig. 5a and b). Moving away from the source emissions, the negative sensitivity to daytime NO_x emission gradually increases and becomes positive in the inland area, indicating the shift of ozone production sensitivity from a VOC-limited to NO_x-limited regime. The negative sensitivities to daytime NO_x emissions extend to the coastal Pingtung county (green borderline in Fig. 5), reflecting the downwind transport of NO_x from the source region by the steering of northeasterly winds to westerly winds from the terrain effect and local circulation. Areas of positive sensitivity to VOC emissions overlap with areas of negative sensitivity to NO_x emissions, confirming the VOC-limited regime, but also extend to NO_x-saturated areas. The larger area of VOC sensitivities could also be explained by the longer atmospheric lifetime of VOCs (hours to days for high-OFP species such as toluene, ethylbenzene, xylene, ethane, and ethylene) as compared to NO_x (1.0–4.5 h) (Laughner and Cohen, 2019).

Figure 5 CMAQ-HDDM first-order sensitivity coefficient of O₃ to (a) NO_x emissions, (b) VOC emissions, second-order sensitivity coefficient of O₃ to (c) NO_x emissions, (d) VOC emissions, (e) second-order cross-sensitivity coefficients of O₃ to NO_x, VOC emissions, averaged during the entire simulation period at daytime, 09:00–15:00 LST. Magenta and green highlighted borderline represents the Xiaogang district and Pingtung region, respectively.

Second-order sensitivity indicates the responsiveness of ozone concentration to broader changes in emissions and also delineates the nonlinear sensitivities of ozone concentration to emissions. Most of the urban areas (except XG, an industrial park) exhibited negative second-order sensitivities to daytime NO_x emissions, which when coincident with negative first-order sensitivities yields a negative O₃ convex response: urban O₃ increases linearly but decreases nonlinearly with decreasing NO_x emissions. Thus, in these areas, the large negative second-order sensitivities countered the linear negative sensitivity of O₃ creating urban O₃ production conditions that were less sensitive than otherwise to changes in NO_x emissions. In contrast, positive second-order sensitivities to daytime NO_x emissions in XG were coincident with negative first-order sensitivities, which yielded a negative O₃ concave response with linear and nonlinear O₃ increases during decreasing NO_x emissions (Fig. S2). This concave O₃ response is also observed in the Chicago O'Hare (ORD) area, which is highly polluted and severely impacted by aircraft emissions (Arter and Arunachalam, 2021); it is similar to XG where it also hosted an international airport and several industrial parks. In XG, the NO_x titration effect that consumes ozone is relatively higher than other urban areas of Kaoping due to the high-NO_x$/$VOC conditions (see Fig. S6 in the Supplement), and hence a rapid reduction in NO_x level could further suppress the titration effect and increase the ozone formation rate, while for urban areas outside of XG, rapid NO_x reductions push the condition into transition or the NO_x-limited regime, and suppression of the titration is countered by less efficient O₃ production. Cross sensitivities indicate the interaction effects of the sensitivities when both NO_x and VOC emissions are changed simultaneously (limited to NO_x and VOC changes in the same direction, positive or negative). In most cases, decreasing the NO_x (VOC) emissions enhanced the sensitivities of ozone to the limiting precursors, NO_x (VOC), due to the positive cross sensitivity of ozone to NO_x and VOCs (see Fig. 5e). Cross sensitivities also indicate how results from emission control strategies may change from summing the results of decreasing individual NO_x and VOC emissions (Arter and Arunachalam, 2021). Cross sensitivities usually have only a small contribution in highly VOC-limited or highly NO_x-limited conditions because the associated net ozone production rates are driven singly either by VOC or NO_x, respectively (Sillman, 1999; Cohan et al., 2005). This result is reproduced in our model simulation; low cross sensitivities between NO_x and VOC are simulated in XG (highly VOC-limited condition) as well as remote areas in the west side of the mountain range in southern Taiwan (highly NO_x-limited condition), while high positive cross sensitivities between NO_x and VOC emissions are concentrated in downwind suburban areas (e.g., Pingtung county) (Fig. 5e). The latter result highlights a potential underprediction of O₃ formation sensitivity if the cross-sensitivity term is excluded in suburban areas.

3.2 Taylor-series expansion approximation

We then sought to exploit the sensitivity coefficients determined above to estimate the impact of NO_x and VOC emission reductions on ozone concentrations in southern Taiwan from the time of our study period, 2018 to 2025. To determine the level of emission reduction in the near term, we referred to the long-term trend of projected emissions in Taiwan as reported in the Taiwan Emission Data System (TEDS 11.0; https://air.epa.gov.tw, last access: 30 January 2022), which estimated an overall reduction of 53 % (NO_x) and 14 % (NMHC) over 22 years (2007–2028). In recent years, the decreasing trends for both NO_x and NMHC emissions have slowed down (see Fig. S7 in the Supplement). We then used 5 % emission reductions over the 8-year near-term period from 2018–2025 as a conservative estimate. Although this reduction is small, the Taylor-series approximation (Eqs. 4 and 5) considers the first-order, second-order, and cross-sensitivity coefficients of O₃ for this change in NO_x and VOC emissions and provides a clearer picture on how O₃ will change with a parallel reduction in emissions in the near future.

We estimated the ozone concentration at different scenarios represented by different emission control strategies for the typical seasonal transition period already simulated above, which often causes high-O₃ episodes in the study area. For simplicity, a reduction of domain-wide NO_x and VOC emissions is considered. Given the high nonlinearity of the ozone response to emissions of NO_x and VOC, the Taylor-series approximation is necessary to consider the higher-order terms involved in second-order and cross-term sensitivity of NO_x and VOC emissions. To compare a series of targeted emission control strategies, four experiments are considered: (1) baseline, (2) NO_x control, (3) VOC control, and (4) NO_x and VOC control. All experiments reduce the targeted emission by 5 % except for the baseline experiment representing the current emission with no perturbations. Figure 6 shows the baseline O₃ concentration and its corresponding response to NO_x control, VOC control, and NO_x and VOC control at daytime 12:00 LST. Benefits of NO_x control (i.e., reduction of O₃ by 0.6–0.8 ppbv) are seen to dominate over the inland area. However, an adverse effect of increased O₃ concentration >1.8 ppbv from NO_x control is simulated near the western coastal urban area where high NO_x emissions are concentrated. This agrees with our previous analysis that urban O₃ concentrations are estimated to increase from decreasing NO_x emissions for most of the urban areas that follow O₃ concave sensitivity behaviors and increase more for highly polluted urban areas (i.e., XG) that follow O₃ convex sensitivity behaviors. In contrast, no adverse effect is observed from VOC control but rather resulted in decreased O₃ concentration over much of the study region. The largest ozone response to VOC control is −1.4 ppbv in the urban area but also extended to the inland area. Control of both NO_x and VOC emissions resulted in reduced inland O₃ concentration, a marginal difference in urban areas, and slightly increased O₃ concentration near the large point source area in XG.

Figure 6 Spatial distribution of O₃ concentration in (a) the baseline with no perturbations in NO_x and VOC emissions, and changes in O₃ concentration under the (b) NO_x control scenario, (c) VOC control scenario, and (d) NO_x and VOC control scenario at daytime 12:00 LST. All scenarios reduced the targeted emissions by 5 % except for the baseline.

The decreasing trend of NO_x emission in the future is likely to continue owing to adoption of electric vehicles (EVs), which may impact the urban and inland areas differently due to the respective sensitivity regime. Using the Taylor-series approximation, we estimated ∂O₃ against ∂VOC emissions at different levels of projected ∂NO_x emissions for Xiaogang, XG (urban), and Chaozhou, CZ (inland), at 12:00 LST (Fig. 7). Projected ∂NO_x emissions at 0 % (black line) and −5 % (red line) represent the current year of 2018 and the near future of 2025, respectively, while other projected ∂NO_x emissions at arbitrary values of −50 % (blue), −25 % (magenta), and +25 % (green) are presented for alternative scenarios. At CZ, our results show that inland O₃ becomes less sensitive to varying ∂VOC emissions (i.e., fewer benefits of VOC control) at greater reductions of ∂NO_x emission at −25 % and −50 % (i.e., moving into the NO_x-limited regime). On the other hand, ∂O₃ at XG is always sensitive to varying levels of ∂VOC emissions at all levels of ∂NO_x scenarios, and ∂O₃ increases in parallel following the ∂NO_x order −50 %, −25 %, and −5 %. This is expected in a highly polluted area that favors VOC-limited conditions where O₃ could have the reverse effect to NO_x due to the reduced titration effect for decreasing NO_x. Besides, at all projected levels of ∂NO_x emissions, urban O₃ continues to decrease with decreasing VOC emissions. We showed in the Taylor-series approximation analysis when we hypothetically reduce NO_x emission at the arbitrary −50 % and −25 % scenario, which should expect an improvement in the overestimation of NO₂, urban O₃ at XG remains in a VOC-limited condition. Therefore, it is certain that the VOC-limited condition at urban areas of southern Taiwan is mainly due to the high anthropogenic NO_x emissions, and VOC controls are beneficial in reducing O₃ concentration in highly polluted urban areas.

Figure 7 Taylor-series approximation of ∂O₃ against ∂VOC emissions at different levels of projected ∂NO_x emissions for Xiaogang (urban) and Chaozhou (inland) averaged at 12:00 LST. Projected ∂NO_x emissions at 0 % (black line) and −5 % (red line) represent the current year of 2018 and the near future of 2025, respectively; other projected ∂NO_x emissions at arbitrary values of −50 % and ±25 % are presented for far-future comparison purposes.

3.3 Sensitivity of individual modeled VOC species

In this subsection, we further evaluate the sensitivity of the ozone response to individual modeled VOC emissions using the CMAQ direct decoupled method (DDM-3D). Note that only the first-order sensitivity coefficient of individual VOC species in CB6 is calculated, considering the high computational cost, and the second-order sensitivity coefficients of ozone to total VOC emissions have a low contribution. A total of 20 individual sensitivity tests were performed to target each modeled lumped VOC class or individual species (refer to Table 1 for a complete list of modeled VOC species). Figure 8a shows the daytime CMAQ-DDM sensitivity coefficient of ozone to the VOC emissions, arranged in ascending order and separated for the urban and inland area. The urban area sensitivity coefficients are substantially higher than that in the inland area for all species, which is expected due to the VOC-limited regime in the urban area. Among the 20 VOC species, the six most important VOC species were identified as XYL (xylene), OLE (terminal C–C olefins), PAR (paraffins), ETH (ethene), TOL (toluene), and IOL (internal C–C olefins).

Figure 8 (a) Daily averaged CMAQ-DDM first-order sensitivity coefficient of O₃ concentrations calculated per number of grids for each modeled VOC species arranged in ascending order for the urban and inland area. Sensitivity of O₃ and ozone formation potential (OFP) to (b, e) alkene emissions (OLE + ETH + IOL), (c, f) aromatic emissions (XYL + TOL), and (d, g) alkane emissions (PAR) at 12:00 LST.

XYL and TOL are aromatic hydrocarbons that are included in the BTEX group, along with benzene and ethylbenzene. The ranking of BTEX with respect to ozone formation potential in a high-NO_x environment is often done using maximum incremental reactivity (MIR), which uses unitless coefficients and indicates the amount contributed to ozone formation in the air mass by individual compounds. Based on the MIR scale adopted from the literature (Atkinson, 1997; Carter, 1994), m-xylene and p-xylene (8.20) are the most dominant contributors to ozone formation among BTEX, followed by toluene (2.70), ethylbenzene (2.70), and benzene (0.42). Our result agrees with the MIR scale, which yielded classifications of XYL and TOL as the second and fifth most sensitive VOC species to O₃ among the 20 VOC species considered in the sensitivity tests. Other sensitive VOC species to O₃ formation include alkanes (i.e., PAR (MIR=0.32)) and alkenes (i.e., OLE (8.24), ETH (4.37), IOL (13.11)).

To gain additional insights about the ozone response to the six most sensitive VOC species, we grouped them into three categories (i.e., alkenes: OLE, ETH, and IOL; aromatics: TOL and XYL; and alkanes: PAR) and plotted the spatial distributions of the sensitivities at 12:00 LST (Fig. 8b–d). The sensitivity of ozone to alkenes and alkanes has a similar pattern with higher sensitivity concentrated near the emission source and slightly shifted towards inland/mountainous areas and Pingtung county due to the local circulation and the prevailing northeasterly winds. For the sensitivity of ozone to aromatics, the spatial distribution has a smaller area coverage, and higher sensitivity is more concentrated near the emission source. To render a clearer picture on the spatial distribution of O₃ sensitivity for each VOCs component, the spatial distribution of some highly sensitive VOC component emissions (i.e., XYL, OLE, PAR, ETH, TOL, IOL) is provided in the Supplement (see Fig. S8). The role of land–sea-breeze circulation in transporting the local O₃ precursors (i.e., VOC) from the urban to inland area is well reflected by the diurnal profile of the sensitivity coefficients (Fig. 9). An obvious shift is simulated in the inland area as compared to the urban area; the urban sensitivity coefficients peak at 10:00 LST, while the peak in the inland area is at 13:00 LST. A lag of 2–3 h is mainly attributed to the daytime sea-breeze penetration that pushes the urban polluted air towards the inland area due to the differential heating between the land and sea surfaces during the noon hours. This result also tells us that VOC emissions from the urban area are largely responsible for the inland high-ozone episodes.

Figure 9 Diurnal variations of the CMAQ-DDM sensitivity of O₃ to the six most sensitive VOC modeled species (i.e., terminal olefins, OLE; xylene, XYL; paraffin, PAR; ethene, ETH; toluene; TOL; internal olefins, IOL) averaged for the (a) urban and (b) inland area, during the entire simulation period.

Using an approach similar to OFP, we quantified the relative ground-level ozone impacts of the six most sensitive VOC species by multiplying the MIR coefficients to its corresponding concentrations at 12:00 LST (see Fig. 8e–g). Due to the abundance of PAR, OFP calculated for alkane emission is significant over the western coastal region of southern Taiwan. The longer atmospheric lifetime of alkanes makes them more susceptible to transportation and impacting a wider area. Despite alkanes having the lowest sensitivity coefficient of the six groups of VOCs in the study, the significance of alkane OFP extends over a larger area, highlighting the need for controlling alkane-related emissions to reduce the ozone pollution problem over the study area. OFP of alkenes features a point source distribution pattern, indicating it is more likely related to large point source emissions. OFP of aromatics has the smallest area coverage and is highly concentrated near the populated industrial park of the study area in southern Taiwan. Among the three main contributors, aromatics have the highest maxima OFP value (56 ppbv), followed by alkenes (49 ppbv) and alkanes (40 ppbv).

Figure 10 shows the relative sensitivities for each modeled VOC to O₃ (DDM – Fig. 10a) and its corresponding emission rate (EMIS – Fig. 10b) averaged in urban and inland areas during the entire simulation period. Similar to the above analysis, the three main contributors were identified as alkenes (42.2 % for urban; 37.1 % for inland), aromatics (29.1 %; 25.1 %), and alkanes (18.2 %; 23.4 %). Other less important contributors were formaldehyde (4.1 %; 3.9 %), isoprene (2.3 %; 5.4 %), alcohols (1.9 %; 2.6 %), and aldehydes (0.9 %; 1.6 %). Due to the abundant biogenic emissions surrounding the inland area, inland ozone is more sensitive to isoprene (5.4 %) as compared to the urban area. Besides, inland ozone is also more sensitive to alcohols (MEOH + ETOH) and aldehydes (ALD2 + ALDX) when compared to the urban area. This is also consistent with the higher contributions of alcohols (10.0 %) and aldehydes (1.4 %) inland than in urban areas (4.3 % and 0.8 %, respectively) (Fig. 10b).

Figure 10 Relative distribution in percentage of (a) CMAQ-DDM O₃ sensitivity and (b) CMAQ-EMIS hydrocarbon emissions to individual VOCs summed across different groups averaged for the urban and inland area during the entire simulation period.

Overall, the DDM sensitivity tests identified three main speciated VOC groups contributing to ozone formation, namely alkenes, aromatics, and alkanes. The six most sensitive VOC surrogates in CB6 to ozone formation in the study area are XYL, OLE, PAR, ETH, TOL, and IOL in descending order. Sources of these targeted VOC species are highly variable, and previous studies have attributed these species to a wide range of emission sources not limited to traffic, industrial, petrochemical, manufacturing, solvent usage, etc. Therefore, in order to identify the sources of these VOC species for better emission control policy, source apportionment is necessary.

3.4 Descriptive statistics of PAMS-VOC data and PMF optimal solution

In this subsection, we first describe the statistics of the PAMS-VOC data measured at Chaozhou (CZ), Qiaotou (QT), and Xiaogang (XG) which were used to drive the PMF source apportionment. The CZ station is an inland station located far away from the urban core, while QT and XG are urban stations. The top 15 major VOC species measured at the three PAMS stations are summarized in Table 3, and other species are shown in Fig. S9 in the Supplement. These top 15 major VOC species accounted for 65.9 %, 70.7 %, and 68.7 % of the total VOC concentrations in CZ, QT, and XG, respectively. Overall, aromatics contributed the largest proportion, accounting for 117.3 ppbv (43.5 %) in CZ, 312.2 ppbv (53.9 %) in QT, and 361.7 ppbv (51.2 %) in XG. Toluene was the most abundant VOC species in all three stations, accounting for 46.3 ± 45.0 ppbv (17.2 % of total) in CZ, 140.4 ± 96.2 ppbv (24.3 %) in QT, and 143.0 ± 97.7 ppbv (20.2 %) in XG. Toluene is also commonly used in the chemical industry (Mo et al., 2015), the iron and steel industry (Nogueira et al., 2011), fuel evaporation (Wu et al., 2016), and organic solvent usage (Chen et al., 2019), and it is a by-product of vehicle exhaust (Dai et al., 2013). The second largest contributor to total VOC concentrations were alkanes, accounting for 42.7 % in CZ, 40.0 % in QT, and 40.6 % in XG. Among the alkanes, n-butane, propane, and isopentane comprised the highest proportions, each reaching around 5 % of the total VOC concentrations. n-Butane isomers are emitted from a variety of sources including petroleum production and natural gas emissions (Rossabi and Helmig, 2018). In the alkene category, ethylene and propylene are major contributors which are commonly used as raw materials in the petrochemical industry and refineries. Among the three stations, the highest proportion of alkene groups is observed in XG, accounting for 7.0 % of total VOC concentrations, with 3.6 % attributed to propylene.

Table 3 Concentration (mean ± SD) and proportions (%) of the top 15 PAMS-VOC species in ascending order at CZ, QT, and XG during 1–31 October 2018. Bold/italic species represent the unique species that present in the top 15 at each station. All units in parts per billion (ppb).

Determining the number of factors is a critical step in receptor-based source apportionment methods like PMF analysis. Combinations of three to eight sources were used to determine the optimum number of VOC sources for each PAMS station. We also calculated the uncorrelated bootstrap (BS) mapping to assist in identifying the optimum number of factors in each station. A high number of uncorrelated BS mapping indicates excessive factors are fitted in the model; therefore, this number should be kept as low as possible to avoid over-fitting. In CZ, the $Q_{\text{true}}/Q_{\text{expected}}$ was at a minimum in the six-factor solution, and the uncorrelated BS mapping was low (n=3), indicating this was the optimal solution for that site (Fig. S10 in the Supplement). Using the same protocol, the optimal solution was determined as seven factors for QT and five factors for XG. Overall, a total of eight unique factors were identified at the three stations: biogenic sources, solvent usage, vehicle emissions, plastic industry, manufacturing industry, mixed industry, aged air mass, and motorcycle exhausts (Fig. 11). The diurnal variations of the PMF factor contribution at Chaozhou, Qiaotou, and Xiaogang PAMS stations are also shown in Fig. 12. Among the eight resolved factor profiles, vehicle emissions (22 %) were the largest contribution to total VOC concentrations, followed by mixed industry (21 %), solvent usage (17 %), biogenic sources (12 %), plastic industry (10 %), and other factor profiles (i.e., aged air mass, motorcycle exhausts, and manufacturing industry). The diurnal pattern of factor distribution of the biogenic sources displays a clear diurnal cycle peak at noontime, 10:00–12:00 LST, in all three stations. Meanwhile, the motor exhaust and vehicle emission sector peak during the morning and evening traffic rush hours. The aged air mass sector is only identified in CZ, and its hourly factor contribution mostly existed at more stable levels when compared to other factors except for an obvious peak at 12:00 LST due to the sea-breeze penetration that pushes the urban polluted aged air mass towards inland area. The diurnal pattern related to industrial activity exhibits different peak hours depending on sector and station. For instance, the hourly factor contribution of solvent usage in XG has a clear bimodal peak at 10:00 and 16:00 LST, plastic industry in QT and CZ has a clear diurnal cycle peak in around noontime at 11:00 and 13:00 LST, and mixed industry in QT and XG has a clear bimodal peak at 07:00 and 18:00 LST. Details of each source profile and comparison with other PMF studies are discussed in the Supplement – “Source Profiles of PMF Model”.

Figure 11 Source profiles calculated using the PMF model at Chaozhou, Qiaotou, and Xiaogang PAMS stations.

Figure 12 Diurnal variations of factor contribution calculated by the PMF model at Chaozhou, Qiaotou, and Xiaogang PAMS stations.

Figure 13 Normalized factor contribution and composite index of source profile to each PAMS-VOC species at Chaozhou, Qiaotou, and Xiaogang PAMS station.

3.5 Dominant sources of highly sensitive VOC species

We then mapped the CMAQ-CB6-modeled VOC species to the PMF apportioned sources in order to identify the dominant sources of highly sensitive VOC species (i.e., XYL, OLE, PAR, ETH, TOL, IOL). The dominant sources of alkene species (OLE, ETH, IOL) were identified using ethylene and propylene as tracers; sources of TOL were tracked by toluene and m-ethyltoluene; sources of XYL were tracked by all three xylene isomers; and sources of PAR were identified using isobutane, n-butane, isopentane, n-pentane, n-hexane, and 2,2,4-trimethylpentane.

We also developed a composite index, ψ, to quantitatively combine the CMAQ-DDM sensitivity coefficient and PMF-resolved factor contribution in order to identify the key VOC species for an effective ozone abatement strategy. The index was calculated according to the jth species and kth source:

$$\begin{array}{}\text{(9)}& {\mathit{\psi }}_{jk}={\displaystyle \frac{\partial C}{\partial {\mathit{\epsilon }}_j}}\times f_{jk},\end{array}$$

where C is the trace gas (O₃) concentration, ε is the unitless scaling factor for the j emission (see Eqs. 2 and 3), f_kj is the jth species concentration to the kth source. The first term in Eq. (9) is obtained from the CMAQ-DDM-calculated first-order sensitivity coefficient where a higher value denotes greater sensitivity (extreme low sensitivity <0.01 ppb is masked out); the second term is the PMF-resolved factor contribution where higher values of a particular species, j, indicate greater apportionment to the source, k. Higher values of the composite index indicate a greater priority should be given to that source species.

According to the PMF source apportionment (Fig. 13), the dominant source of alkene-related emissions (OLE, IOL, ETH) was attributed to mixed industry (i.e., petrochemical industry, printing industry, and metal industry) (Chen et al., 2022; Pinthong et al., 2022). This source profile contributed 61.5 % propylene and 47.6 % ethylene in XG, as well as 52.4 % propylene and 37.9 % ethylene in QT. Meanwhile, propylene and ethylene in CZ were mainly attributed to vehicle emissions and aged air masses. Given that alkene-related emissions are closely linked to large point sources (see Fig. 8e), control of alkene emissions should focus on the mixed industry, particularly petrochemical plants. Among the three stations, the highest composite index associated with ethylene and propylene summed across all factors was obtained at XG. Though ethylene and propylene have higher sensitivities compared to other species, their composite index (ψ_{ETHY, all k}=5.7 and ψ_{PROP, all k}=7.4 – see Fig. 13) is relatively low due to their low abundance, 11.5 ppb (1.6 % of total VOCs) and 25.5 ppb (3.6 %), respectively, making them less important in terms of priority.

The highest composite index was assigned to aromatic compounds such as toluene (ψ_{TOL, all k}=12.1, 16.4, and 53.5 in CZ, QT, and XG, respectively), m,p-xylene (${\mathit{\psi }}_{m,p-\mathrm{XYL},\mathrm{all}k}=\mathrm{12.0}$, 19.7, and 53.5), and o-xylene (ψ_{o-XYL, all k}=4.5, 6.7, and 19.0) due to their high-abundance toluene, 46.3–143.0 ppb (17.2 %–24.3 % of total VOCs), m,p-xylene, 13.6–63.8 ppb (5.0 %–9.0 %), and o-xylene, 5.2–23.3 ppb (1.9 %–3.3 %), in all three stations (see Fig. 13 and Table 3). Although toluene has much higher abundance than xylene, the composite index for both species summed across all sources is similar because the sensitivity coefficient of xylene ($\partial C/\partial \mathit{\epsilon }>\mathrm{12}$ ppb) is relatively higher than toluene ($\partial C/\partial \mathit{\epsilon }>\mathrm{3}$ ppb) (see Fig. 8a). This indicates that emission control of both of these compounds is of high priority considering their high abundance and high sensitivity in ozone formation. The dominant sources of aromatic-related emissions (XYL, TOL) are mainly attributed to solvent usage and mixed industry. Solvent usage, including building coatings, paint thinners, and other products thinners (Li et al., 2021; Chen et al., 2019; Wu et al., 2016; Huang and Hsieh, 2020), contributed 64.5 % xylene and 32.2 % toluene in CZ, 40.6 % xylene and 29.7 % toluene in QT, and 38.5 % xylene and 43.3 % toluene in XG. Mixed industry contributed 26.2 % xylene and 21.9 % toluene in QT, as well as 43.5 % xylene and 45.3 % toluene in XG. Vehicle emissions, including motorcycle exhaust, contributed less than 10 % of xylene and toluene in CZ and QT. Therefore, emission control of XYL and TOL should target solvent usage, particularly painting, coating, synthetic fragrances, adhesives and cleaning agents, and industrial sources in the region.

The dominant source of alkane-related emissions (i.e., PAR) is attributed to vehicle emissions in all three stations. The main tracers identified in the PMF source apportionment related to vehicle emissions are isobutane, n-butane, isopentane, and n-pentane; these compounds accounted for 54.9 %, 53.0 %, 47.0 %, and 47.0 % of the PMF normalized factor contribution in CZ; 54.8 %, 58.7 %, 25.9 %, and 29.1 % in QT; and 44.5 %, 43.8 %, 23.9 %, and 26.8 % in XG, respectively. Motorcycle exhaust also contributed to these compounds, but the average contribution was relatively lower (e.g., <20 % in QT). High amounts of C4–C5 alkanes (i.e., isobutane, n-butane, isopentane, n-pentane) are known indicators of traffic-related sources (S. Yu et al., 2021; Huang and Hsieh, 2020; Chen et al., 2019). Among these compounds, n-butane had the highest composite index due to its high abundance, 15–34 ppb (5.1 %–5.5 % of total VOCs). Other than freshly emitted vehicle emissions, aged air masses also contributed a significant amount (∼30 %) to these compounds at CZ, reflecting the role of land–sea-breeze circulation transporting urban vehicle plumes inland. In addition, acetylene is a well-known indicator for combustion sources and related to liquefied petroleum gas (LPG) leakage. The ratio of propane, n-butane, or i-butane to acetylene is often used to assess the domination of the LPG-vehicle sources across different cities (An et al., 2014). These ratios were 0.7, 0.8, and 0.4 in CZ; 6.0, 8.2, and 3.9 in QT; and 2.9, 3.8, and 1.7 in XG, respectively, which were much lower than 11.5, 1.8, and 2.6 in Guangzhou (Zhang et al., 2013) and 11.4, 5.0, and 2.3 in Mexico City (Blake and Rowland, 1995), which are heavily impacted by LPG-vehicle emissions. This indicates that there were relatively fewer LPG-fueled vehicles in southern Taiwan, whereas LPG is merely used to fuel a gas stove for domestic household use or catering. Therefore, emission control of PAR species should focus on the gasoline-fueled vehicle emissions, particularly in heavy-traffic cities.

4 Conclusion

In this work, we used the CMAQ-DDM-3D sensitivity tool in conjunction with PMF model source apportionment to provide a comprehensive analysis of the major contributors to VOC species and the ozone formation potential (OFP) over an area of southern Taiwan. We developed a composite index that quantitatively combines the CMAQ-DDM sensitivity coefficient and PMF-resolved factor contribution to identify the key VOC source species for the effective ozone abatement strategy, which is applicable to other urban areas that are characterized by the VOC-limited condition. A representative case in October 2018 with a daily 8 h maxima O₃>75 ppb occurring often was the focus of this study. Low NO_x levels and high biogenic volatile organic compound (BVOC) levels in the inland areas yielded a NO_x-limited regime. Although reducing NO_x emissions can reduce the inland O₃ concentrations, it could adversely increase urban O₃ due to the reduced titration effect. This is because the O₃ sensitivity production in the urban area is mainly dominated by VOC-limited regime, and further reducing NO_x emissions can suppress the titration effect and eventually increase the urban O₃ concentration. In contrast, control of VOC emissions is effective in the urban area to reduce O₃ concentrations and has no adverse effect in the rural area. Our DDM sensitivity analysis identified the six most sensitive VOC species or groups: XYL (xylene), OLE (terminal olefins), PAR (paraffins), ETH (ethene), TOL (toluene), and IOL (internal olefins) in descending order, which are mainly attributed to the petrochemical industry, painting and printing industry, and vehicle emissions. Based on a composite index, effective ozone abatement strategies should prioritize (1) solvent usage such as the painting, coating, and printing industry that emits abundant toluene and xylene; (2) petrol vehicular emissions with high compositions of n-butane, isopentane, isobutane, and n-pentane; and (3) the petrochemical industry with emphasis on ethylene and propylene.

Besides, our DDM results also highlighted the important role of formaldehyde in OFP over the study area given that it is ranked the seventh highest VOC species for O₃ sensitivity and on a per-molecule basis for O₃ sensitivity is similar to alkene (i.e., OLE, ETH, IOL) and aromatic (XYL, TOL) species. Mapping of formaldehyde to the PMF source apportionment was not done due to the lack of oxygenates in the PAMS measurement inventories. Thus, oxygenates such as FORM, ALD2, ALDX, ACET, and KET along with alcohols (ETOH and MEOH) should be included in the PMF source apportionment for more complete source identification in future studies. In the current work, WRF and CMAQ are resolved at a high resolution of 1 km to best represent the features of local circulations (i.e., land–sea breeze, urban heat island effect) at the urban scale. However, simulation at 1 km is obviously too expensive for large domain or regional modeling. The differences in HDDM and PMF analysis between D03 (3 km) and D04 (1 km) due to grid resolution remain an open question which deserves future in-depth investigations. Although the performance of the simulated meteorological parameters (T2, WS, and WD) at both urban and rural stations is acceptable in the benchmark recommended by the US EPA, notable differences in temperature (underestimation) and wind speed (overestimation) are still observable in our simulation work. These biases could be susceptible to underestimation in photochemical ozone production due to the fictitious cold bias and enhanced dispersion. Therefore, careful treatment on the urban-scale data assimilation in temperature, wind field, and relative humidity are recommended in the future to improve the model prediction.

Other important findings obtained in this study are summarized below.

• Negative (positive) first-order sensitivities to daytime NO_x emissions are dominant in the urban (inland) area, indicating ozone production sensitivity favors the VOC-limited regime (NO_x-limited regime). Negative sensitivities are also extended to some parts of the coastal area of Pingtung county, reflecting the downwind transport of the urban NO_x by the steering northeasterly winds due to the terrain effects and local circulations.

• Most of the urban areas (except Xiaogang, an industrial park) exhibited negative first-order and second-order sensitivities to daytime NO_x emissions, indicating a negative O₃ convex response. However, due to the large negative second-order sensitivities, the urban O₃ increases from the linear effect are largely attenuated by the nonlinear effect. As a result, urban O₃ becomes less sensitive to changes in NO_x emissions but more sensitive to VOC emissions, which favors the VOC-limited conditions.

• Based on the PMF model source apportionment, a total of eight factor profiles are identified as mixed industry (21 %), vehicle emissions (22 %), solvent usage (17 %), biogenic sources (12 %), plastic industry (10 %), aged air mass (7 %), motorcycle exhausts (7 %), and manufacturing industry (5 %).

• The benefit of VOC control in inland areas is expected to decrease gradually when NO_x emissions continue to decrease over the long-term, but control of VOCs in highly polluted urban areas remains effective despite a large reduction in NO_x in the future.