Simulated reductions in heterogeneous isoprene epoxydiol reactive uptake from aerosol morphology in the contiguous United States using the Community Multiscale Air Quality Model (CMAQv5.3.2)

Aerosol particles contain complex mixtures of polar and non-polar species that can undergo organic-inorganic phase separation. In phase-separated aerosol particles, the phase state of the outer organic coating can modulate heterogeneous chemistry. Heterogeneous chemistry leading to isoprene epoxydiol (IEPOX)-derived secondary organic aerosol (IEPOX-SOA) is encoded in the Community Multiscale Air Quality (CMAQ) model and has been the focus of previous aerosol phase separation and phase state work. In a previous study, a constant ratio of water in the organic coating (w_s) was assumed in modeling phase separation and state. Recent studies, however, have highlighted w_s as an important modulator of phase state. This work uses a later CMAQ version (version 5.3.2) with capabilities to model dynamic water uptake to the organic coating – to better predict w_s and its impact on the organic coating phase state. In addition, new parameterizations for estimating organic aerosol phase state were encoded into CMAQ, and were compared with respect to their impacts on phase state and IEPOX-SOA predictions. These evaluations were completed simulating a summertime episode over the continental United States. Simulated diurnal profiles of aerosol phase state agreed within one standard deviation of observationally-derived field measurements. The implementation of phase separation and phase state parameterizations resulted in times and grid cells where IEPOX reactive uptake is completely suppressed. While modelled positive bias in 2-methyltetrol concentrations were decreased with phase separation and phase state updates, modelled methyltetrol sulfates and total IEPOX-SOA concentrations further underpredicted field observations in comparison to Base CMAQ.

1 Introduction

Fine particulate matter (PM_2.5) is a critical component of atmospheric pollution that can impact the climate and human health both directly and indirectly (Cziczo et al., 2013; Kulmala et al., 2013; Pye et al., 2021; Tong et al., 2017). Mitigating these impacts requires knowledge of PM_2.5 composition. The organic fraction of PM_2.5 or organic aerosol (OA) is of particular interest, which, in the Northern Hemisphere, can account for up to 70 % of PM_2.5 mass (Hallquist et al., 2009). Approximately 19 %–93 % of OA is secondary organic aerosol (SOA) (Jimenez et al., 2009) which can form from the partitioning (condensation) of gas-phase organic species to existing aerosol, nucleation mechanisms, or by multiphase chemical processes (Nozière et al., 2015). SOA formation via condensation and nucleation relies on functionalization chemistry that reduces the volatility of gas-phase precursors (Donahue et al., 2006; Odum et al., 1996; Pankow, 1994). Gas-phase SOA precursor species can also dissolve into existing aqueous aerosols and cloud droplets and undergo aqueous-phase oxidation leading to higher or lower volatility species – known as multiphase and/or heterogeneous chemistry (Eddingsaas et al., 2010; Jang et al., 2002; Kurtén et al., 2016; McNeill, 2015; Surratt et al., 2010; Zhang et al., 2018b; Zhang et al., 2019a).

Isoprene epoxydiols (IEPOX) are an isoprene-oxidation product that has been found to participate in heterogeneous formation of SOA (Gaston et al., 2014; Riedel et al., 2016; Surratt et al., 2010). IEPOX-SOA is an important source of SOA (Jo et al., 2021) as its precursor, isoprene, is estimated to be the most abundantly emitted non-methane volatile organic compound (VOC) globally (Guenther et al., 2006; Guenther et al., 2012; Sindelarova et al., 2014). In air masses that have higher hydrogen oxide radicals (HO_x= OH + HO₂) than nitrogen oxide radicals (NO_x= NO + NO₂) concentrations, the formation of IEPOX is favored in comparison to other isoprene SOA intermediates (Lin et al., 2013; Paulot et al., 2009; Surratt et al., 2010). IEPOX-SOA formation has been found to be self-limiting (Riva et al., 2019; Zhang et al., 2019a). Scanning electron microscopy images have indicated that OA, including IEPOX-SOA, can phase separate into organic-rich and inorganic-rich phases (Riva et al., 2019; Zhang et al., 2019a), and this phase separation can limit IEPOX reactive uptake (Gaston et al., 2014; Zhang et al., 2018b). Field campaigns and chamber studies have shown that IEPOX-SOA species predicted by structure-activity relationships to be sufficiently volatile to partition partially out to the gas phase (2-methyltetrols and C₅-alkene triols) remain nevertheless largely in the particle phase (D'Ambro et al., 2019), which could be attributed to the formation of viscous organic coatings (Armstrong et al., 2022; Hu et al., 2016; Yan et al., 2023). Hu et al. (2016) noted that IEPOX-SOA species with higher volatilities (like 2-methyltetrols and C5 alkene triols) should partition partially out to the gas phase, however when relative humidity (RH) was less than 60 %, IEPOX-SOA mass loss to heterogeneous OH oxidation at this RH was largely less than at higher RHs unless OH concentrations were ≥10¹² mol-s cm⁻³. Armstrong et al. (2022) also found that with heterogeneous OH oxidation, there was minimal loss of 2-methyltetrols from IEPOX-SOA, suggesting that their limited evaporation could be attributed to viscous organic coatings. Accounting for phase separation and the phase state of outer organic coatings is paramount to constraining the formation of IEPOX-SOA and its lifetime against heterogeneous OH oxidation.

Modeling of IEPOX-SOA heterogeneous formation has evolved over the years (Budisulistiorini et al., 2017; Pye et al., 2013; Schmedding et al., 2019), with recent attention turned to exploring the impacts of phase separation and phase state on IEPOX reactive uptake (Chen et al., 2024; Gaston et al., 2014; Octaviani et al., 2021; Pye et al., 2017; Schmedding et al., 2020; Zhang et al., 2023). Recently, Schmedding et al. (2020) included and tested the parametrization of a phase separation algorithm in CMAQ and found associated reductions in heterogeneous IEPOX formation coinciding with diffusional limitations in the organic coating. Recent box modeling of chamber experiments that include phase separation found the impacts of phase separation to be sensitive to diffusional limitations in the organic coating (Chen et al., 2024). In both of these studies diffusional limitations were attributed mostly to the phase state of the organic coating (Chen et al., 2024; Schmedding et al., 2020) consistent with experimental findings (Zhang et al., 2018b, 2019a).

To account for phase state of the outer organic coating, Schmedding et al. (2020) made use of Shiraiwa et al. (2017)'s glass transition temperature (T_g) equation (which has been used as a proxy for phase state in previous studies (DeRieux et al., 2018; Koop et al., 2011)). In this study, CMAQ estimated T_g for each individual OA species was calculated from encoded oxygen-to-carbon ratios (O:C) and molar masses (M) of modeled OA species (Schmedding et al., 2020). Individual modeled OA T_g values were then aggregated by mass fraction along with the mass fraction of water associated with the organic coating (w_s) to calculate the overall phase state of the organic coating (T_g(w_org)) assuming ideal mixing (Chen et al., 2023). Although it has been reported that the amount of w_s has a significant influence on the phase state (Lilek and Zuend, 2022; Rasool et al., 2021), in this study, the authors had to assume the unlikely condition that 10 % of all water uptake onto fine aerosol was in the organic coating using CMAQv5.2.1 – as this version only tracked aerosol liquid water (ALW) associated with inorganic aerosols (Schmedding et al., 2020). Despite this assumption, the predicted particle phase state agreed within uncertainty with observationally-derived phase state data from the Centreville, Alabama, supersite during the 2013 SOAS campaign, although, predicting a less viscous phase state than observed (Schmedding et al., 2020; Zhang et al., 2018). Recently, starting with CMAQv5.3.2, OA hygroscopicity parameters (κ_org) – which dictate the amount of water that can be taken up by oxygenated OA – have been encoded allowing for more precise predictions of OA phase state (Pye et al., 2017). Another factor that can potentially impact the estimation of the organic coating phase state is the equation used to determine individual OA T_g – which has yet to be systematically compared.

New relationships between T_g and OA M, O:C, and saturation concentrations (C⁰) – which are properties that commonly define OA species in air quality models – have recently been derived and are summarized in Table 1 (Li et al., 2020; Zhang et al., 2019b). In Li et al. (2020), a T_g equation was fit using a multi-linear least squares regression with O:C and C⁰ used as independent variables, which are used within the two-dimensional volatility basis set framework (Donahue et al., 2011). This T_g parameterization utilized a larger training dataset in comparison to Shiraiwa et al. (2017) which included larger molecular weight OA species, along with sulfated and nitrated species. When this parameterization was evaluated against measured T_g values, the authors reported a correlation coefficient, R=0.93 (Li et al., 2020). Another important process influencing the phase state of an aerosol is the rate at which it is cooled (Zhang et al., 2019b). Aerosols can strengthen vertical updrafts which could potentially result in aerosols cooling at different rates (Abbott and Cronin, 2021). In Zhang et al. (2019b), T_g was measured for 13 OA species at different cooling rates and used to derive a relationship between T_g, M, temperature (T), and C⁰. In comparison to the T_g equation formulated in Shiraiwa et al. (2017), these two new T_g parameterizations account for heavier and more functionalized OA species and therefore warrant exploration in CMAQ (Li et al., 2020; Shiraiwa et al., 2017; Zhang et al., 2019b).

Table 1 The implemented glass transition temperature (T_g) parameterization and variables used such as molar mass (M), oxygen-to-carbon ratios (O:C), temperature (T), and saturation concentration at 298 K (C⁰). Also shown are the experimental data and fit methodology used to develop each parameterization based on carbon-hydrogen (CH), carbon-hydrogen-oxygen (CHO), carbon-hydrogen-oxygen-sulfur (CHOS), and carbon-hydrogen-oxygen-nitrogen (CHON) species.

This study aims to explore the impacts of phase separation and phase state of the organic coating on modeled IEPOX heterogeneous reactive uptake, taking advantage of recently published T_g parameterizations that were derived taking into consideration commonly encoded chemical transport model OA properties (M, O:C, and C⁰) (Table 1) (Li et al., 2020; Shiraiwa et al., 2017; Zhang et al., 2019b). The phase state algorithms, previously used to calculate bulk OA phase state, can now take advantage of water uptake to the organic coating via κ_org parameters released in CMAQv5.3 (Pye et al., 2017). Furthermore, this study explores the influence of modeled aerosol physical properties, such as w_s, C⁰ and O:C on organic coating phase state across parameterizations, along with how each compare with observational data, to provide ranges of OA phase state predicted by CMAQ and its impact on heterogeneous SOA formation.

2 Methods

2.1 Air Quality Model Setup

CMAQ version 5.3.2 was used for this analysis with meteorological inputs and emissions files developed to model the SOAS field campaign from 1 June to 15 July 2013, with a spatial resolution of 12 km by 12 km (Appel et al., 2021) and similar inputs used in Schmedding et al. (2020). Meteorological inputs were predicted using the Weather Research Forecasting Model (WRF) version 3.8 (Skamarock et al., 2008) with lightning assimilation (Appel et al., 2017; Heath et al., 2016). Anthropogenic emissions were sourced from the EPA's National Emissions Inventory (NEI) 2011 version 2. Biogenic emissions were predicted using the Biogenic Emissions Inventory System (BEIS) version 3.6.1, with biogenic isoprene emissions scaled up by a factor of 1.5 to better match isoprene measured at the Centreville, AL SOAS site (CTR) as detailed in (Pye et al., 2017). The State Air Pollution Research Center version 07tic with extended isoprene chemistry and aero7i treatment of SOA (SAPRC07tic_ae7i) was used as the chemical mechanism (Xie et al., 2013) as it explicitly tracks 2-methyltetrols (corresponding to AIETET in the SAPRC07tic_ae7i mechanism) and methyltetrol sulfates (corresponding to AIEOS in the SAPRC07tic_ae7i mechanism, also known as IEPOX organosulfate) (Pye et al., 2013). These are the predominant IEPOX-derived SOA species that permit the tracking of sulfate aerosol's influence on IEPOX reactive uptake (Budisulistiorini et al., 2015, 2017; Pye et al., 2013).

2.2 IEPOX Heterogeneous Reactive Uptake

In CMAQ, IEPOX reactive uptake has been parameterized by the following heterogeneous rate constant (Eddingsaas et al., 2010; Pye et al., 2013):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathrm{IEPOX}}_{\left(g\right)}\to {\mathrm{IEPOXSOA}}_{\left(\mathrm{aerosol}\right)}\text{(2)}& {\displaystyle }{k}_{\mathrm{het}}={\displaystyle \frac{\mathrm{SA}}{\frac{r_{\mathrm{p}}}{D_{\mathrm{g}}}+\frac{\mathrm{4}}{\mathit{\nu }{\mathit{\gamma }}_{\mathrm{IEPOX}}}}}\end{array}$$

