Evaluating simulations of organic aerosol volatility and degree of oxygenation in eastern China

Volatility and oxygen-to-carbon ($\mathrm{O}/\mathrm{C}$) molar ratios are critical properties of organic aerosols (OA), influencing their viscosity, hygroscopicity, and light absorption thereby resulting in impacts on air quality and climate. While atmospheric models often track these properties to simulate OA evolution, their performance remains insufficiently evaluated. This study assessed OA volatility and $\mathrm{O}/\mathrm{C}$ simulations by comparing CMAQ model outputs using official AERO7i and community-contributed two-dimensional volatility basis set (2D-VBS) schemes, against two field measurements in eastern China. Apart from baseline modelling, two additional simulations using AERO7i incrementally incorporated low-volatility/semi-volatile/intermediate-volatility organic compound (L/S/IVOC) emissions and enhanced anthropogenic secondary organic aerosol (SOA) yields. An optimized 2D-VBS simulation further constrained $\mathrm{O}/\mathrm{C}$ ratios of primary organic aerosol (POA) emissions using observations. The results showed that OA mass concentrations were underestimated by 24 % in 2D-VBS and 27 %–34 % with updated AERO7i, likely due to underrepresented vehicular POA emissions and nighttime SOA formation. All simulations captured the substantial contribution of low-volatility products ($C^{*}<\mathrm{0.1}$ µg m⁻³) but failed to reproduce the detailed volatility distributions within this range. Simulated $\mathrm{O}/\mathrm{C}$ ratios were biased low in aged air masses (notably with 2D-VBS) and slightly overestimated in areas with more local emissions using updated AERO7i. Misrepresentations of OA volatility significantly led to biases in viscosity predictions, while the hygroscopicity parameter played a more important role. These findings highlight the need to better constrain OA volatility and $\mathrm{O}/\mathrm{C}$ in models to improve projections of OA air quality and climate impacts.

1 Introduction

Organic aerosols (OA) are a major component of fine particulate matter (PM_2.5), accounting for 20 %–90 % of its mass (Jimenez et al., 2009; Zhang et al., 2007), and play a critical role in global air quality, climate, and public health (Nault et al., 2021; Shrivastava et al., 2017; Wu et al., 2018). OA can be classified as either primary organic aerosol (POA), emitted directly from sources such as combustion, or secondary organic aerosol (SOA), formed through the oxidation of volatile organic compounds (VOCs) and semi-volatile and intermediate volatility organic compounds (S/IVOCs) (Hallquist et al., 2009; Heald and Kroll, 2020). The chemical aging of OA involves functionalization (incorporation of oxygen-containing groups), fragmentation, and oligomerization. These processes alter OA's composition and physicochemical properties, including volatility, oxidation state, viscosity, hygroscopicity, and light absorption (Chacon-Madrid and Donahue, 2011; Tritscher et al., 2011; Massoli et al., 2010; Rothfuss and Petters, 2017; Hems et al., 2021). Understanding these evolving properties is crucial for accurately predicting the impacts of OA on climate, air quality, and public health in atmospheric and climate models (Tsigaridis and Kanakidou, 2018).

Among the various properties of OA, volatility and the oxygen-to-carbon molar ratio ($\mathrm{O}/\mathrm{C}$, an indicator of the extent of oxygenation) are pivotal for constraining their atmospheric fate and impacts, and they are therefore important parameters in chemical transport models (CTMs) (Donahue et al., 2013; Rao and Vejerano, 2018). The two-dimensional volatility basis set (2D-VBS), a widely adopted framework in OA modeling, tracks volatility and $\mathrm{O}/\mathrm{C}$ to simulate chemical evolution (Donahue et al., 2011). Volatility primarily governs the partitioning of organic compounds between particulate and gas phases at equilibrium, thereby influencing their atmospheric behavior (e.g., chemical aging, transport, and deposition) and overall OA concentrations (Donahue et al., 2014; Shiraiwa and Seinfeld, 2012). Li et al. (2020) recently parameterized the glass transition temperature (T_g, at which a phase transition between amorphous solid and semi-solid states occurs, and viscosity changes dramatically) as a function of volatility and $\mathrm{O}/\mathrm{C}$. This parametrization highlights the link between volatility and phase state (or viscosity), which influence the kinetics of gas-particle interactions, with implications for diffusion, partitioning, and heterogeneous reaction rates (Zaveri et al., 2014; Reid et al., 2018; Marshall et al., 2018; Li and Shiraiwa, 2019). Furthermore, OA volatility affects cloud condensation nuclei (CCN) activity due to its connection with particle hygroscopicity and growth dynamics (Liu and Matsui, 2022; Zhang et al., 2023).

The $\mathrm{O}/\mathrm{C}$ ratios of OA influence their viscosity and hygroscopicity, similar to the effects of volatility (Massoli et al., 2010; Koop et al., 2011). The $\mathrm{O}/\mathrm{C}$ ratios and volatility of OA components are largely coupled, with more volatile components typically exhibiting lower $\mathrm{O}/\mathrm{C}$ ratios. Lower $\mathrm{O}/\mathrm{C}$ ratios favor phase separation under specific relative humidity conditions (Pye et al., 2017). Conversely, elevated $\mathrm{O}/\mathrm{C}$ ratios indicate increased hygroscopicity resulting in enhanced water uptake, thereby increasing their potential to act as cloud condensation nuclei and ice nuclei (Mahrt et al., 2022; Malek et al., 2023; Song et al., 2012; Tian et al., 2022). Furthermore, $\mathrm{O}/\mathrm{C}$ ratios substantially impact the optical properties of OA (Xu et al., 2024). During aging, bulk OA generally darkens with increasing $\mathrm{O}/\mathrm{C}$, while excessive oxidation at higher $\mathrm{O}/\mathrm{C}$ levels diminishes light absorption (Jiang et al., 2022; Duan et al., 2024). These findings suggest that $\mathrm{O}/\mathrm{C}$ could serve as a critical parameter in radiative forcing estimation in climate models.

A few prior studies have revealed significant uncertainties in OA volatility and $\mathrm{O}/\mathrm{C}$ simulations with CTMs. For instance, Saha et al. (2017) demonstrated that the VBS module in WRF-Chem significantly underestimated OA concentrations, and failed to reproduce low-volatility OA components with effective saturation concentration (C^∗) between 10⁻⁴ and 10⁻¹ µg m⁻³. This discrepancy could be partially explained by underestimated SOA formation due to wall losses of condensable vapors and missing low-volatility products, with the parameters empirically derived from aerosol growth experiments. In contrast, parameters derived from dual thermodenuder (TD) measurements predicted 2–4× higher SOA yields, and produced more low-volatility products under atmospherically relevant conditions (Saha and Grieshop, 2016). Additionally, the lack of S/IVOC emissions in the model likely contributed to biases in simulated OA volatility distributions (Xu et al., 2019). Regarding $\mathrm{O}/\mathrm{C}$ simulations, Tsimpidi et al. (2018) introduced the ORACLE 2-D module into a global chemistry-climate model, which tended to overpredict OA $\mathrm{O}/\mathrm{C}$ ratios in urban downwind areas. Overall, the model exhibited a 5 %–7 % overestimation of $\mathrm{O}/\mathrm{C}$ for OA and SOA, with the most pronounced positive biases in summer. The regional model Community Multiscale Air Quality (CMAQ) also overestimated $\mathrm{OM}/\mathrm{OC}$ ratio (closely related to $\mathrm{O}/\mathrm{C}$) of OA compared to observations in the southeastern US (Pye et al., 2017). However, accounting for interactions between OA and aerosol water, which enhance semi-volatile partitioning by increasing the available partitioning medium, could reduce model biases in $\mathrm{OM}/\mathrm{OC}$. The air quality model Polyphemus underestimated $\mathrm{O}/\mathrm{C}$ in the northwestern Mediterranean region, even after implementing multi-generational oxidation processes, highlighting insufficient representation of aging processes in the model (Chrit et al., 2018). While model performance varies across different regions and models globally, simulations of OA volatility and $\mathrm{O}/\mathrm{C}$ ratios in the polluted atmospheres of China remain insufficiently explored.