where SA is the surface area of the aerosol that IEPOX partitions to, r_p is the aerosol particle radius, D_g is the gas-phase diffusion of IEPOX ($\mathrm{1.9}\times (M_{\mathrm{IEPOX}}{)}^{-\mathrm{2}/\mathrm{3}}$ m² s⁻¹) where M_IEPOX is 118 g mol⁻¹, ν is the mean molecular speed of IEPOX ($\sqrt{\mathrm{8}RT/\mathit{\pi }{M}_{\mathrm{IEPOX}}})$, where R is the ideal gas constant (0.08206 L atm mol⁻¹ K⁻¹) and γ_IEPOX is the following reactive uptake probability coefficient (Eddingsaas et al., 2010; Gaston et al., 2014; Pye et al., 2013; Schmedding et al., 2019; Schmedding et al., 2020):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathrm{IEPOX}}}}& ={\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}+{\displaystyle \frac{\mathit{\upsilon }\cdot r_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{4}{H}_{\mathrm{inorg}}\cdot R\cdot T\cdot D_{\mathrm{a}}\cdot r_{\mathrm{core}}}}\\ & \cdot {\displaystyle \frac{\mathrm{1}}{q\cdot \coth \left(q\right)-\frac{\mathrm{1}}{q}}}\\ & +{\displaystyle \frac{\mathit{\nu }\cdot l_{\mathrm{org}}\cdot r_{\mathrm{p}}}{\mathrm{4}\cdot H_{\mathrm{org}}\cdot R\cdot T\cdot D_{\mathrm{org},\mathrm{eff}}\cdot r_{\mathrm{core}}}}\end{array}\end{array}$$

where α is the unitless accommodation coefficient (0.02), H_inorg is the Henry's law coefficient of IEPOX into the inorganic aqueous core of the aerosol particle (3×10⁷ M atm⁻¹) derived in Nguyen et al. (2014) and used in Budisulistiorini et al. (2017) (Budisulistiorini et al., 2017; Nguyen et al., 2014), D_a is the diffusivity of IEPOX into the inorganic aqueous portion of the aerosol particle (10⁻⁹ m² s⁻¹), and q is the diffuso-reactive length represented by the following equation:

$$\begin{array}{}\text{(4)}& q=r_{\mathrm{p}}\sqrt{{\displaystyle \frac{k_{\mathrm{particle}}}{D_{\mathrm{a}}}}}\end{array}$$

where k_particle is the pseudo-first-order rate constant (s⁻¹) that takes into account the opening of the epoxydiol group on IEPOX by a proton (acid) followed by a nucleophilic attack (by sulfate, water, or monomer IEPOX-SOA species) on the free carbo-cation formed from this ring opening reaction, and is represented by the following equation (Eddingsaas et al., 2010; Pye et al., 2013):

$$\begin{array}{}\text{(5)}& k_{\mathrm{particle}}=\sum _{i=\mathrm{1}}^N\sum _{j=\mathrm{1}}^M{k}_{i,j}\left[{\mathrm{nuc}}_i\right]\left[{\mathrm{acid}}_j\right]\end{array}$$