In this study, we evaluated the performance of the CMAQ model in simulating OA mass concentration, volatility distribution and $\mathrm{O}/\mathrm{C}$ ratio, by comparing to the observations at two sites in eastern China. The evaluation focused on two OA schemes that differed in their level of complexity regarding the representation of volatility and $\mathrm{O}/\mathrm{C}$. Additionally, two adjustments to the default CMAQ scheme were evaluated: (1) adding S/IVOC emissions with SOA formation from these precursors, and (2) updating the SOA yields of aromatic species (benzene, toluene, xylene) and polycyclic aromatic hydrocarbons (PAHs). The assessment also explored how uncertainties in OA volatility and $\mathrm{O}/\mathrm{C}$ propagated to T_g and viscosity predictions. These efforts aimed to provide valuable insights into improving OA simulations (in terms of both mass concentrations and properties), and their implications for viscosity and phase state modeling in CTMs.

2 Methods and Data

2.1 Model configuration

In this study, two versions of the CMAQ model (https://epa.gov/cmaq, last access: 6 November 2023), the official v5.3.2 (Appel et al., 2021) (available from https://doi.org/10.5281/zenodo.4081737, US EPA Office of Research and Development, 2020) and v5.4 with a community contribution (available from https://doi.org/10.5281/zenodo.7218076, US EPA Office of Research and Development, 2022), were used to perform nested simulations with horizontal resolutions of 36 and 12 km. The outer domain encompassed most of China, while the inner domain focused on eastern China (Fig. S1 in the Supplement). The simulation period covered 17 March–21 April and 29 September–21 November 2018, aligning with the campaigns conducted in Dongying (DY) and Guangzhou (GZ), respectively. Meteorological fields were generated using the Weather Research and Forecasting (WRF) model version 4.2.1. Anthropogenic emissions for China were represented using the 2018 high-resolution (0.25° × 0.25°) Multi-resolution Emission Inventory for China (MEIC) v1.4 (http://www.meicmodel.org, last access: 9 September 2023), while the Regional Emission Inventory in Asia (REAS) v3.2.1 (https://www.nies.go.jp/REAS/, last access: 12 September 2023) was used for the rest of Asia. Emissions from open burning were obtained from the Fire Inventory from NCAR (FINN) v1.5 (https://www2.acom.ucar.edu/modeling/finn-fire-inventory-ncar, last access: 13 January 2021), and biogenic emissions were calculated using the Model of Emissions of Gases and Aerosols from Nature (MEGAN) v2.1 (Guenther et al., 2012). The simulations utilized the SAPRC07tic gas-phase mechanism within CMAQ (Xie et al., 2013). The capabilities of WRF and CMAQ in simulating meteorological factors and major pollutants (NO₂, SO₂, O₃ and PM_2.5) were evaluated (Tables S1–S2 in the Supplement).

2.2 OA representations

All the simulations modeled SOA formation from isoprene, glyoxal, and methylglyoxal, as well as NO₃-initiated oxidation of monoterpenes, using consistent parameterizations from CMAQ's official AERO7i module (see Fig. S2). Specifically, isoprene SOA formed via the aqueous uptake of isoprene epoxydiols (IEPOX)/methacrylic acid epoxide (MAE), and monoterpene-derived organic nitrates (ONs) via NO₃ oxidation were explicitly represented (Pye et al., 2013; Pye et al., 2015). Semi-volatile isoprene SOA, including contributions from non-aqueous pathways, was parameterized with the two-product model (Carlton et al., 2010), which utilized mass yields of two products (α_i, i=1, 2) and their effective saturation concentrations ($C_i^{\ast }$, $i=\mathrm{1},\mathrm{2}$) to simulate aerosol growth observed in experiments (Pankow, 1994). All the simulations accounted for multiphase SOA formation from glyoxal and methylglyoxal both in cloud (Carlton et al., 2008) and on wet aerosol surfaces (Pye et al., 2015).

Five simulations were conducted to evaluate model performance in reproducing OA mass concentrations, volatility distributions, and $\mathrm{O}/\mathrm{C}$ ratios against observational data (Table 1). All simulations retained the same AERO7i heterogeneous chemistry and selected other systems (see above), and focused on different representations of key volatility-based systems. The simulation using the standard AERO7i treatment including SOA from anthropogenic precursors (e.g., aromatics, long-chain alkanes) (Qin et al., 2021), monoterpenes via ${\mathrm{O}}_{\mathrm{3}}/\mathrm{OH}$ oxidation (Xu et al., 2018), and sesquiterpenes (Carlton et al., 2010) was referred to as a 1-D VBS for simplicity. In another simulation, the 2D-VBS framework (Zhao et al., 2015; Zhao et al., 2016; Chang et al., 2022) was used for anthropogenic VOCs, monoterpenes (excluding NO₃ oxidation pathways), and sesquiterpenes. Additionally, in 2D-VBS, the oxidation of POA (or L/S/IVOCs) was represented differently compared to the 1D-VBS simulations (Murphy et al., 2017).

Table 1 Description of the five simulations conducted in this study.

The 1D-VBS simulation reflected the default CMAQ configuration without modifications to emissions or chemistry. The 1D-VBS_E simulation incorporated additional IVOC emissions, and updated emissions of S/LVOCs (historically classified as POA), based on MEICv1.4 (Sect. S1 in the Supplement). L/SVOC and IVOC emissions in each volatility bin were scaled from POA and total VOC emissions, respectively, using source-specific scaling factors (Table S3). Since the lowest volatility bin for emissions in the standard 1D-VBS parameterization was set to $C^{\ast }={\mathrm{10}}^{-\mathrm{1}}$ µg m⁻³, the difference between original POA emissions and L/SVOC emissions within the model-resolved volatility range ($-\mathrm{1}\le \log C^{\ast }\left(\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}\right)\le \mathrm{2}$) was reclassified as non-volatile ($C^{\ast }={\mathrm{10}}^{-\mathrm{10}}$ µg m⁻³, see Fig. S3) and treated with heterogeneous aging chemistry (Simon and Bhave, 2012). This adjustment conserved POA mass while addressing the model's limitations in volatility coverage. The estimated nation-wide L/SVOC and IVOC emissions, which were 3.18 and 6.68 Tg yr⁻¹, respectively, were higher than those reported by Zheng et al. (2023), but agreed well with Chen et al. (2024) in magnitude (Table S4, Figs. S4–S5), as both this study and Chen et al. (2024) applied a ratio-based methodology combined with the MEIC emission inventory. In contrast, Zheng et al. (2023) employed a different approach, using emission factors and activity data obtained from ABaCAS (Air Benefit and Cost and Attainment Assessment System) and Chang et al. (2022), which led to larger discrepancies relative to our estimates. The source contributions and volatility distributions of L/S/IVOCs differed slightly between studies. For example, our estimates indicated a lower contribution of solvent use to IVOCs compared with Zheng et al. (2023) (40 % vs. 57 %), and a higher contribution from residential sources to S/LVOCs (49 % vs. 30 %). In addition, the differences in emission magnitudes were primarily in the IVOC volatility range, whereas discrepancies in L/SVOCs were smaller. In the 1D-VBS_E simulation, the aging of L/SVOCs followed the POA treatment established in prior studies (Donahue et al., 2012; Murphy et al., 2017), where a fraction of L/SVOCs was oxidized in the gas phase to form SOA and the volatility distribution continually updated. The IVOC-derived SOA formation adopted the parameterization of Lu et al. (2020). In the 1D-VBS_EY simulation, SOA yields for benzene, toluene, xylene, and naphthalene oxidation were updated (Table S5) to account for vapor wall losses and formation of highly oxygenated organic molecules (HOMs) via autoxidation (Bilsback et al., 2023).