where k_i,j are individual acid-nucleophile rate constants defined by Table S3. A slight alteration was made to Eq. (5) in Base CMAQ and all of the phase state sensitivity model runs in that we use the individual third-order rate constant for AIEOS from Riedel et al. (2016), instead of the default from Piletic et al. (2013), to be consistent with Schmedding et al. (2020).

When phase separation occurs, r_core is calculated by the following equation (Schmedding et al., 2020):

$$\begin{array}{}\text{(6)}& l_{\mathrm{org}}=r_{\mathrm{p}}-r_{\mathrm{core}}\end{array}$$

The first two terms of the γ_IEPOX reactive uptake probability equation (Eq. 3) represent that of an assumed homogeneously mixed aerosol. With assuming phase separation, a third term (the organic coatings resistor term) was added to this equation to represent the potential resistance to IEPOX reactive uptake posed by the organic coating in a phase-separated aerosol particle where l_org is the organic coating thickness, H_org is the Henry's Law coefficient dictating the dissolution of IEPOX into the organic coating (2×10⁶ M atm⁻¹) (Zhang et al., 2018b), D_org,eff is the diffusivity of IEPOX through organic coating (Fig. 1, Eq. 7), and r_core is the radius of the inorganic aqueous core. When there is no phase separation l_org=0, cancelling out the third resistor term and r_core=r_p thus reverting Eq. (3) to the original parameterization of IEPOX reactive uptake currently in CMAQ (Eddingsaas et al., 2010; Lin et al., 2013). It is important to note that we do not simulate water movement between the aqueous inorganic core and the organic coating, and therefore are unable to capture the impacts of dilution of acids and nucleophiles involved in IEPOX reactive uptake.

2.3 Implementation of phase separation and phase state

A summary of the CMAQ algorithm used for determining phase separation and phase state are shown in Fig. 1, mirroring the implementation of the “PhaseSep2” model setup documented in (Schmedding et al., 2020), with the exception of using individual T_g from multiple studies (Li et al., 2020; Shiraiwa et al., 2017; Zhang et al., 2019b). As shown in Fig. 1, phase separation was determined based on the separation relative humidity (SRH) and occurs when the SRH ≥ RH (Bertram et al., 2011; You et al., 2014). The SRH was determined based on aggregated aerosol O:C (O:C_avg) and organic matter to inorganic sulfate ratios (OM:IN_sulf) (Bertram et al., 2011; Schmedding et al., 2020; Song et al., 2018; Zuend and Seinfeld, 2012). The O:C were derived from organic matter-to-organic carbon ratio (OM:OC) using the relationship published in Simon and Bhave (2012). If phase separated, either the Shiraiwa, Zhang, or Li T_g equations (explained in Sect. 2.4) are used to calculate the T_g of individual OA species. The individual T_gs are aggregated by source type (anthropogenic or biogenic) and then Eq. (8) (shown in Fig. 1) is used to calculate the overall T_g of the organic coating accounting for water, referred to from here-on-out as T_g(w_org). The T_g(w_org):T ratios were then used to determine the viscosity of the organic coating (η_org) (Eqs. 10–12 shown in Fig. 1) using a modified Vogel-Tamman-Fulcher equation (Angell, 1991; DeRieux et al., 2018; Fulcher, 1925; Schmedding et al., 2020; Tammann and Hesse, 1926; Vogel, 1921). The calculated η_org was then used to determine the diffusivity of IEPOX through the organic coating (D_org,eff, Eq. 7).

Figure 1 Algorithm used to determine phase separation (Bertram et al., 2011; You et al., 2013, 2014; Zuend and Seinfeld, 2012) and phase state (Schmedding et al., 2020), including the parameters used to determine whether aerosol particles were in a liquid, semi-solid or solid phase state. Also shown are the equations used to determine viscosity (η_org) (Fulcher, 1925; Tammann and Hesse, 1926; Vogel, 1921) and the effective diffusivity of IEPOX through the organic coating (D_org,eff) (Miller, 1924).

When w_s=0, the aggregated OA phase state equation (Eq. 8) collapses into the mass fraction weighted aggregated dry phase state equation, represented by T_g,org (Dette et al., 2014; Li et al., 2021, 2020):

$$\begin{array}{}\text{(14)}& T_{\mathrm{g},\mathrm{org}}=w_{\mathrm{a}}{T}_{\mathrm{g},\mathrm{a}}+w_b{T}_{\mathrm{g},\mathrm{b}}\end{array}$$

Given there was never an occurrence during our modeling episode where w_s=0, we use Eq. (14) when w_s≤0.1 (Rasool et al., 2021), signifying a dry aggregate glass transition temperature that does not consider water uptake to the organic shell. While our implementation of phase separation and phase state in CMAQ did not implement in-model conditionals to switch between using T_g(w_org) and T_g,org based off of these w_s thresholds, we identify dry aerosol phase state instances offline in our analysis

2.4 CMAQ Implementation of Glass Transition Temperature

The details of the Shiraiwa T_g equation and its formulation is provided in Schmedding et al. (2020) and Shiraiwa et al. (2017):

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{T}_{\mathrm{g},i,\mathrm{S}}& =-\mathrm{21.57}+\mathrm{1.51}M-\mathrm{0.0017}{M}^{\mathrm{2}}+\mathrm{131.4}\left(\mathrm{O}:\mathrm{C}\right)\\ & -\mathrm{0.25}M\left(\mathrm{O}:\mathrm{C}\right)\end{array}\end{array}$$

The details of the Zhang et al. (2019b) and Li et al. (2020) T_g equation and their formulation are provided below.

2.4.1 Zhang parameterization

A logarithmic relationship between the relaxation time of an aerosol species and inverse temperature has been established in Zhang et al. (2018a). Using this relationship, Zhang et al. (2019b) explored the implication of cooling rates on T_g values. Ultimately, a cooling rate of 5 K min⁻¹ was determined to be the most atmospherically-relevant cooling rate for calculating the T_g values of individual aerosol species using the Zhang et al. (2019b) parameterization ($T_{\mathrm{g},i,\mathrm{Z}}$) at average atmospheric conditions. Thirteen measured $T_{\mathrm{g},i,\mathrm{Z}}$ values were then related to saturation concentrations predicted using the EVAPORATION model (Compernolle et al., 2011) and the following equation :

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{T}_{\mathrm{g},i,\mathrm{Z}}& =\mathrm{480.071}\\ & -{\displaystyle \frac{\mathrm{54395}}{({\log }_{\mathrm{10}}\left(\frac{RT}{M{\mathrm{10}}^{\mathrm{6}}}{C}^{\mathrm{0}}\right)-\mathrm{1.7929}{)}^{\mathrm{2}}+\mathrm{116.49}}}\end{array}\end{array}$$

where R is the ideal gas law constant ($\mathrm{8.2}\times {\mathrm{10}}^{-\mathrm{5}}$ m³ atm mol⁻¹ K⁻¹), T is the temperature (K), M is molar mass of the organic compounds, and C⁰ is the saturation concentration (µg m⁻³).

2.4.2 Li parameterization

The Li et al. (2020) parameterization used a database of 2448 CH, CHO, CHON, and CHOS compounds, of which 943 are sulfated and 276 are nitrogenated species. Measured T_g values are available for a total of 337 compounds, with the majority (of 259) being CHO compounds. When T_g measurements are not available, T_g is estimated from the melting temperature (T_m) by applying the Boyer-Kauzmann rule of $T_{\mathrm{g}}=g\times T_{\mathrm{m}}$ with g≈0.7 (Koop et al., 2011), referred to as the “estimated T_g” (Li et al., 2020). The authors also used a test dataset of 654 CHO and 212 CHON compounds that were not included in the training dataset. Using experimentally measured or estimated T_g values for aerosol species in a multi-linear least squares analysis was completed with C⁰ values and O:C ratios, and the authors reported the following equation (Li et al., 2020):

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{T}_{\mathrm{g},i,\mathrm{L}}& =\mathrm{289.10}-\mathrm{16.50}\times {\log }_{\mathrm{10}}\left(C^{\mathrm{0}}\right)-\mathrm{0.29}\\ & \times {\left[{\log }_{\mathrm{10}}\left(C^{\mathrm{0}}\right)\right]}^{\mathrm{2}}+\mathrm{3.23}\times {\log }_{\mathrm{10}}\left(C^{\mathrm{0}}\right)\left(\mathrm{O}:\mathrm{C}\right)\end{array}\end{array}$$