The 2D-VBS scheme (Zhao et al., 2015; Zhao et al., 2016; Chang et al., 2022), was implemented into CMAQv5.4 and evaluated using emission inputs largely consistent with the 1D-VBS_E simulation. The major differences in emissions included (1) reallocating non-volatile POA to the LVOC species with C^∗ of 10⁻² µg m⁻³ (i.e., the least volatile primary emission category in 2D-VBS); and (2) mapping SVOC emissions with C^∗ of 10⁻¹ µg m⁻³ to the $C^{\ast }={\mathrm{10}}^{-\mathrm{2}}$ µg m⁻³ bin as well (due to the absence of a corresponding volatility bin in the 2D-VBS scheme). L/SVOC aging products were classified as SOA. Initial $\mathrm{O}/\mathrm{C}$ ratio distributions for L/S/IVOC emissions across volatility bins followed 2D-VBS default profiles. However, these default assumptions led to underestimated POA $\mathrm{O}/\mathrm{C}$ ratios (see Sect. 3.3). Therefore, we conducted the 2D-VBS_A simulation, wherein sector-specific $\mathrm{O}/\mathrm{C}$ ratios for POA emissions from gasoline/diesel vehicles, non-road mobile sources, power plants and industrial sources were constrained using observational averages in eastern China (Table S6).

2.3 Volatility distribution and $\mathrm{O}/\mathrm{C}$ ratio modelling

The volatility (expressed in terms of C^∗) and $\mathrm{O}/\mathrm{C}$ ratios for OA surrogate species across all the simulations were summarized in Tables S7–S8. Volatilities for all species, except those in VBS-based parameterizations, were aligned with surrogate data in chamber experiments. Non-volatile POA and SOA formed via pathways such as aromatic oxidation under low-NO_x conditions, aqueous-phase uptake of IEPOX/MAE/glyoxal/methylglyoxal, ON hydrolysis and particle-phase oligomerization were assigned a fixed C^∗ of 10⁻¹⁰ µg m⁻³ for simplicity. The volatility of primary L/S/IVOCs in the condensed phase and aged OA, was described in detail by Murphy et al. (2017). In the 1D-VBS_EY simulation, a new SOA species with $C^{\ast }=\mathrm{0.1}$µg m⁻³ (unrepresented in the base 1D-VBS in CMAQ) was explicitly incorporated. The volatility distribution of total SOA (POA) was then derived using simulated concentrations of individual SOA (POA) species.