This parameterization has been recently applied in a global chemical transport model GEOS-Chem and a regional air quality model WRF-Chem, providing consistent results with phase state measurements at the surface level (Luu et al., 2025; Zhang et al., 2024)

2.5 CMAQ Modifications for Calculation of T_g

The T_g implementations described above require C⁰ (at 298 K) for all CMAQ OA species. These values are shown in Table S1 in the Supplement. It is important to note that some of the CMAQ low-volatility organic species (LVOS) are assigned a nonvolatile saturation concentration ($C^{\mathrm{0}}={\mathrm{10}}^{-\mathrm{9}}$ (µg m⁻³)) to ensure that these species condense entirely, consistent with the initial default implementation (Pye et al., 2015, 2023). With the exception of dimers and oligomers, there has been evidence to suggest that some of these LVOS species, like monoterpene hydrolysis products (AMTHYD) and AIETET, may still partition back out to the gas phase (Budisulistiorini et al., 2017; Kurtén et al., 2016). To avoid biasing modeled T_g which depended on C⁰ values, all LVOS that had a C⁰ of 10⁻⁹ µg m⁻³ were updated in CMAQ as shown in Table S2. The process for updating C⁰ values started with a literature review for any reported C⁰ or saturation vapor pressure (p₀) that could replace the existing CMAQ species' C⁰. If neither a C⁰ nor p₀ was found, the chemical structure of the LVOS surrogates were referenced and their structures were translated into the Simplified Molecular Input Line Entry System (SMILES) coding notation (Anderson et al., 1987). Once SMILES were coded, they were input into the OPEn structure–activity/property Relationship App (OPERA) (Mansouri et al., 2018), where p₀s were estimated and then used to calculate C⁰s using the following relationship (Donahue et al., 2006; Pankow, 1994; Zhang et al., 2019b):

$$\begin{array}{}\text{(18)}& C^{\mathrm{0}}={\displaystyle \frac{M{\mathrm{10}}^{\mathrm{6}}{p}_{\mathrm{0}}}{RT}}\end{array}$$

where T is assumed to be 300 K to correspond with temperature used to derive all CMAQ C⁰s (Pye et al., 2017). CMAQ M's were not updated to correspond with the proxy species' p₀ and instead reflect that of derived values from the SAPRC07tic_ae7i mechanism (for conservation of mass) and represent a potential limitation of this study.

2.6 Observational Data

Observed O:Cs and Ms were obtained from 800+ OA species measured at the Centreville, Alabama (CTR) site using a high-resolution time-of-flight chemical ion mass spectrometer (HR-ToF-CIMS) coupled with a filter inlet for gases and aerosols (FIGAERO) and a two-dimensional gas chromatography time-of-flight mass spectrometer (GCxGC-ToF-MS) (Zhang et al., 2018). Observed ($T_{\mathrm{g},\mathrm{org}}:T)$, signifying the dry phase state, were estimated using Shiraiwa et al.'s (2017) T_g parameterization from O:Cs and Ms, and were collected at the CTR site (Zhang et al., 2018). The dry $T_{\mathrm{g},\mathrm{org}}:T$ was estimated from OA observations given the absence of ALW in measurements, and the Shiraiwa parameterization was used due to the lack of C⁰ measurements for OA at this site (Zhang et al., 2018). IEPOX-SOA measurements were taken from the CTR and Look Rock, Tennessee (LRK) 2013 SOAS field campaign sites (Budisulistiorini et al., 2017; Shiraiwa et al., 2017; Zhang et al., 2018). IEPOX-SOA was measured at the CTR site using an aerosol mass spectrometer and positive matrix factorization (Hu et al., 2016; Ulbrich et al., 2009). IEPOX-SOA at the LRK site was measured using gas chromatography/electron ionization mass spectrometry (GC/EI-MS) and ultra-performance liquid chromatography/diode array detection-electrospray ionization-high-resolution quadrupole time-of-flight mass spectrometry (UPLC/DAD-ESI-HR-QTOFMS) with resolved AIETET and AIEOS measurements (Budisulistiorini et al., 2017).

3 Results

3.1 Phase separation frequency across space and altitude

The frequency of aerosol phase separation predicted at different model layers across the U.S. is shown in Fig. 2. At the surface (Fig. 2a), phase separation occurs 90 %–100 % of the time in the Southeastern U.S. and Western U.S. which can be attributed to lower O:C_avg (ranging from 0.4–0.6) in the Southeast and lower RH (ranging between 10 %–40 %) in the West. Over the oceans and Great Lakes, phase separation frequency decreases due to both high RH (≥ 80 %) and high O:C (≥ 0.7). Phase separation frequencies in layer 18 (representing the lower troposphere) (Fig. 2b) decrease over the eastern U.S. largely corresponding to increases in O:C_avg (by ∼ 0.1 on average) where RH increases on average by ∼ 5 %. Phase separation continues to increase in layers 28 (representing the upper troposphere) (Fig. 2c) and 35 (representing the lower stratosphere) (Fig. 2d) corresponding with continued decreases in RH with altitude.

Figure 2 Frequency of modelled aerosol phase separation for the SOAS 2013 Field Campaign episode (1 June–15 July 2013) at layer 1 (surface layer; 1–0.9975 atm (a)), layer 18 (corresponding to ∼ 1.8 km; 0.82–0.8 atm (b)), layer 28 (corresponding to ∼ 8 km; 0.4–0.35 atm (c)), and layer 35 (corresponding to ∼ 17 km; 0.05–0 atm (d)) across hourly model estimations. Black open circles in panel (a) correspond to Look Rock, Tennessee (LRK) and Centreville, Alabama (CTR) SOAS observation sites.

The aggregated phase separation frequency across the first layer in this work was estimated to be 71.4 % which compares well with 68.5 % predicted in Schmedding et al. (2020) though slightly higher (Schmedding et al., 2020). This difference can be attributed to a slightly lower O:C_avg with the SAPRC07tic_ae7 mechanism compared to the CB6R3 mechanism used in Schmedding et al. (2020). At the CTR SOAS site, the phase separation frequency is ∼ 97.9 %, significantly higher than that reported in Schmedding et al. (2020) (i.e., 65.4 %) and in Pye et al. (2017) (i.e., 79.1 %). This suggests that phase separation predicted in this study at the CTR SOAS site may represent an upper bound (Pye et al., 2017; Schmedding et al., 2020).

3.2 Impacts of modeling dynamic water uptake to the organic coating on phase state

One major advancement in CMAQv5.3.2 is the addition of hygroscopicity parameters to track the uptake of water to modelled organic species (Appel et al., 2021; Pye et al., 2017). In comparison to Schmedding et al. (2020) we find that this update increases w_s by up to 0.05–0.25 (Fig. S2) over the Eastern and northern part of the domain at the surface (layer 1), resulting in decreases in T_g(w_org):T by 0.05–0.15 (Fig. S3). The w_s largely decreases over the oceans in this study compared to Schmedding et al. (2020) and can be partially attributed to increased ALW over the oceans (which result in higher w_s in Schmedding et al. (2020); w_s≥0.6) and relatively lower OA concentrations over the oceans (compared to land) that have resulted in less hygroscopic water uptake ($\mathrm{0.2}\le w_{\mathrm{s}}\le \mathrm{0.6}$) to the organic coating in our model runs. Differences in w_s decrease with height, given the lack of humidity at higher altitudes (Fig. S2c–d), however, are still higher in Schmedding et al. (2020) resulting in increases in T_g(w_org):T in our model runs in comparison (Fig. S3c–d).

3.3 Phase state of the organic coating across space and altitude

Following criteria used in Schmedding et al. (2020), a liquid phase organic coating is predicted when $T_{\mathrm{g}}\left(w_{\mathrm{org}}\right):T<\mathrm{0.8}$ or $T_{\mathrm{g},\mathrm{org}}:T<\mathrm{0.8}$, a semi-solid phase organic coating when $\mathrm{0.8}\le T_{\mathrm{g}}\left(w_{\mathrm{org}}\right):T<\mathrm{1}$ or $\mathrm{0.8}\le T_{\mathrm{g},\mathrm{org}}:T<\mathrm{1}$, and a solid phase organic coating when $T_{\mathrm{g}}\left(w_{\mathrm{org}}\right):T\ge \mathrm{1}$ or $T_{\mathrm{g},\mathrm{org}}:T\ge \mathrm{1}$. The frequency of time that the organic coating of an aerosol is in a liquid phase in layer 1 (i.e., ground level, GL) for the duration of this modeling episode when the amount of ALW present in the organic coating is not capped (all w_s) and when w_s≤0.1, corresponding to T_g(w_org):T and $T_{\mathrm{g},\mathrm{org}}:T$, is shown in Fig. 3. The Southeast U.S. has the highest frequency of liquid phase organic coatings in all three parameterizations when all fractions of ALW are present (all w_s) (Fig. 3a–c) in the organic coating for all model simulations, with a distinct transition going from east to west (Fig. 3a–c). The Southeast U.S. is known to experience hot and humid summers with both increased temperatures and increased ALW, which result in the expected plasticizing of the organic coating. In this study we only simulate Summer 2013, however, expect that this spatial trend in organic coating phase state calculated may be different during other seasons with changes in emissions, RH, and T.

Figure 3 Frequency that the organic coating is in the liquid phase state during the SOAS 2013 Field Campaign episode (1 June–15 July 2013) for the (a, d) Shiraiwa, (b, e) Zhang and (c, f) Li parameterizations for all w_s in the organic coating (a–c) and when w_s≤0.1 in the organic coating (d–f) across hourly model estimates. Black open circles in each panel correspond to Look Rock, Tennessee (LRK) and Centreville, Alabama (CTR) SOAS observation sites.

When w_s>0.1 the Li and Zhang simulations predicted a higher frequency of semi-solid and solid organic coatings across the U.S in comparison to the Shiraiwa simulation, with the liquid-to-solid transition zone moving further west (Fig. 3a–c). Zhang et al. (2018b) found that IEPOX uptake is phase state limited ∼ 40 % of the time in the Southeast U.S, which agrees qualitatively with the Zhang and Li simulations. The western U.S. shows distinctly less frequent liquid organic coatings (when w_s>0.1) for all simulations (Fig. 3a–c) corresponding with decreased RH.

The T_g(w_org):T of the organic coating also increases with altitude, as shown in Fig. 4. Across all model runs the percentage of time when w_s <0.1 was ∼ 31 %, 47 %, 76 %, and 100 %, corresponding with decreases in RH with increases in altitude. In the surface layer of the model (Layer 1; @GL), across all grid cells and hourly timesteps, the Shiraiwa simulation predicts the highest frequency of the liquid phase state, in agreement with Fig. 3a. The Zhang and Li simulations show a slightly lower liquid phase state frequency. In layer 18, the Shiraiwa simulation also predicts the highest frequency of the liquid phase state (∼ 50 % of the time) in comparison to the Zhang and Li simulations. While the phase state of the organic coating at the surface and lower troposphere (layer 18) are important for heterogeneous chemistry, the phase state of the organic coating in the upper troposphere (layer 28) and lower stratosphere (layer 35) are important for cirrus cloud formation (Berkemeier et al., 2014; Murray et al., 2010; Wolf et al., 2020). In all simulations at layers 28 and 35, the phase state of the organic coating is mostly semi-solid to solid given lower RH and lower T at higher altitudes. This phase state can potentially provide a surface for heterogeneous ice nucleation via deposition (Berkemeier et al., 2014; Murray et al., 2010; Wolf et al., 2020).

Figure 4 The normalized frequency of occurrence of organic coating phase states (when w_s>0.1; (T_g(w_org):T)) across the surface spatial domain and for the entire simulation period (hourly model estimations) for modeled layers 1 (blue, GL = ground level), 18 (green, 1.8 km a.g.l. = above ground level), 28 (yellow, 8 km a.g.l.), and 35 (red, 17 km a.g.l.). Shown are the estimates using the Shiraiwa (a), Zhang (b), and Li (c) simulations. Distributions were normalized so that the area under each curve was equal to 1. Dotted blue lines representing the bounds for the three different phase states: liquid, semi-solid, and solid.

Recent studies have indicated that ALW is the most influential factor on aerosol phase state (Li et al., 2021; Rasool et al., 2021). The impact that ALW in the organic coating has on phase state is demonstrated by the comparison to the frequency of the liquid phase state when w_s≤0.1 (Fig. 3d–f). While all simulations predict some liquid organic coatings in the Eastern U.S. when w_s >0.1 (Fig. 3a–c), both the Zhang and Li parameterizations almost never predict a liquid phase organic coating at w_s≤0.1 (Fig. 3d–f). However, the frequency of the liquid phase state for the Shiraiwa model run is still ≥80 % in the Southeastern U.S., whereas the Li and Zhang simulations are 90 %–100 % in the semi-solid phase state. While ALW is a significant modulator of the organic coating's phase state, the composition and T_g equation used can also impact the phase state.

3.4 The influence of T_g equations on organic coating phase state

Phase state differences across the three model parameterizations are partially due to differences in individual OA component's T_g values. The T_g values predicted by the Shiraiwa equation for each OA species along with differences between the Shiraiwa, Zhang and Li equations are shown in Table 2. Since the Zhang equation takes into consideration ambient temperatures in predicting T_g, the comparison was made at 298 K. The largest differences in individual T_g where those of the isoprene organic nitrates (AISOPNN) between the Shiraiwa and Zhang equations (Table 2) differing by 143 K. The Shiraiwa equation predicts a solid T_g at normal surface temperatures, while the Zhang parameterization predicts a semi-solid or liquid phase. AISOPNN has a higher M and O:C, which would increase its T_g, however a higher C⁰ relative to most OA species, which would decrease its T_g, and therefore would contribute to these differences.

Table 2 Glass transition temperatures (T_g) predicted for each OA species in the SAPRC07tic_ae7i mechanism using the Shiraiwa T_g equation, differences between the Zhang and Shiraiwa equations, the Li and Shiraiwa equations, and between the Zhang and Li equations. Species descriptions can be found in Table S1 in the Supplement.

For most species, the Zhang and Li equations predict a higher T_g than the Shiraiwa equation, and the Li equation predicts a higher T_g than the Zhang equation (Fig. S4). There are some exceptions to this trend where the Shiraiwa T_g equation predicts a higher T_g than both Zhang and Li for AISOPNN, methacrolein epoxide derived organosulfate (AIMOS), and AIEOS (Table 2), due to their higher O:C (Table 2, Fig. S4a–b). The relationship between the Zhang and Li T_g predictions are linear (Fig. S4c, f, and i) in agreement with findings in Li et al. (2020) with higher volatility species having a higher T_g predicted by Li in comparison to that predicted by Zhang (Fig. S4f). Two CMAQ species of note that have a much higher predicted T_g by Zhang and Li than by Shiraiwa are potential combustion SOA (APCSO) and glyoxal- and methylglyoxal-derived aerosol (AGLY). The M of APCSO is 170 g mol⁻¹ and the O:C is 0.67, however it has the lowest C⁰ of all the OA species (0.00001 µg m⁻³) which would be the main influence on predicting T_g in the Zhang and Li equations. With the Shiraiwa equation, AGLY has a T_g of 160 K (meaning the ambient temperature would need to be 200 K for it to be in a semi-solid phase and 160 K for it to be in a solid phase state) and therefore at normal surface temperatures for this modeling episode and domain, AGLY would always be in a liquid phase state. It should be noted that the M used for AGLY is 66.4 g mol⁻¹, while the surrogate species used to estimate the C⁰ for AGLY is ammonium oxalate – a common particle-phase product of both glyoxal and methylglyoxal – which has a published low vapor pressure of $\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{6}}$ Pa (Paciga et al., 2014). If the M of ammonium oxalate's true molar mass (124 g mol⁻¹) were to be used instead, the T_g predicted for AGLY in the Shiraiwa parameterization would increase to 216 K, which could increase the T_g(w_org) or the T_g,org for this simulation.