The overall $\mathrm{O}/\mathrm{C}$ ratios for SOA, POA, and OA were calculated using the mole-weighted averages of the $\mathrm{O}/\mathrm{C}$ values for individual species. The $\mathrm{O}/\mathrm{C}$ ratios of OA species were calculated based on mass-based $\mathrm{OM}/\mathrm{OC}$ ratios provided with the CMAQ code. The $\mathrm{OM}/\mathrm{OC}$ ratios for SOA were primarily derived from chamber experiments (Pye et al., 2017; Carlton et al., 2010; Xu et al., 2018). For example, SOA species with relatively well-known structures (e.g., IEPOX-/MAE-derived SOA, isoprene dinitrates and monoterpene nitrates), along with seven VBS bins representing monoterpene SOA formed via O₃ and OH oxidation, adopted $\mathrm{OM}/\mathrm{OC}$ values and other molecular properties from surrogate compounds identified in chamber studies. For other SOA species, once their C^∗ values were determined (as described above) and the number of carbons (nC) was assumed, the $\mathrm{OM}/\mathrm{OC}$ ratios were either inferred based on plausible structures (Pankow et al., 2015) or estimated using the relationship between volatility, nC and the number of oxygens (nO) used in 2D-VBS ((Pye et al., 2017). Additionally, $\mathrm{OM}/\mathrm{OC}$ ratios for primary L/S/IVOCs and their aging products were constrained by laboratory work and filed observations as documented in previous studies (Murphy et al., 2017; Lu et al., 2020). The non-volatile POA emissions were assumed to have an $\mathrm{OM}/\mathrm{OC}$ ratio of 1.6 (Turpin and Lim, 2001), which was equivalent to $\mathrm{O}/\mathrm{C}$ of 0.35.

In the base 2D-VBS, the $\mathrm{O}/\mathrm{C}$ distributions for L/S/IVOC emissions in each volatility bin followed the default settings in the community-contributed 2D-VBS. These $\mathrm{O}/\mathrm{C}$ ratios were primarily based on Chang et al. (2022), which used emission test data for specific sectors when available and otherwise adopted values from prior studies. However, in the 2D VBS, the $\mathrm{O}/\mathrm{C}$ ratios for emissions were capped at 0.4. For SOA products overlapping with AERO7i, such as isoprene SOA, the $\mathrm{O}/\mathrm{C}$ ratios were set to match those in AERO7i. For other SOA surrogate species, the $\mathrm{O}/\mathrm{C}$ ratios were determined either by adjusting the $\mathrm{O}/\mathrm{C}$ distribution of first-generation products to align with experimental data (OH/O₃-intitiated monoterpene and sesquiterpene SOA), or based on explicit chemical mechanisms for initial oxidation (anthropogenic SOA), followed by aging through functionalization and fragmentation within 2D-VBS (Zhao et al., 2015). As a result, the $\mathrm{O}/\mathrm{C}$ ratios for OH/O₃-intitiated monoterpene and sesquiterpene SOA in 2D-VBS generally ranged from 0 to 1.0 in 0.1 increments. In contrast, anthropogenic SOA (ASOA) exhibited a broader range, with a maximum $\mathrm{O}/\mathrm{C}$ ratio of 2.0, to account for the high degree of oxygenation observed in toluene SOA in chamber experiments.

2.4 Predictions of T_g and viscosity

Previous measurements have demonstrated a close relationship between volatility and viscosity. For instance, Champion et al. (2019) found that SOA with higher fractions of EL/LVOCs showed increased viscosity. Similarly, an inverse correlation between T_g and vapor pressure was observed for isoprene SOA components (Zhang et al., 2019). Although the dependence of T_g on atomic $\mathrm{O}/\mathrm{C}$ ratios is generally weaker than on vapor pressure (Koop et al., 2011), strong correlations between T_g and $\mathrm{O}/\mathrm{C}$ ratios have been observed for oxidation products formed from specific precursors, such as α-pinene (Dette et al., 2014), n-heptadecane, and naphthalene (Saukko et al., 2012). Therefore, the uncertainties in volatility and $\mathrm{O}/\mathrm{C}$ ratios may impact T_g and viscosity predictions when using the parameterizations that relate T_g and viscosity to volatility and $\mathrm{O}/\mathrm{C}$ ratios (Li et al., 2020; Zhang et al., 2019).

Here we calculated the T_g of individual OA surrogate species i (T_g,i) using Eqs. (1), (2) (Li et al., 2020) or (3) (Zhang et al., 2019):

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{T}_{\mathrm{g},i}& =\mathrm{288.70}-\mathrm{15.33}\times {\log }_{\mathrm{10}}\left(C_i^{\mathrm{0}}\right)\\ & -\mathrm{0.33}\times {\left[{\log }_{\mathrm{10}}\left(C_i^{\mathrm{0}}\right)\right]}^{\mathrm{2}}\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{T}_{\mathrm{g},i}& =\mathrm{289.10}-\mathrm{16.5}\times {\log }_{\mathrm{10}}\left(C_i^{\mathrm{0}}\right)-\mathrm{0.29}\\ & \times {\left[{\log }_{\mathrm{10}}\left(C_i^{\mathrm{0}}\right)\right]}^{\mathrm{2}}+\mathrm{3.23}\times {\log }_{\mathrm{10}}\left(C_i^{\mathrm{0}}\right)\times \left(\mathrm{O}/\mathrm{C}\right)\end{array}\text{(3)}& {\displaystyle }{T}_{\mathrm{g},i}=\mathrm{480.07}-{\displaystyle \frac{\mathrm{54}\mathrm{395}}{{\left({\log }_{\mathrm{10}}\left(\frac{RT}{M_i}{C}_i^{\mathrm{0}}\right)-\mathrm{7.7929}\right)}^{\mathrm{2}}+\mathrm{116.49}}}\end{array}$$

Where Eq. (1) relies solely on volatility ($C_i^{\mathrm{0}}$ represents the saturation concentration at 298 K, and is equal to $C_i^{\ast }$ assuming ideal thermodynamic mixing), and was developed for coupling into the 1-D VBS framework. Equation (2) incorporates $\mathrm{O}/\mathrm{C}$ as an additional factor, which is used in 2D-VBS. Both equations yield similar predictions, particularly for compounds with low $\mathrm{O}/\mathrm{C}$ ratios (Li et al., 2020). Equation (3) is a semi-empirical formula derived from the measurements of isoprene SOA components, which relates T_g to volatility and includes molar mass M_i (Zhang et al., 2019). The overall T_g of OA mixtures under dry conditions (T_g,org) can be calculated using the Gordon-Taylor equation (see Sect. S2 in the Supplement).

OA viscosity depends on aerosol water content, since water significantly influences the phase state of aerosols (Koop et al., 2011). Under humid conditions, the T_g of organic-water mixtures (T_g,ωorg) is calculated based on the mass fraction of organics (ω_org) in the mixtures, along with the T_g values of pure water (136 K) and dry OA (using Eqs. 1, 2 or 3). Notably, ω_org varies with ambient relative humidity (RH) and the effective hygroscopicity parameter of OA (κ_org). Further details can be found in previous studies (DeRieux et al., 2018; Shiraiwa et al., 2017) and Sect. S2 in the Supplement. The temperature-dependent viscosity (η) can be estimated using the modified Vogel-Tammann-Fulcher (VTF) equation, as given in Eq. (4) (Angell, 1991), when the ambient temperature T is at or above T_g,ωorg; otherwise η is fixed at 10¹² Pa s.

$$\begin{array}{}\text{(4)}& \mathit{\eta }={\mathit{\eta }}_{\mathrm{\infty }}{e}^{\frac{T_{\mathrm{0}}D}{T-T_{\mathrm{0}}}}\end{array}$$

where η_∞ and D denote the viscosity at infinite temperature and the fragility parameter, with the values of 10⁻⁵ Pa s and 10, respectively (DeRieux et al., 2018). T₀ represents the Vogel temperature, and can be estimated using Eq. (5) (Angell, 1991; DeRieux et al., 2018):

$$\begin{array}{}\text{(5)}& T_{\mathrm{0}}={\displaystyle \frac{\mathrm{39.17}{T}_{\mathrm{g},\mathit{\omega }\mathrm{org}}}{D+\mathrm{39.17}}}\end{array}$$

2.5 Measurement Data

Two sets of observational data from field campaigns were used to evaluate the simulations. One campaign was conducted in the spring (17 March–21 April) of 2018 in DY, at a site located in the Yellow River Delta National Nature Reserve and strongly influenced by urban outflows from the North China Plain (NCP) (Feng et al., 2023). The observational data include the mass concentration and elemental ratios (e.g., $\mathrm{O}/\mathrm{C}$, $\mathrm{H}/\mathrm{C}$) of OA, measured using a high-resolution time-of-flight aerosol mass spectrometer (HR-AMS). The volatility distribution of OA (expressed as VBS) was estimated using an empirical method based on data from a thermodenuder (TD) combined with AMS. The sources of OA components were characterized using the positive matrix factorization (PMF) method, which resolved OA factors, including hydrocarbon-like OA (HOA), biomass-burning OA (BBOA), transported OOA and background OOA. The other campaign, led by the same team, was conducted in the autumn (29 September–21 November) of 2018 in urban GZ surrounded by transportation and residential districts, with significant impacts from local emissions (Chen et al., 2021). The resolved OA factors in GZ identified three POA groups (cooking OA (COA), HOA, and nitrogen-containing OA (NOA)) along with two types of SOA (low-volatility oxygenated OA, LV-OOA, and semi-volatile oxygenated OA, SV-OOA).

3 Results

3.1 Evaluation of mass concentration

Figure 1 compared the observed and simulated diurnal variations of OA and its components, POA and SOA, in DY and GZ, respectively. The simulation with updated $\mathrm{O}/\mathrm{C}$ ratios of POA emissions in the 2D-VBS framework (i.e., 2D-VBS_A) did not alter the emission magnitudes or volatility, and therefore had a minimal impact on the predicted POA mass concentration or volatility distribution. As a result, this case is not discussed in Sect. 3.1 and 3.2. The observations indicated that SOA dominated OA at both sites, contributing 72 % and 64 % at DY and GZ, respectively. The model generally underestimated the SOA contribution in DY (except for the 2D-VBS case), while slightly overestimating it in GZ (Table S9). In general, the 2D-VBS case predicted higher contributions of SOA to total OA (> 80 %) than in the 1D-VBS simulations.

Figure 1 Simulated diurnal variations of OA, POA, and SOA with different SOA schemes compared to the observations in DY (a–c) and GZ (d–f). The light grey shading indicates the standard deviations (±1σ) from the mean concentration in the observations.

The base 1D-VBS significantly underestimated SOA production, with NMBs of −72 % in DY and −71 % in GZ (Table S10). Adding L/S/IVOC emissions and their contributions to SOA in the 1D-VBS_E case resulted in approximately a two-fold enhancement in SOA mass concentration. This finding was consistent with earlier studies showing that 64 % to 100 % of observed SOA in GZ can be explained when both conventional precursors and S/IVOCs are considered (Hu et al., 2022), and that simulated S/IVOC contributions to SOA exceeded 50 % across most of China (Li et al., 2021a; Miao et al., 2021). The 1D-VBS_EY case substantially increased SOA production from aromatics and PAHs with updated yields (Fig. S6), and the diurnal variation was similar to that in the 2D-VBS simulation in DY. In GZ, the notable daytime differences in SOA between the 2D-VBS and 1D-VBS_EY cases were driven by variations in ASOA representations (Fig. S7). Mass increases in ASOA as a result of gas-phase aging were likely overestimated with the 4 km grid spacing in the 1D-VBS_EY simulation (Bilsback et al., 2023). While both simulations significantly improved SOA predictions, they still underestimated SOA concentrations compared to observations by 11 %–14 % (2D-VBS) and 30 %–33 % (1D-VBS_EY), respectively. Diagnosing the causes of these underestimations remains challenging due to the lack of chemically resolved SOA measurements. Although POA was underestimated (as discussed below), the sensitivity simulations indicated that its impact on SOA was limited (Fig. S8). At the GZ site, for example, increasing POA emissions by 70 % in the 1D-VBS_EY simulation to match the observed POA concentrations led to only 4 % increase in SOA, suggesting that missing or incomplete SOA formation pathways are likely the primary drivers of the SOA low biases. We infer that the nighttime low biases in both cases were likely attributed to insufficient formation of organic nitrates in the presence of NO₃ radicals. This might result from missing anthropogenic terpene and phenolic emissions from sources such as biomass burning and volatile chemical products (Wang et al., 2022; Coggon et al., 2021; Xie et al., 2025; Liu et al., 2024). The underestimated nocturnal SOA formation could also be attributed to the underprediction of aqueous-phase formation pathways, which are enhanced under high RH conditions at night (Wang et al., 2019; Gu et al., 2023).

Simulated POA also exhibited underestimation at both locations. In the comparison, BBOA and COA were excluded from the observations, due to the complexity and considerable uncertainties associated with cooking and biomass burning (particularly open burning) emission estimation (Li et al., 2023; Zhou et al., 2017). Therefore, the POA underestimation was most likely attributable to uncertainties in mobile emissions. Observed POA concentrations at both sites showed peaks during traffic rush hours, which were not well captured by the simulations. At the DY site, the resolved HOA in observations was partially aged, as noted in Feng et al. (2023) and included some SOA transported over long distances, which may also explain the POA underestimation. At the GZ site, which is more influenced by local emissions, the simulated peak concentration exhibited an earlier shift compared to the observed peak, suggesting that the diurnal patterns of POA emissions may differ from those specified in the model. Additionally, the POA volatility distributions (see Methods) and aging schemes in the 1D-VBS and 2D-VBS simulations were different: (1) The 1D-VBS involved heterogeneous aging of non-volatile POA. (2) In 2D-VBS, POA aging led to its complete transformation into SOA, whereas in the 1D-VBS, some oxidation products remained as POA. These differences led to distinct POA levels between the two cases.

Overall, total mass concentrations of OA were underestimated. The 1D-VBS_EY case, which had the smallest negative biases among the 1D-VBS simulations, underestimated OA by approximately 27 % in DY and 34 % in GZ. These discrepancies were dominated by SOA underestimation at both sites (Fig. S9). The 2D-VBS simulation demonstrated better performance in OA predictions than 1D-VBS_EY, with NMBs of −24 % in both DY and GZ, primarily due to higher predicted SOA mass concentrations.

3.2 Volatility simulations

The mass distributions of OA, POA, and SOA across different volatility bins (i.e., volatility distribution) were examined. This study primarily focused on compounds with C^∗ equal or below 10⁻¹ µg m⁻³ (log$C^{\ast }\le -\mathrm{1}$, referred to as low-volatility OA hereafter), and on SVOCs (i.e., $-\mathrm{1}<\log C^{\ast }\le \mathrm{2}$), which exist in both gas and condensed phases. Despite variations in OA mass concentrations across different cases, CMAQ generally reproduced the observed volatility distributions of OA, SOA, and POA in all the simulations (Fig. 2). Specifically, the observations indicated that more than 70 % of OA was low-volatility OA, consistent across the simulations. The base 1D-VBS significantly underestimated low-volatility OA concentrations, leading to underestimation in total OA mass (Fig. S10). The 1D-VBS_E simulation, which included L/S/IVOC emissions, increased low-volatility SOA formation by 1.1 µg m⁻³ in DY and 1.5 µg m⁻³ in GZ, although the contributions of low-volatility SOA to total SOA decreased. It also replaced the default volatility distribution of POA in CMAQv5.3 with source-specific gas/particle partitioning, resulting in a significant amount of non-volatile POA (Fig. S11). Collectively, these adjustments led to slightly higher fractions of low-volatility OA with the inclusion of L/S/IVOC emissions, as proposed by Xu et al. (2019), albeit for a different reason (mainly due to changes in POA rather than SOA volatility distributions). The updated SOA yields for aromatics and PAHs had minor effects on POA but increased SOA mass concentrations across three of the four volatility bins in Fig. S10, slightly changing the SOA volatility distributions.

Figure 2 Simulated volatility distributions of OA, POA, and SOA across volatility bins with C^∗ up to 10² µg m⁻³, compared to the observations in DY (a–c) and GZ (d–f). The compounds with C^∗ equal or below 10⁻¹ µg m⁻³ were grouped into a single volatility bin, and referred to as low-volatility OA in this study.

The simulated OA volatility distributions, particularly for POA, differed between the 2D-VBS and 1D-VBS_EY simulations. The observations presented a higher fraction of low-volatility POA in GZ than in DY. The spatial variation in POA volatility was captured by the 1D-VBS_EY case. However, both 1D-VBS_EY and 2D-VBS overestimated the contributions of low-volatility POA. In 1D-VBS_EY, non-volatile POA was substantial in mass, and increased with aging (Fig. 3); semi-volatile POA underwent aging, with a portion remaining as POA (rather than fully converting into SOA as in 2D-VBS) but becoming more volatile. This ultimately led to a higher low-volatility POA mass and a lower contribution compared to 2D-VBS. The SOA volatility distributions in the 1D-VBS_EY and 2D-VBS cases both agreed well with the observations, which were similar at the two locations. However, the 2D-VBS simulation predicted a lower fraction of low-volatility BSOA in GZ compared to the 1D-VBS simulations (over 73 % vs. less than 86 %), primarily due to differences in the treatment of monoterpene-derived SOA through ${\mathrm{O}}_{\mathrm{3}}/\mathrm{OH}$ oxidation pathways. Additionally, the contribution of low-volatility ASOA was also lower in the 2D-VBS simulation (over 62 % vs. over 69 %; see Fig. S12).

Figure 3 Simulated mass concentrations of OA, POA, and SOA across volatility bins with C^* ranging from 10⁻¹⁰ to 10² µg m⁻³, compared to observations in DY (a–c) and GZ (d–f).

It is important to note that LVOCs ($-\mathrm{4}<\log C^{\ast }$ (µg m⁻³) $\le -\mathrm{1}$) and ELVOCs (extremely low-volatility VOCs with logC^∗ (µg m⁻³) $\le -\mathrm{4}$) were not explicitly resolved in this study (Fig. 3). In all simulations, SOA formation was dominated by species within C^∗ bins ranging from 10⁻² to 10⁴ µg m⁻³. In addition, certain assumed “non-volatile” products (IEPOX/MAE-derived SOA, oligomers, etc.) were arbitrarily assigned to the volatility bin of $C^{\ast }={\mathrm{10}}^{-\mathrm{10}}$ µg m⁻³ (see Methods). The abundance of mass in the $C^{\ast }={\mathrm{10}}^{-\mathrm{2}}$ µg m⁻³ bin was particularly high for the 2D-VBS and largely due to aging. The absence of predicted SOA between the C^∗ bins at 10⁻¹⁰ µg m⁻³ and $C^{\ast }={\mathrm{10}}^{-\mathrm{2}}$ µg m⁻³ was in contrary to the observed more uniform distribution across C^∗ bins spanning from 10⁻⁹ to 10⁻¹ µg m⁻³. The 1D-VBS_EY case included formation of mass for the volatility bin at $C^{\ast }={\mathrm{10}}^{-\mathrm{6}}$ µg m⁻³ for HOMs formed from aromatics and PAHs under low-NO_x conditions, but the predicted mass concentration was negligible and biased low relative to the observations. Major ELVOC formation pathways, including autoxidation and bimolecular peroxy radical reactions for monoterpenes, have been incorporated into the Community Regional Atmospheric Chemistry Multiphase Mechanism (CRACMM) (Pye et al., 2023). However, the volatility distribution of the resulting products remained insufficiently resolved. Furthermore, anthropogenic ELVOC formation is not yet well understood (Shrivastava et al., 2024; Yin et al., 2024). As ambient observations consistently reported significant amounts of LVOCs and ELVOCs (Chen et al., 2024; Huang et al., 2024), future work should refine and expand the representation of autoxidation and other chemical processes contributing to LVOC and ELVOC formation from both anthropogenic and biogenic precursors. These improvements will enhance model accuracy in predicting SOA volatility and new particle formation.

3.3 Oxygen-to-carbon ratio simulations

Both observed and simulated OA showed higher $\mathrm{O}/\mathrm{C}$ ratios during the day and lower ratios at night (Fig. S13), aligning with the temporal patterns of SOA mass concentrations. The higher $\mathrm{O}/\mathrm{C}$ ratios in the afternoon were linked to elevated oxidant levels, which facilitated SOA formation (typically characterized by higher $\mathrm{O}/\mathrm{C}$ ratios compared to POA, see Table 2) and the photochemical aging of OA. In contrast, during the night and early morning, the increased contributions of fresh POA led to lower $\mathrm{O}/\mathrm{C}$ ratios.

Table 2 Simulated and observed $\mathrm{O}/\mathrm{C}$ ratios of OA, POA and SOA in DY and GZ.

Most simulations underestimated the $\mathrm{O}/\mathrm{C}$ ratios of OA, with a more pronounced underestimation at the DY site, where both POA and SOA $\mathrm{O}/\mathrm{C}$ ratios were lower than observed (Table 2). In DY, OA underwent prolonged aging during transport, resulting in a higher $\mathrm{O}/\mathrm{C}$ ratio (0.84) than the national average (0.3–0.65) (Feng et al., 2023). Although the 1D-VBS_EY case predicted the highest $\mathrm{O}/\mathrm{C}$ ratio of OA (due to potentially overestimated multi-generational oxidation as discussed in Sect. 3.1), it was still lower than the observed value (0.61 vs. 0.83). This discrepancy suggests that accurately simulating the evolution of OA $\mathrm{O}/\mathrm{C}$ ratios with aging remains a significant challenge for current CTMs (despite satisfactory mass simulations). The observed POA contained aged HOA, which, although essentially SOA, could not be separated from HOA in DY (Feng et al., 2023), partly explaining the underestimation of POA $\mathrm{O}/\mathrm{C}$ ratios. In the 1D-VBS_E case, the inclusion of L/S/IVOC emissions increased the $\mathrm{O}/\mathrm{C}$ ratio of POA but decreased that of SOA, consistent with the impacts on volatility, i.e., increased (decreased) the contribution of low-volatility POA (SOA) (see Figs. 2 and S11). The 1D-VBS_EY case promoted SOA formation at $C^{\ast }=\mathrm{0.1}$ µg m⁻³ through the newly added volatility bin for ASOA, which was highly oxygenated ($\mathrm{O}/\mathrm{C}$ = 2.53, see Table S7), thereby increasing the $\mathrm{O}/\mathrm{C}$ ratios of SOA. However, the impacts of HOMs were negligible due to their minimal contribution to mass.

At the GZ site, the overestimation of POA $\mathrm{O}/\mathrm{C}$ ratios in certain cases (i.e., 1D-VBS_E and 1D-VBS_EY) partially offset the underestimation of SOA $\mathrm{O}/\mathrm{C}$ ratios. As a result, the $\mathrm{O}/\mathrm{C}$ ratio of total OA in 1D-VBS_E was closest to the observed value (both were 0.59), while the 1D-VBS_EY case overpredicted the $\mathrm{O}/\mathrm{C}$ ratio (0.66 vs. 0.59). The more oxygenated POA in both simulations than observed could be attributed to overestimations in non-volatile POA mass concentrations. Additionally, COA, with an observed $\mathrm{O}/\mathrm{C}$ ratio of 0.19 in GZ, were likely underestimated in the simulations given uncertainties in cooking emissions. If COA were well represented, it could slightly lower the simulated POA $\mathrm{O}/\mathrm{C}$ ratios and reduce the positive biases in the 1D-VBS_E and 1D-VBS_EY simulations in GZ. It was also found that wind speeds were overestimated at GZ (Fig. S14), which could lead to excessive dilution and consequently reduce the simulated concentrations of locally emitted POA and freshly formed SOA (both of which typically exhibit lower $\mathrm{O}/\mathrm{C}$ ratios than aged POA or SOA). This effect would be expected to increase the simulated $\mathrm{O}/\mathrm{C}$. However, the low biases in $\mathrm{O}/\mathrm{C}$ at GZ indicate that meteorological biases are unlikely to be the dominant factor, and that limitations in the SOA representations may play a more important role.

The $\mathrm{O}/\mathrm{C}$ ratios predicted by the 2D-VBS case were lower than those in other simulations, primarily due to the underestimation in SOA $\mathrm{O}/\mathrm{C}$ ratios. In particular, the distribution of ASOA $\mathrm{O}/\mathrm{C}$ ratios differed significantly from those in the 1D-VBS simulations (see Fig. S15), highlighting the constraint of ASOA $\mathrm{O}/\mathrm{C}$ ratios in future work. Simulated POA $\mathrm{O}/\mathrm{C}$ ratios were mainly influenced by how emissions were specified in the model. As the default settings in 2D-VBS underestimated POA $\mathrm{O}/\mathrm{C}$ ratios, an additional simulation, 2D-VBS_A, was conducted with updated $\mathrm{O}/\mathrm{C}$ ratios for emissions, constrained by those of POA factors in observations from prior studies (Table S6 and Fig. S16). For instance, in the default configuration, mobile sources were assumed to emit primarily hydrocarbons and low-oxygenated compounds, with over 90 % of emissions distributed in the $\mathrm{O}/\mathrm{C}$ bins of 0 and 0.1. The optimized parameters in the 2D-VBS_A case reflected a higher degree of oxygenation, with more than 60 % of emissions allocated to the $\mathrm{O}/\mathrm{C}$ bins of 0.2 and 0.4. As a result, the POA $\mathrm{O}/\mathrm{C}$ increased from 0.26 (DY) and 0.20 (GZ) in 2D-VBS to 0.29 and 0.25 in 2D-VBS_A, respectively, and showed good agreement with the observation in GZ. However, the persistent low bias in POA $\mathrm{O}/\mathrm{C}$ ratios in DY was likely attributable to uncertainties in emissions from sources with higher $\mathrm{O}/\mathrm{C}$, such as biomass burning (with an observed $\mathrm{O}/\mathrm{C}$ ratio of 0.37 for BBOA) in the simulations, as well as potential overestimation of observed POA $\mathrm{O}/\mathrm{C}$ ratios due to the influence of aged HOA.

Apart from underestimations in OA $\mathrm{O}/\mathrm{C}$ ratios in most cases, the simulations struggled to capture the spatiotemporal variability on an hourly basis (Fig. 4). The observations revealed significant differences in the $\mathrm{O}/\mathrm{C}$ ratios between DY and GZ, reflecting distinct aging processes associated with site characteristics (i.e., regional transport at DY versus dominant local emissions at GZ), rather than differences in oxidant levels, as comparable O₃ concentrations were observed at both sites. As discussed above, current SOA parameterizations inadequately represent chemical aging processes and thereby fail to reproduce the observed spatial contrasts in oxidation state. The observations also exhibited considerable temporal variability, particularly at the GZ site with abundant fresh emissions. Specifically, the interquartile range (IQR) of hourly $\mathrm{O}/\mathrm{C}$ ratios was 0.77–0.89 in DY and 0.52–0.68 in GZ in the observations. However, the simulated $\mathrm{O}/\mathrm{C}$ ratios of OA showed less variability than observed in GZ, with narrower IQRs and Standard Deviation Ratios (SDR) lower than 1.0.

Figure 4 Density distribution of simulated versus observed hourly $\mathrm{O}/\mathrm{C}$ ratios of OA in DY (a–e) and GZ (f–k). The edges of the grey shading represent the 25th and 75th percentiles in both the observations and simulations, with the shading indicating the interquartile range (IQR). The Index of Agreement (IOA, reflecting both bias and variability) and Standard Deviation Ratio (${\mathit{\sigma }}_{\mathrm{sim}}/{\mathit{\sigma }}_{\mathrm{obs}}$, comparing the magnitude of variability in simulations and observations) are calculated.

The responses of daily OA $\mathrm{O}/\mathrm{C}$ ratios to ambient O₃ concentrations were also examined (Fig. 5). All the simulations demonstrated an increased in $\mathrm{O}/\mathrm{C}$ ratios as O₃ levels rose (indicating a higher atmospheric oxidative capacity), with the rates in the range of 0.0023–0.0045 per ppb O₃. The 2D-VBS simulations exhibited a steeper slope, suggesting stronger sensitivity to O₃ levels. This implied that oxidant abundance in the gas phase was an important driver of the increase in $\mathrm{O}/\mathrm{C}$ ratios in the simulations. However, the observations did not show a clear correlation between $\mathrm{O}/\mathrm{C}$ ratios and O₃. In DY, this could be due to the decoupling of $\mathrm{O}/\mathrm{C}$ ratios from local O₃ in aged air masses. Other factors, such as fresh emissions in GZ (see Fig. S17) or aqueous formation pathways, may also contribute to the observed variations in $\mathrm{O}/\mathrm{C}$ ratios.

Figure 5 Simulated changes in OA $\mathrm{O}/\mathrm{C}$ ratios with daily mean ozone (O₃) concentrations, compared to observations in DY (a) and GZ (b). The paired data of $\mathrm{O}/\mathrm{C}$ ratios and O₃ were grouped into three bins with equal intervals of O₃ levels. Dashed lines represent fitted trends, and the values of N denote the sample size within each bin.

3.4 Implications for T_g and viscosity predictions

Several widely used CTMs, including WRF-Chem and GEOS-Chem, have been utilized to predict OA viscosity and phase state based on volatility and $\mathrm{O}/\mathrm{C}$ (Zhang et al., 2024; Luu et al., 2025). In this study, three parameterizations of T_g were employed to investigate how model capability in representing volatility and $\mathrm{O}/\mathrm{C}$ influences the predictions of T_g (and the resulting viscosity using the VTF equation, see Methods): (1) one incorporating both $\mathrm{O}/\mathrm{C}$ and volatility (Eq. 2), (2) one accounting for both volatility and molecular mass (Eq. 3), and (3) one based solely on volatility (Eq. 1). Given that Eq. (3) was developed from laboratory experiments focused on isoprene SOA, its applicability to POA and ASOA warrants further evaluation (which is beyond the scope of this study). Nevertheless, it was included here as a comparative parameterization to evaluate impacts of volatility representations. The simulation at the DY site was selected for detailed analysis as a representative case.

In DY, the mean T_g,org values for POA and SOA showed differences across parameterizations, ranging from 15.5 to 24.4 K for POA and 19.3 to 23.0 K for SOA (Fig. 6). The parameterization based on Eq. (1) generally predicted the highest T_g,org values, closely aligning with those from Eq. (2) as has been reported in Li et al. (2020). This suggests that simulated volatility distributions exerted a stronger influence on T_g,org predictions than $\mathrm{O}/\mathrm{C}$ ratios, as further evidenced by the T_g,org values of individual SOA surrogates (Fig. S18). Equation (3) predicted the lowest T_g,org values, including for nearly all SOA components except a few isoprene SOA surrogates. For isoprene SOA components, T_g estimates remained constant using either Eqs. (1) or (3), due to the use of a fixed C^* of 10⁻¹⁰ µg m⁻³ in CMAQ (see Methods), which resulted in unrealistically high T_g values, e.g., ∼ 400 K with Eq. (1). However, 2-methyltetrol (AIETET), despite its unexpectedly high ambient concentrations, might be semi-volatile, with predicted C^* of 10² µg m⁻³, while 2-methyltetrol sulfate (AIEOS) has a relatively low volatility, with C^* of 10⁻¹ µg m⁻³ (Budisulistiorini et al., 2017). These values resulted in much lowerT_g estimates of approximately 230 and 276 K, respectively (Zhang et al., 2019). Therefore, the C^* values for these compounds should be revisited when applied in T_g predictions.

Figure 6 Predicted T_g,org of (a) POA and (b) SOA using three parameterizations across different simulation cases. The solid line represents T_g,org derived from observed volatility and $\mathrm{O}/\mathrm{C}$ ratios. Symbols denote T_g,org estimated from: volatility and $\mathrm{O}/\mathrm{C}$ ratio using Eq. (2) (circles), volatility and molecular weight using Eq. (3) (triangles), and volatility only using Eq. (1) (diamonds). Error bars indicate one standard deviation.

The discrepancies in POA T_g,org derived from observed versus simulated volatility and $\mathrm{O}/\mathrm{C}$ ratios using Eq. (2) ranged from 12.4 to 33.7 K across the five simulations. Compared to observation-based T_g,org, the 2D-VBS, 2D-VBS_A, and 1D-VBS simulations underestimated POA T_g,org, whereas the 1D-VBS_E and 1D-VBS_EY simulations overpredicted it, primarily due to a substantial contribution from non-volatile POA (Fig. S11). Although the 2D-VBS simulation exhibited the highest fraction of low-volatility POA (Fig. 2), these species were mainly allocated to the volatility bin of $C^{\ast }={\mathrm{10}}^{-\mathrm{2}}$ µg m⁻³, which led to lower T_g,org values than those inferred from observations. The biases in SOA T_g,org were smaller, ranging from 2.5 to 17.2 K. Most simulations exhibited a slight underestimations of SOA T_g,org relative to observation-derived values. Conversely, the 1D-VBS case overestimated SOA T_g,org, with an excessive proportion of SOA at $C^{\ast }={\mathrm{10}}^{-\mathrm{10}}$ µg m⁻³ (Fig. S11). Since Eqs. (1) and (2) were developed for use within the 1D-VBS and 2D-VBS frameworks, respectively, the evaluation presented here suggests that the differences in T_g,org estimates were driven by the volatility representations in the 1D-VBS versus 2D-VBS, rather than by the specific T_g parameterizations. Notably, the 2D-VBS implementation in CMAQv5.4 generally produced lower SOA T_g,org values.

The values of OA viscosity η, which determine their phase state, were calculated using the VTF equation, with T_g parameterized as a function of volatility and $\mathrm{O}/\mathrm{C}$ ratios (Eq. 2), and incorporating additional variables including κ_org, RH, and T (see Methods). Figure 7a and e compared viscosity estimated from simulations and observations, with all input variables (i.e., T_g,org, κ_org, RH and T) obtained consistently from either model output or observations. In the simulations, κ_org for each OA species was parameterized as a function of $\mathrm{OM}/\mathrm{OC}$ (Eq. S7), assuming a constant density of 1.4 g cm⁻³. Bulk κ_org for SOA and POA was calculated as a mass-weighted average across species. For observations, κ_org was inferred from f₄₄ (i.e., the fraction of $m/z$ 44 signal in total organic signals) following Feng et al. (2023). The simulations generally reproduced the diurnal variability in OA viscosity, predicting higher η values during the daytime than at night, in line with the observed diurnal patterns. Field measurements indicated that POA was predominantly semi-solid (${\mathrm{10}}^{\mathrm{2}}\le \mathit{\eta }\le {\mathrm{10}}^{\mathrm{12}}$ Pa s) with a transition to the solid phase (η>10¹² Pa s) between 12:00 and 17:00 LT. However, the model overestimated viscosity. Consequently, the 2D-VBS and 2D-VBS_A simulations predicted prolonged solid phase for POA, while 1D-VBS_E and 1D-VBS_EY predicted POA to remain solid throughout the day. For SOA, observations showed that it was semi-solid, while the simulations reproduced this phase state during nighttime hours only.

Figure 7 Diurnal variations in predicted OA viscosity using Eq. (2), based on observations and simulations. Black points indicate observation-based values, while colored lines represent predictions using combinations of T_g,org, κ_org, RH, and T from observations or simulations as defined in Table 3. Red points indicate viscosity estimated using observation-based T_g,org and κ_org values, along with simulated RH and T. The OA phase state is determined based on viscosity: liquid (η<10² Pa s), semi-solid (${\mathrm{10}}^{\mathrm{2}}\le \mathit{\eta }\le {\mathrm{10}}^{\mathrm{12}}$ Pa s), and solid (η>10¹² Pa s) (Reid et al., 2018). Viscosity values in the figure are capped at 10¹² Pa s, beyond which values are not physically meaningful for solid-phase OA.

To identify the dominant source of overestimation in predicted viscosity, three sensitivity experiments were conducted (Table 3), each isolating the influence of a specific input variable. The results from sensitivity experiment A (Fig. 7b and f) revealed that the evident model-observation discrepancies in POA and SOA viscosity were positively correlated with the biases in T_g,org. For instance, in the 1D-VBS case, POA viscosity was significantly underestimated (by up to five orders of magnitude) due to a ∼ 40 K negative bias in T_g,org (Fig. S19a), while an overestimation of T_g,org for SOA (∼ 20 K, Fig. S19d) led to viscosity exceeding observation-based values by up to two orders of magnitude. In contrast, simulations such as 1D-VBS_E and 1D-VBS_EY, which better captured SOA T_g,org, yielded viscosity values in good agreement with observations. Thus, by comparing Fig. 7e and f, we conclude that the considerable overestimation of SOA viscosity in Fig. 7e, except in 1D-VBS, cannot be attributed to T_g,org bias. O'Brien et al. (2021) also reported a moderate sensitivity of predicted OA viscosity to T_g,org, e.g., under conditions of RH = 60 %, a reduction in T_g,org by ∼ 200 K resulted in a decrease in predicted viscosity from 10¹² to 10 Pa s (with κ_org of 0.14).

Table 3 Combinations of T_g,org, κ_org, RH and T used in sensitivity cases for OA viscosity estimation.

On the other hand, the OA hygroscopicity parameter (κ_org) was an important parameter in viscosity estimation. Sensitivity experiment B predicted viscosity comparable to those in Fig. 7a and e, indicating that the discrepancies between simulation- and observation-based viscosity values was due to the variation in κ_org. Specifically, the κ_org values used in the simulations underestimated the hygroscopicity of both POA and SOA relative to the values derived from $m/z$ 44 (Fig. S20), e.g. 0.16 versus 0.048–0.069 (1D-VBS)$/$0.064–0.068 (2D-VBS) for POA, and 0.35 verus 0.013–0.014(1D-VBS)$/$0.12 (2D-VBS) for SOA. Earlier sensitivity simulations showed that a perturbation of ±0.05 in κ_org resulted in changes of T_g,ωorg by ∼ 5 %–15 % over high RH areas, supporting a critical role of κ_org in viscosity estimation (Shiraiwa et al., 2017). Consistently, the underestimation in κ_org led to an overestimation in simulated viscosity in this study (given the inverse relationship between viscosity and hygroscopicity), and even offset the underestimations of viscosity that were caused by low biases in T_g,org (see Sensitivity experiment A). In CMAQ, κ_org was parameterized as a function of $\mathrm{OM}/\mathrm{OC}$. In light of the notable underestimation of $\mathrm{O}/\mathrm{C}$ ratios at the DY site (Table 2), this could be a key driver of the bias in viscosity predictions. While the impact of relative humidity on viscosity has also been emphasized in other studies (Maclean et al., 2021; Rasool et al., 2021; Li et al., 2021b), its effect is minimal in this study due to the strong agreement between simulated and observed meteorological variables (Fig. S21). This was further supported by the similarity in viscosity estimates derived using identical κ_org and T_g,org but different RH and T inputs (Fig. 7d and h).

4 Conclusions and Discussion

This study evaluated the CMAQ model's performance in predicting OA mass concentrations, volatility distributions, and $\mathrm{O}/\mathrm{C}$ ratios, and examined their implications for T_g and viscosity estimates, using observations at two representative sites in Eastern China, i.e., DY, influenced by aged air masses, and GZ, more impacted by local emissions. The major findings include:

1. The base 1D-VBS simulation with AERO7i underestimated OA by ∼ 70 %. However, the optimized 1D-VBS simulation, which incorporated additional L/S/IVOC emissions and enhanced SOA production from aromatics and PAHs (accounting for vapor wall losses and autooxidation), significantly improved OA predictions with NMBs of −27 % (DY) and −34 % (GZ). The 2D-VBS simulation exhibited similar performance with NMBs of −24 %. The OA biases were mainly driven by SOA underestimation, with the optimized 1D-VBS underestimating SOA by ∼ 32 % and 2D-VBS underestimating it by ∼ 13 %.

2. All the simulations accurately captured the observed fraction of low-volatility products in OA (> 70 %), although the detailed distributions within these volatility bins (10⁻⁹ µg m${}^{-\mathrm{3}}\le C^{\ast }\le \mathrm{0.1}$ µg m⁻³) were not well represented. Simulated SOA was abundant in volatility bins of $C^{\ast }={\mathrm{10}}^{-\mathrm{2}}$ µg m⁻³ or 10⁻¹⁰ µg m⁻³. The volatility distributions of POA showed notable variations across the simulations, influenced by the gas/particle partitioning of L/SVOC emissions specified in the model.

3. OA $\mathrm{O}/\mathrm{C}$ ratios were generally underestimated across all the simulations, particularly in DY, where aged air masses prevailed. In the 1D-VBS simulations, the inclusion of L/S/IVOC emissions increased $\mathrm{O}/\mathrm{C}$ ratios with more non-volatile POA, while updated SOA yields added SOA at $C^{\ast }=\mathrm{0.01}$ µg m⁻³ , contributing to higher SOA $\mathrm{O}/\mathrm{C}$ ratios. The 2D-VBS simulations presented lower $\mathrm{O}/\mathrm{C}$ ratios than 1D-VBS overall, whereas constraining POA with observations improved its $\mathrm{O}/\mathrm{C}$ representation.

The simulated volatility distribution strongly influenced T_g estimates. In particular, the fraction of POA and SOA in the lowest volatility bin ($C^{\ast }={\mathrm{10}}^{-\mathrm{10}}$ µg m⁻³) played an important role. The assumption of $C^{\ast }={\mathrm{10}}^{-\mathrm{10}}$ µg m⁻³ resulted in unrealistically high T_g values for some OA components (e.g., isoprene SOA). Although the high (low) biases in T_g led to overestimation (underestimation) of viscosity in the case of DY, uncertainties in OA hygroscopicity parameter, parameterized as a function of $\mathrm{O}/\mathrm{C}$ ratios in CMAQ, were the dominant source of model-observation gaps in viscosity.

These findings highlight the limitations of the CMAQ model in predicting the volatility distribution of OA and the associated $\mathrm{O}/\mathrm{C}$ ratios. While the allocation of OA in low-volatility bins may have a limited impact on simulated mass concentrations under atmospherically relevant conditions, it is closely linked to key physicochemical properties, such as T_g, viscosity, and phase state. These properties, which have received insufficient attention in previous studies, play a crucial role in the kinetics of multiphase SOA formation. Additionally, estimates of OA hygroscopicity and light absorptivity – both of which influence climate effects – may be subject to uncertainties when parameterized based on $\mathrm{O}/\mathrm{C}$ ratios. In some cases, inaccurate treatment of these properties could also influence OA mass concentrations. For example, OA hygroscopicity affects water uptake, which can, in turn, enhance SOA mass through interactions between water and organics. Therefore, we summarize several recommendations based on this study for future improvements in OA modeling:

Given the substantial contribution of L/S/IVOCs to SOA, emission inventories should be further refined not only in terms of total magnitudes but also volatility-resolved distributions, and constrained using ambient measurements. In addition, primary emissions from sources such as cooking, open biomass burning, and mobile sources require improved representation. Better treatment of nighttime SOA formation pathways, particularly NO₃ oxidation and aqueous-phase chemistry, is also needed to reduce SOA mass underestimation.

While updated SOA yields can partially improve model performance, explicitly accounting for autoxidation processes and the formation of HOMs (including both biogenic and anthropogenic origins) would provide a more physically based description of $\mathrm{O}/\mathrm{C}$ evolution. Moreover, SOA aging schemes should be better constrained by chamber experiments, particularly with respect to the relationship between the degree of oxygenation and multigenerational aging. Constraining POA $\mathrm{O}/\mathrm{C}$ ratios using source-specific measurements also represents a promising approach for improving the modeled elemental composition of OA.

The linkage between OA volatility and T_g,org (and viscosity) requires revisiting the volatility assignments of existing SOA surrogates (e.g., isoprene-derived SOA in CMAQ) and developing more accurate, dynamic parameterizations of κ_org.