As shown in Fig. S5, the model predicts an aerosol composition mainly influenced by CMAQ biogenic species in the Southeastern U.S. and the Western U.S. Figure S6 shows that the biogenic species with the highest contribution to biogenic aerosol mass in these regions are monoterpene-derived hydration products (AMTHYD), and low-volatility monoterpene SOA species (AMT1 and AMT2). As shown in Table 2, while the T_g predicted by the Shiraiwa and Zhang equations agree within 2 K for AMTHYD, the Li equation predicts a T_g 34–36 K higher. Nonetheless, for this modeling episode and domain, T_g's of 226–262 K would still likely indicate a liquid phase state at the surface. Differences in T_g for low volatility monoterpene derived SOA (AMT1 and AMT2) may impose more differences in phase state given the range in predicted T_g for AMT1 is 299–318 K and the T_g for AMT2 ranges from 243–305 K.

In the Northern Midwest, the leading contributor to aerosol mass were anthropogenic species (Fig. S5), and the leading anthropogenic species were semi-volatile oxidized combustion organic products (ASVOO1) and low-volatility oxidized combustion organic products (ALVOO2) (Fig. S6a and b). For ASVOO1, T_g ranges from 207–289 K, while for ALVOO2, T_g ranges from 222–303 K (Table 2). These species and their ranges of T_g are likely the cause for the differences in the liquid-to-solid transition in the Northern Midwest (Fig. 3a–c) across all model simulations.

3.5 Comparison of modelled $T_{\mathrm{g},\mathrm{org}}:T$ to derived observations

Figure 5 shows the 2013 SOAS CTR ground site estimated $T_{\mathrm{g},\mathrm{org}}:T$ diurnal profile and the predicted $T_{\mathrm{g},\mathrm{org}}:T$ (when w_s≤0.1) diurnal profiles for each parameterization at that grid cell, with vertical fliers representing one standard deviation. It should be noted that observations did not account for ALW, and thus, it would be inappropriate to compare them with T_g(w_org):T. The Shiraiwa parameterization predicted the lowest $T_{\mathrm{g},\mathrm{org}}:T$ at the CTR site. The Zhang and the Li predictions of $T_{\mathrm{g},\mathrm{org}}:T$ were higher with similar diurnal profiles. Modeled O:C_avg are biased negatively, having a normalized mean bias (NMB) of ∼ −0.25, compared to observed O:C_avg at the SOAS CTR site, and are lower than O:C_avg predicted in Schmedding et al. (2020) (Zhang et al., 2018). While this may influence model-measurement comparisons, O:C ratios were found to have a minor influence of phase state in comparison to M and w_s (Koop et al., 2011; Li et al., 2020; Shiraiwa et al., 2017). As noted previously (Sect. 2.4), there remains a disconnect between modeled M and modeled C⁰ in CMAQ which could explain why the Shiraiwa predicted $T_{\mathrm{g},\mathrm{org}}:T$ were biased low and lower than the other two parametrizations. Estimations of measured $T_{\mathrm{g},\mathrm{org}}:T$ at the CTR SOAS site (Fig. 5) may also be biased low given that C⁰ values were not measured.

Figure 5 Diurnal comparisons of model predicted and observed $T_{\mathrm{g},\mathrm{org}}:T$ for the Shiraiwa, Zhang and Li model simulations at the CTR SOAS site for 1 June 2013–15 July 2013. Vertical flyers on observations represent 1 standard deviation.

3.6 Impacts on IEPOX-SOA Modeled Formation

Episode-averaged IEPOX reactive uptake coefficients, γ_IEPOX, were compared across simulations (Fig. 6). The episode-averaged base model predictions of γ_IEPOX ranged from 10⁻⁶–10⁻³, with the average value at the LRK SOAS site being $\mathrm{6.21}\times {\mathrm{10}}^{-\mathrm{4}}$, which is in agreement with episode-averaged value from Budisulistiorini et al. (2017) (i.e. $\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{4}}$), though slightly higher. The diffusional hindrance to heterogeneous IEPOX reactive uptake due to existing viscous organic coatings had the largest impact in the Western part of the U.S., reducing γ_IEPOX by 80 %–99 % (Fig. 6). In most of the Western U.S., the abundance of IEPOX is limited (Fig. S1), with the exception of coastal Northern California where these decreases can potentially impact model performance in estimating IEPOX-SOA and total OA (Fig. 6). IEPOX concentrations are also appreciable in parts of Canada (Fig. S1) where γ_IEPOX reductions can be >0.001 (representing up to a 99.99 % reduction). In the Southeast U.S. (where IEPOX concentrations peak (Fig. S1)), predicted reductions in γ_IEPOX range from 20 %–80 % depending on the phase state simulation used (Fig. 6). At the LRK SOAS site, γ_IEPOX were reduced to values ranging from $\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{4}}$ to $\mathrm{4.84}\times {\mathrm{10}}^{-\mathrm{4}}$ and are in closer agreement to values in Budisulistiorini et al. (2017). The D_org,eff across all simulations ranges between 10⁻¹²–10⁻¹³ m² s⁻¹ and are in agreement with box modeling of resistance to IEPOX reactive uptake (Chen et al., 2024). Increases in γ_IEPOX occur over the oceans, are up to ∼ 0.0001, and account for up to a ∼ 7 % increase. These bounds can be considered a consequence of numerical noise along with decreases of the same magnitude and are not likely to have a major impact on IEPOX-SOA formation. The phase state and phase separation model run-times are ∼ 5 % slower in all of the sensitivity runs compared with the base.

Figure 6 Average difference in IEPOX-SOA reactive uptake coefficients (γ_IEPOX) from the base (no organic coatings resistor term) simulation compared to γ_IEPOX predicted by Shiraiwa (a), Zhang (b), and Li (c) model simulations. Average percent changes in γ_IEPOX from the base model simulation compared to γ_IEPOX predicted by Shiraiwa (d), Zhang (e), and Li (f) model simulations. White spaces represent areas where the average differences are below −0.001.

With the implementation of phase separation, water equilibrium between the inorganic core and the organic coating was not considered. In reality, water taken up by the organic shell could make its way into the inorganic core given its polar attraction, diluting both acidity and nucleophile concentrations, and therefore, reducing k_particle. This is important to note given γ_IEPOX's sensitivity to k_particle. When increasing k_particle from 10⁻³ to 10⁻² s⁻¹, γ_IEPOX increases by ∼ 740 % of magnitude (assuming an organic coating radius of 50 nm and T=298 K) (Eq. 3). An increase in D_org,eff from 10⁻¹³ to 10⁻¹² m² s⁻¹ results in an increase in γ_IEPOX by ∼ 0.02 % (assuming an organic coating radius of 50 nm and T=298 K) (Eq. 3). The impacts of internal water movement in phase separated aerosol requires further exploration, as it may be another important modulator of IEPOX heterogeneous reactive uptake (Zhang et al., 2019a), and other heterogeneous reactive uptake reactions (Farrell et al., 2025).

Total predicted IEPOX-SOA concentrations were compared between each parameterization (including base CMAQv5.3.2) and the observations obtained at the CTR SOAS site (Hu et al., 2016) shown in Fig. 7a. The base CMAQ run had the best model performance in simulating total IEPOX-SOA, yet still underpredicted observations, with an NMB of −0.37. Model performance in predicting IEPOX-SOA concentrations worsened when including phase separation and phase state with NMB's of −0.48, −0.53, and −0.56 for the Shiraiwa, Zhang, and Li model simulations.

Figure 7 Observational and predicted concentrations of (a) total IEPOX-SOA at CTR and (b) AIETET and (c) AIEOS at the LRK SOAS site for 1 June–15 July 2013. Predicted concentrations are shown for base CMAQ (grey), Shiraiwa (green), Zhang (red), and Li (blue) model simulations.

Observations of AIETET and AIEOS (Budisulistiorini et al., 2015) were compared with corresponding CMAQ species concentrations (Fig. 7b–c) The LRK site is situated in the Great Smoky Mountains, where the altitude (∼ 800 m above sea level) can affect both T and RH. CMAQ predicts Ts ranging from 290–302 K at this site for the duration of the SOAS campaign, with predicted w_s values ranging from 0.1–0.6. Differences in IEPOX-SOA predictions across model runs at the LRK site coincide with both low Ts and low w_s. While we do not evaluate T and RH for our models in this paper, both variables were indirectly nudged in WRF simulations to observations of soil moisture (Gilliam et al., 2006; Pleim and Gilliam, 2009; Pleim and Xiu, 2003). AIETET is overpredicted most of the time for all parameterizations (including the base run), with the base run having a NMB of 1.54. For the model runs accounting for phase separation and phase state, the bias improves with NMBs of 1.26, 1.19, and 1.14 for the Shiraiwa, Zhang, and Li parameterizations, respectively. It is important to note that recent studies have suggested AIETET can re-volatilize (Riedel et al., 2015; Su et al., 2025); however, heterogeneous oxidation experiments of IEPOX-SOA have also suggested aged IEPOX-SOA (which produces more functionalized oligomer species) may limit this re-volatilization (Armstrong et al., 2022; Hu et al., 2016; Yan et al., 2023). The base CMAQ run underpredicts AIEOS concentrations with a NMB of −0.66. All model runs accounting for phase separation and phase state underpredict AIEOS concentrations further with NMBs of −0.73, −0.75, and −0.76 for the Shiraiwa, Zhang, and Li parameterizations, respectively. These underpredictions of AIEOS coincide with underpredictions in sulfate concentrations (Fig. S7), which may decrease modeled branching of IEPOX-SOA towards AIEOS. High AIEOS concentrations at the LRK site (Fig. 7c) coincide with high inorganic sulfur concentrations (Fig. S7). and the formation of IEPOX-derived organosulfates has previously been shown to depend on availability of inorganic sulfate (Brüggemann et al., 2020; Chen et al., 2021; Piletic et al., 2013). There is also uncertainty in branching ratios between AIETET and AIEOS – arising from the reported variability of sulfate ion activity (Petters et al., 2021). Updating sulfate formation mechanisms to include heterogeneous sulfur chemistry in ALW (Farrell et al., 2025), along with increasing the AIEOS branching fraction, may help resolve this bias for all model runs (including the Base) (Budisulistiorini et al., 2017; Chen et al., 2024). Organosulfates have been recently predicted and demonstrated to be surface active species, and thus, their incorporation into total ion activity could change previous branching ratio parameterizations (Hyttinen et al., 2020; Olson et al., 2019; Riva et al., 2019). Overall, we were confined to measurements from two sites, and thus our model performance metrics are statistically underpowered. More measurements of IEPOX-SOA, and measurements of speciated IEPOX-SOA would help better inform model performance in the future.

4 Discussion

4.1 Modulators of organic coating phase state

Model predictions show that both the w_s and the choice in T_g parameterization impact the phase state of the organic coating on phase-separated aerosols. Previous studies have shown that ALW has the largest influence on the phase state of organic aerosols (DeRieux et al., 2018; Rasool et al., 2021). In comparison to assuming 10 % of ALW is in the organic coating, we find in using a w_s that is based on OA chemical properties, that organic coatings are more liquid in comparison to Schmedding et al. (2020) over the Eastern and Northern parts of the contiguous U.S. domain and therefore may not pose as much of a diffusional barrier to heterogeneous reactive uptake (Schmedding et al., 2020). The assumption that water uptake to the organic coating, as parameterized, ignores the phase state may introduce bias, however, to our knowledge a consistent relationship between organic aerosol phase state and hygroscopic growth has yet to be explored (Diveky et al., 2021; Hodas et al., 2015; Lienhard et al., 2015; Pajunoja et al., 2016; Parsons et al., 2004; Tan et al., 2024), as opposed to the relationship between O:C ratio and hygroscopic growth (Jimenez et al., 2009; Lambe et al., 2011; Massoli et al., 2010; Pajunoja et al., 2015).

This study is the first to compare the impacts of different T_g parameterizations on phase state in a chemical transport air quality model (to our knowledge), providing insight into the ranges of organic coating phase states that can occur. To disentangle the effects of T_g parameterizations and the effect of ALW, the phase state of the organic coating was analyzed when w_s ≤0.1 (Fig. 3d–f). When w_s ≤0.1 the Shiraiwa run still predicted a liquid phase state in the Southeast U.S., whereas the Zhang and Li simulations predicted a semi-solid organic coating. This difference in predictions can partially be explained by inconsistencies in CMAQ species properties and differences in individual T_g parameterizations (leveraging different species properties (Table 1)). For example, in the case of AGLY, the CMAQ modeled species M implies a much smaller compound than the volatility (C⁰) suggests. As a result, the estimated AGLY T_g varies by a factor of 2 across the different model simulations. With the updates to C⁰ in this work, the Ms were not updated to reflect the proxy species used, and therefore, may introduce bias in the prediction of T_g in the Shiraiwa and Zhang model runs. Species that did not require updates to C⁰ still had wide ranges of predicted T_g between parameterizations, particularly species that are leading contributors to OA mass (AMT1, AMT2, ASVOO1, and ALVOO2) (Fig. S6).

The methods used to derive T_g parameterizations also had an impact on predicted total $T_{\mathrm{g},\mathrm{org}}:T$. When w_s≤0.1, the Li and Zhang model simulations almost never predicted a liquid organic coating (Fig. 3d–f), and rather largely predicted a semi-solid phase state in the Eastern U.S. and a solid phase state in the Western U.S. When compared with derived simulations of $T_{\mathrm{g},\mathrm{org}}:T$ (Fig. 5), predictions from the Zhang and Li parameterizations were also more viscous than that predicted in the Shiraiwa parameterization. In comparison to the Shiraiwa T_g equation, the Zhang and Li T_g parameterizations were fit with the inclusion of much lower-volatility OA, although not biased in that direction, and therefore may lead to a more informed estimate of phase state. In fitting T_g equations, the derivations in C⁰ between the Zhang and Li parameterizations differed with the potential to cause differences in phase state predictions (Li et al., 2020; Zhang et al., 2019b). In the Zhang parameterization, 13 OA T_g values were measured (Zhang et al., 2019b), and the p₀ values were either sourced from publications or via structural-activity relationships using the EVAPORATION model (Compernolle et al., 2011). In the Li parameterization, a similar process was followed in first searching the literature for established p₀, however, due to the lack of observed T_g values of some sulfur-containing compounds, T_g values for CHOS compounds (∼ 35 % compounds in the training dataset) were estimated based on the Boyer-Kauzmann rule (Li et al., 2020). With using the Boyer-Kauzmann rule, melting temperatures were estimated using EPI Suite along with p₀s also estimated using EPI Suite (Li et al., 2020; USEPA, 2012). In comparison to other structural-activity relationship models, EVAPORATION predicts higher volatilities than both OPERA and EPI Suite (Mansouri et al., 2018). In addition, EPI Suite can predict higher melting temperatures than OPERA (Mansouri et al., 2018). Associating a T_g with a higher volatility has the potential to lower T_g and to predict a less viscous phase state, and associating a higher melting temperature with T_g has the potential to predict a more viscous phase state. Ultimately, model-measurement $T_{\mathrm{g},\mathrm{org}}:T$ performance is underpowered in this study and warrants the need for future field measurements of aerosol phase state.

4.2 Impacts on IEPOX reactive uptake

Across all simulations that considered the effect of phase separation and diffusivity limitations, γ_IEPOX was reduced in comparison to the base simulation and to different degrees based on the phase state of the organic coating predicted (Fig. 6). These reductions resulted in furthering the negative bias in predicting total IEPOX-SOA and IEPOX organosulfates (Fig. 7a and c); however, it improved the positive bias in predicting 2-methyltetrols (Fig. 7b). While the impacts of phase state on D_org,eff and subsequently γ_IEPOX were explored in this study, there are other parameters that γ_IEPOX has been known to be more sensitive to that warrant further exploration. Without phase separation and phase state implemented, when increasing the k_particle from 0.0001 to 0.001 s⁻¹ (consistent with bounds seen in Budisulistiorini et al. (2017), γ_IEPOX increases by one order of magnitude. The AIEOS specific third-order rate constant used here (Table S3) and in previous studies (Budisulistiorini et al., 2017; Riedel et al., 2016) is slower than that determined by Piletic et al. (2013) (based on AIEOS:AIETET from laboratory studies) (Budisulistiorini et al., 2017; Riedel et al., 2016), and may be partially responsible for the underpredictions in AIEOS and total IEPOX-SOA. As the AIEOS:AIETET ratio is considered in determining the third-order rate constant for AIEOS (Budisulistiorini et al., 2017), it has been suggested that this rate constant should be higher (Chen et al., 2024), which could work to increase γ_IEPOX to a level consistent with chamber studies (D'Ambro et al., 2019). The total k_particle for the formation of AIEOS also depends on particle acidity (Cooke et al., 2024; Gaston et al., 2014; Riedel et al., 2016; Riva et al., 2016) and particulate sulfate concentrations (Cooke et al., 2024; Jo et al., 2021; Riva et al., 2019) (Eq. S5, Table S3) (Eddingsaas et al., 2010), and an underprediction in AIEOS and subsequently IEPOX-SOA, can possibly be attributed to model underpredictions of sulfate (Fig. S7). The addition of missing heterogeneous sulfate formation pathways increased the mean bias in sulfate concentration predictions at the LRK site by ∼ 0.2 µg m⁻³ for July 2016 (Farrell et al., 2025) and could potentially help resolve model-measurement differences for AIEOS. A combination of increased particle-phase reaction rates due to increased sulfate and volatilization of 2-methytetrols could improve both the ratio of predicted 2-methyltetrols to organosulfates as well as their overall abundance in model predictions. A limitation of this study is that water did not move between the organic coating and the inorganic aqueous core (Schmedding and Zuend, 2023), which may pose another impact on the k_particle – increasing it or decreasing it based on dilution or concentration of acids and nucleophiles (Cooke et al., 2024; Gaston et al., 2014; Riva et al., 2016).

4.3 Atmospheric Implications

In this study we explore the impacts of aerosol organic-inorganic phase separation and outer organic coating phase state on IEPOX heterogeneous reactive uptake, however, aerosol phase separation and phase state could potentially impact other heterogeneous aerosol formation pathways including secondary sulfur and secondary organic aerosol formation (Farrell et al., 2025; Huang et al., 2023). In more humid regions of the U.S. at the surface and at altitudes where liquid clouds occur, the contribution of water to the outer organic shell resulted in its phase state being liquid, which had relatively little impact on IEPOX diffusivity into the inorganic aqueous core. In drier regions of the U.S. and at higher altitudes, the outer organic shell's more solid phase state could pose more of a diffusional barrier to IEPOX reactive uptake, however IEPOX concentrations in these locations are not as abundant.

The effect of phase separation has been experimentally quantified for liquid cloud droplet activation (Davies et al., 2019; Ovadnevaite et al., 2017; Schmedding and Zuend, 2023; Vepsäläinen et al., 2022). Accounting for phase separation in predicting the hygroscopic growth of cloud condensation nuclei (CCN) resulted in an increase in the super-saturation required for CCN activation (Davies et al., 2019); however, the presence of organics and liquid-liquid phase separation has been seen to lower the surface tension of liquid particles and reduces underpredictions in CCN (Ovadnevaite et al., 2017). Further theoretical work has found that liquid-liquid phase separation with a thin organic film can result in substantial CCN activation (Schmedding and Zuend, 2023; Vepsäläinen et al., 2022) for smaller particles (Schmedding and Zuend, 2023). All of the above point to the impacts of liquid-liquid phase separation on low-lying cloud formation and opacity, which are also recent avenues for model improvement in predicting the radiative budget (Szopa et al., 2023). In addition, phase separation has recently been shown to also impact the radiative properties of some aerosols (Fard et al., 2018; Zhang et al., 2022). Fard et al. (2018) found that liquid-liquid phase separation in brown carbon aerosols increased their scattering cross-sections and decreased their absorbing cross-sections. Zhang et al. (2022) found that with increasing organic coatings in a phase separated aerosol, that black carbon absorption decreased. Accounting for phase separation can impact the potency of aerosols as short-lived climate forcers (Fard et al., 2018; Zhang et al., 2022) have the potential to influence meteorological predictions, and consequently pollutant concentrations in meteorological-chemical coupled transport models (Forkel et al., 2015; Gao et al., 2024; Hogrefe et al., 2015; Wang et al., 2021; Wong et al., 2012).

The phase state of phase separated aerosols can impact the formation of cirrus clouds (Berkemeier et al., 2014; Cziczo et al., 2013; Maclean et al., 2021; Murray et al., 2010; Wagner et al., 2012; Wolf et al., 2020). For all sensitivity simulations, in higher layers of the modeling grid (> 8 km; i.e., model vertical layers exceeding 28th layer), the phase state of the outer organic coating was in either a semi-solid or solid phase state attributed to modeled RH (less than 60 %) and modeled T (less than 262 K) (Fig. 4). The phase state (influenced by RH and temperature) of phase-separated aerosols have been shown to influence both liquid water uptake to solidified inorganic cores (immersion freezing) (Berkemeier et al., 2014) and the formation of ice from water vapor on the outer organic coating (heterogeneous depositional freezing) (Murray et al., 2010). Heterogeneous depositional freezing would be a particularly important pathway for cirrus cloud formation in the context of furthering this work, as the mixing time of water and organic species at higher altitudes – where air-mass updrafts increase cooling rates and RH is lower – has been shown to be slow (Maclean et al., 2021).

Accounting for phase state and phase separation also has the potential to extend the atmospheric life-times (i.e., semi-solid or solid phase state) and thus transport of harmful air pollutants (HAPs) (Mu et al., 2018; Shrivastava et al., 2017). Shrivastava et al. (2017) found that accounting for aged organic coatings resolved model-measurements gaps of Benzo[a]pyrene (BaP), a HAP with lung carcinogenicity (Boström et al., 2002; Bukowska et al., 2022), suggesting shielding from the organic coating limits photo-oxidant degradation. Mu et al. (2018) found that the lifetime of BaP (against ozone) in phase-separated aerosols decreases with increased RH. Both of these studies highlight the importance of long-range transport of HAPs that can harm human health particularly in the mid-to-upper latitudes (Mu et al., 2018; Shrivastava et al., 2017) and highlight the importance of accurately predicting phase separation and phase state in air quality models.

5 Conclusions

Although inorganic-organic phase separation and phase state of aerosols, is known to happen (Riva et al., 2019; Zhang et al., 2019a), these physical properties are not traditionally included in chemical transport models, yet can impact heterogeneous aerosol formation. In this work, we demonstrate that the phase state of phase separated aerosols can decrease the heterogeneous reactive uptake of IEPOX and subsequent formation of IEPOX-SOA. This decrease is largely attributed to a more solid phase state of the outer organic shell (Figs. 3a–c, 6). While this reduction overall furthers a negative bias in simulating total IEPOX-SOA, and explicitly tracked methyltetrol sulfates (AIEOS), it reduces a positive bias in simulating explicitly tracked 2-methyltetrols (AIETET). We find that with increased ALW associated with the outer organic shell, its phase state is more liquid-like and does not pose as much of a diffusional resistance to heterogeneous reactive uptake. The simulated phase states (also modulated by T_g parameterization used), can be used as bounds for SAPRC07tic_ae7i organic aerosol species in future studies. Within each phase state parameterization, future modeling studies should also aim to update model organic aerosol properties (C⁰, M, and O:C) upon which phase state parameterizations rely. Ultimately phase separation and phase state of aerosols have implications beyond impacting heterogeneous chemistry and more field measurements of aerosol phase state are warranted to further validate the organic aerosol phase states simulated in this study.