Predicted impacts of heterogeneous chemical pathways on particulate sulfur over Fairbanks (Alaska), the Northern Hemisphere, and the Contiguous United States

A portion of Alaska's Fairbanks North Star Borough was designated as nonattainment for the 2006 24 h fine particulate matter 2.5 µm or less in diameter (PM_2.5) National Ambient Air Quality Standards (NAAQS) in 2009. PM_2.5 NAAQS exceedances in Fairbanks mainly occur during dark and cold winters, when temperature inversions form and trap high emissions at the surface. Sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), often the second-largest contributor to PM_2.5 mass during these wintertime PM episodes, is underpredicted by atmospheric chemical transport models (CTMs). Most CTMs account for primary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formed via gas-phase oxidation of sulfur dioxide (SO₂) and in-cloud aqueous oxidation of dissolved S(IV). Dissolution and reaction of SO₂ in aqueous aerosols are generally not included in CTMs but can be represented as heterogeneous reactive uptake and may help better represent the high ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations observed during Fairbanks winters. In addition, hydroxymethanesulfonate (HMS), a particulate sulfur species sometimes misidentified as ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, is known to form during Fairbanks winters. Heterogeneous formation of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS in aerosol liquid water (ALW) was implemented in the Community Multiscale Air Quality (CMAQ) modeling system. CMAQ simulations were performed for wintertime PM episodes in Fairbanks (2008) as well as over the Northern Hemisphere and Contiguous United States (CONUS) for 2015–2016. The added heterogeneous sulfur chemistry reduced model mean sulfate bias by ∼ 0.6 µg m⁻³ during a cold winter PM episode in Fairbanks, AK. Improvements in model performance are also seen in Beijing during wintertime haze events (reducing model mean sulfate bias by ∼ 2.9 µg S m⁻³). This additional sulfur chemistry also improves modeled summertime ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ bias in the southeastern US, with implications for future modeling of biogenic organosulfates.

1 Introduction

Radiative forcing and climate effects attributed to fine particulate matter 2.5 µm or less in diameter (PM_2.5) remain among the most uncertain in climate change assessments (IPCC, 2013). Acute and long-term exposure to PM_2.5 has been associated with negative health outcomes including but not limited to acute myocardial infarction, stroke, and respiratory complications (Chen et al., 2018; Hayes et al., 2020; Silva et al., 2016; USEPA, 2019; Yang et al., 2021; Yitshak-Sade et al., 2018). Some of the most extreme PM_2.5 episodes in history have occurred when inverted vertical temperature profiles cause stable atmospheric conditions that limit pollution dilution (Holzworth, 1972; Malek et al., 2006; Scott, 1953; Wallace and Kanaroglou, 2009). Near-surface temperature inversions (with a minimum modeled planetary boundary layer height of 1.7 m) are characteristic of Fairbanks and North Pole (Alaska) winters and are associated with degraded air quality in this region (Malingowski et al., 2014; Mayfield and Fochesatto, 2013). Wintertime PM_2.5 episodes have impacted human health in these cities (McLaughlin and Castrodale, 2010), with this region exceeding the National Ambient Air Quality Standards (NAAQS) for PM_2.5 since 2009, when portions of the Fairbanks North Star Borough were designated as nonattainment for the 2006 24 h PM_2.5 NAAQS¹ (ADEC, 2017, 2019).

Sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), often a major component of PM_2.5 in Fairbanks and North Pole (ADEC, 2017) as well as globally (Snider et al., 2016), can be emitted directly (primary) or formed secondarily via atmospheric oxidation of sulfur dioxide (SO₂). Known secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation processes include but are not limited to gas-phase oxidation of SO₂ (Calvert et al., 1978), particle surface oxidation of SO₂ (Clements et al., 2013; Wang et al., 2021), and in-cloud aqueous-phase oxidation of inorganic sulfur species with oxidation number 4 (S(IV) = SO₂ ⋅ H₂O + ${\mathrm{HSO}}_{\mathrm{3}}^{-}+{\mathrm{SO}}_{\mathrm{3}}^{\mathrm{2}-}$) (secondary) (Hoffmann and Calvert, 1985; Ibusuki and Takeuchi, 1987; Lagrange et al., 1994; Lee and Schwartz, 1983a; Maahs, 1983; Maaß et al., 1999; Martin and Good, 1991; McArdle and Hoffmann, 1983). Aside from contributing directly to PM_2.5 mass, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ can facilitate the formation of other PM_2.5 species as a reactant (Brüggemann et al., 2020; Huang et al., 2019; Huang et al., 2020; Surratt et al., 2010) by increasing aerosol water uptake (Kim et al., 1994; Nguyen et al., 2014) and by altering aerosol acidity (Li et al., 2022; Pye et al., 2020).

Under heavily polluted haze conditions, such as those common in the North China Plain during the winter, recent studies have suggested that secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ may be efficiently produced in aerosol liquid water (ALW) (Cheng et al., 2016; Fan et al., 2020; Liu et al., 2020). Hygroscopic PM_2.5 (both inorganic and organic) can increase ALW content (Nguyen et al., 2014; Petters and Kreidenweis, 2007; Pye et al., 2017), which can facilitate secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation (Zhang et al., 2021b), enhancing ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in a positive feedback loop – which is particularly important during high-relative-humidity haze events (Cheng et al., 2016; Song et al., 2021a; Wang et al., 2016, 2014).

Chemical transport models generally include secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation via gas-phase oxidation of SO₂ by OH and in-cloud aqueous-phase oxidation of dissolved S(IV) species by oxidants such as hydrogen peroxide (H₂O₂); ozone (O₃); peroxyacetic acid (PAA); methyl hydroperoxide (MHP); and/or oxygen (O₂) catalyzed by transition metal ions (TMIs-O₂), iron (Fe³⁺), and manganese (Mn²⁺). Limited to these formation pathways, CTMs have been unable to reproduce the high levels of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ observed during wintertime PM events in Beijing and regions that experience extremely cold winters. This persistent underprediction suggests that CTMs lack ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation pathways that take place when photochemistry and cloud liquid water are limited (Eckhardt et al., 2015; Gao et al., 2016; Wang et al., 2014). Previous studies have suggested that heterogeneous sulfate formation in deliquesced aerosol particles may account for at least part of this underprediction (Wang et al., 2014; Zheng et al., 2015). Wang et al. (2014) implemented generalized heterogeneous reactive uptake of SO₂ to form ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (where relative-humidity-dependent uptake coefficients were specified rather than calculated) in the GEOS-Chem global model (Chen et al., 2009), and this led to improved model–observation comparisons for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ during a wintertime haze event over North China (Wang et al., 2014). Similarly, when Zheng et al. (2015) implemented generalized heterogeneous sulfur chemistry in WRF-CMAQ, the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ normalized mean bias decreased from −54.2 % to 6.3 % for Beijing haze events (Zheng et al., 2015).

The inclusion of hydroxymethanesulfonate (HMS) chemistry in CTMs may also ameliorate negative model bias. General ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ PM measurement methods struggle to disentangle ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ spectra from those of inorganic S(IV) and HMS (Dovrou et al., 2019). HMS is formed from the aqueous-phase reaction of S(IV) with formaldehyde (HCHO) (Boyce and Hoffmann, 1984; Deister et al., 1986; Kok et al., 1986; Kovacs et al., 2005; Olson and Hoffmann, 1986). Recent field, modeling, and experimental studies have highlighted the importance of secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS formation in aqueous aerosols during wintertime haze events (Campbell et al., 2022; Dovrou et al., 2019; Ma et al., 2020; Moch et al., 2020; Song et al., 2019). Moch et al. (2020) suggested up to 25 % of measured ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ may actually be HMS in heavily polluted regions. Arctic, sub-Arctic, and regions that experience extremely cold and dark winters may favor HMS formation, as colder temperatures increase the solubility of SO₂ and HCHO, and limited sunlight reduces the photo-oxidation of HCHO (Moch et al., 2020; Pandis and Seinfeld, 1989; Sander, 2015; Song et al., 2019; Staudinger and Roberts, 1996). In the work of Song et al. (2021a), including generalized heterogeneous cloud HMS chemistry in GEOS-Chem reduced model–measurement differences in $\mathrm{HMS}/{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ ratios. The inclusion of heterogeneous formation and loss of HMS in deliquesced aerosols increased modeled HMS concentrations, particularly in China (Song et al., 2021b).

Exploring heterogeneous sulfur chemistry in deliquesced aerosols requires retrofitting of multi-phase reactions. Many laboratory studies characterizing rate coefficients and expressions for aqueous-phase ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation have been performed under dilute conditions, characteristic of cloud droplets. The ionic strength of aerosol particles, however, can be several orders of magnitude higher than that of cloud droplets due to aerosols containing much lower concentrations of water (Mekic and Gligorovski, 2021). Ionic strength has been found to impact the aqueous-phase formation of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, increasing the rate of aqueous-phase kinetics for NO₂ and H₂O₂ oxidation of S(IV) (Chen et al., 2019; Liu et al., 2020) and inhibiting the aqueous-phase kinetics of TMI-catalyzed O₂ oxidation (Ibusuki and Takeuchi, 1987; Martin and Good, 1991; Martin and Hill, 1987). Experimental studies have also shown that high ionic strength may increase or decrease (effective) Henry's law coefficients of reactants compared to pure water (Ali et al., 2014; Chen et al., 2019; Kosak-Channing and Helz, 1983; Lagrange et al., 1994; Liu et al., 2020; Millero et al., 1989; Shao et al., 2019).

In this paper we describe the implementation of heterogenous sulfur chemistry in ALW in the Community Multiscale Air Quality (CMAQ) v5.3.2 modeling system (USEPA, 2020), leading to additional ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS formation. We refer to this chemistry as heterogeneous given the use of the heterogeneous framework (Hanson et al., 1994). Heterogeneous sulfur chemistry pathways implemented include the oxidation of dissolved S(IV) species by H₂O₂, O₃, PAA, MHP, TMI-O₂, and NO₂ and the in-aerosol aqueous formation of HMS. In addition to heterogeneous chemistry updates, ionic strength effects were added to condensed-phase rate expressions and Henry's law coefficients of some species. The updated model was applied for several time periods and for different domains and horizonal resolutions. Two historical wintertime PM episodes were simulated for a finely resolved (1.33 km) domain centered over Fairbanks, Alaska; winter and summer periods over the Contiguous United States (CONUS) (12 km) during 2016; and the 2015–2016 winter season over the Northern Hemisphere (108 km) to investigate the impacts of these updates for different chemical regimes, domains, and seasons. Changes to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, HMS, and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+$ HMS (PM_2.5,sulf) predictions were tracked with each update (i.e., for (1) adding heterogeneous sulfur reactions and (2) adding ionic strength effects), and model performance was evaluated with available observations. This study aims to better understand the impacts that heterogeneous sulfur chemistry parameterizations may have on predicted PM_2.5,sulf concentrations and whether the additional chemistry can resolve ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ underpredictions in cold and dark conditions.

2 Methods

2.1 Heterogeneous sulfur chemistry

In this study, reactions that transform SO₂ to SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS are simulated in both clouds and ALW. In ALW, the production of particulate ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS is parameterized as a set of first-order heterogeneous reactions of SO₂ (gas) or reactants (e.g., HCHO, O₃, NO₂) as follows:

$$\begin{array}{}\text{(R1)}& P_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{reactant}}\stackrel{k_{\mathrm{het}}}{⟶}{\mathrm{PM}}_{\mathrm{2.5},\mathrm{sulf}},\end{array}$$

where PM_2.5,sulf refers to both sulfate and HMS, and k_het is the heterogeneous rate constant (Eq. 1), which accounts for gas-to-particle mass transfer processes, aqueous reactive uptake, and sulfur transformations (Jacob, 2000).

$$\begin{array}{}\text{(1)}& k_{\mathrm{het}}\left({\mathrm{s}}^{-\mathrm{1}}\right)={\displaystyle \frac{\mathrm{SA}}{\frac{r_{\mathrm{p}}}{D_{\mathrm{g}}}+\frac{\mathrm{4}}{v\mathit{\gamma }}}}\end{array}$$

Here, SA is the aerosol surface area (m² m⁻³), r_p is the effective particle radius (m), D_g is gas-phase diffusivity of a reactant (m² s⁻¹), and v is the mean molecular speed of the partitioning gas (m s⁻¹). γ is the reactive uptake coefficient and is given by the following equation (Hanson et al., 1994; Jacob, 2000; Schwartz, 1986):

$$\begin{array}{}\text{(2)}& \mathit{\gamma }=[{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}+{\displaystyle \frac{v}{\mathrm{4}HRT\sqrt{D_{\mathrm{a}}{k}_{\mathrm{chem}}}}}\cdot {\displaystyle \frac{\mathrm{1}}{\coth \left(q\right)-\frac{\mathrm{1}}{q}}}{]}^{-\mathrm{1}},\end{array}$$

where α is the mass accommodation coefficient of a species, H is the effective Henry's law coefficient (M atm⁻¹) at temperature T (K), R is the ideal gas constant 0.08206 (L atm mol⁻¹ K⁻¹), D_a is the aqueous-phase diffusivity (assumed here to be 10⁻⁹ m² s⁻¹), k_chem represents the pseudo-first-order condensed-phase rate coefficient (s⁻¹) (Table 1), and q is the diffuso-reactive parameter defined as $q=r_{\mathrm{p}}\sqrt{\frac{k_{\mathrm{chem}}}{D_{\mathrm{a}}}}$ (Schwartz and Freiberg, 1981). The rate expressions in aerosol water for the base configuration (hereafter referred to as “Base_Het”) are the same as those used for cloud chemistry (Appel et al., 2017; Sarwar et al., 2013). Heterogeneous sulfur chemistry was calculated when relative humidity (RH) was greater than or equal to 50 % following Shao et al. (2019), assuming that below a 50 % RH, the aerosol water content would be too low for heterogeneous sulfur chemistry to take place (Sun et al., 2013). The heterogeneous chemistry is also solved simultaneously with gas-phase chemistry.

Table 1 Pseudo-first-order rate constants (k_chem) and ionic strength (I) factors for each simulation.

^a The ionic strength cap chosen corresponds with the activation energy for the rate constant, k₁, which matches the ionic strength cap, b, for sodium perchlorate (NaClO₄). ^b [S(IV)](M)= [SO₂⋅H₂O] + [HSO${}_{\mathrm{3}}^{-}$] + [SO${}_{\mathrm{3}}^{\mathrm{2}-}$] ^c The b value corresponds to the factor associated with Mn-catalyzed O₂ oxidation of SO₂, corresponding to ≥ 10⁻⁴ M sulfur (Martin and Hill, 1987). ^d Also included is the loss of HMS from decomposition and OH oxidation (Song et al., 2021a).

2.2 Accounting for ionic strength effects and alternative chemical rate expressions

Several studies have investigated the impact of high ionic strength on sulfate oxidation rates originally developed for dilute conditions. Depending on the reaction pathway, rates may be enhanced or diminished with increased ionic strength (Chen et al., 2019; Lagrange et al., 1994; Liu et al., 2020; Shao et al., 2019). Ionic strength is calculated in CMAQ as

$$\begin{array}{}\text{(3)}& I={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\sum }_{i=\mathrm{1}}^n{m}_i{z}_i^{\mathrm{2}},\end{array}$$

where m_i is the solute concentration (M) in aerosol or cloud liquid water, and z_i is the charge associated with each modeled ion.

SO₂ can be oxidized in the aqueous phase by O₂ when catalyzed by TMI (specifically Fe³⁺ and Mn²⁺) with synergies existing when both Fe³⁺ and Mn²⁺ are present (Altwicker and Nass, 1983; Ibusuki and Takeuchi, 1987; Martin and Good, 1991). This reaction pathway has been found to be an important secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation pathway, especially at low pH and when photochemistry is limited and Fe is more soluble (Cheng et al., 2016; Li et al., 2020; Liu et al., 2022; Song et al., 2021a). Ibusuki and Takeuchi (1987) found that the rate of S(IV) oxidation by the TMI-O₂ pathway peaked around pH = 4.2, decreased with decreased temperature, and was enhanced by high concentrations of TMI (Ibusuki and Takeuchi, 1987). Martin and Good (1991) investigated the impact of higher S(IV) concentrations on this ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation pathway and found that higher S(IV) concentrations alter the rates of catalysis via Fe³⁺ and Mn²⁺, along with their synergistic catalysis, but did not explore a pH or temperature dependency. However, at similar pH ranges, S(IV) concentrations, and soluble Fe³⁺ and Mn²⁺ concentrations, Martin and Good (1991) found that their rate expression agreed well with that by Ibusuki and Takeuchi (1987). When Shao et al. (2019) implemented the Ibusuki–Takeuchi TMI-catalyzed O₂ oxidation rate expression in ALW in GEOS-Chem, they found that this pathway accounted for 67 %–69 % of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formed over China (Shao et al., 2019). The presence of higher S(IV) concentrations, however, may warrant the use of a ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation rate that takes into consideration faster rates of TMI catalysis (Martin and Good, 1991). The TMI-O₂ oxidation pathway from Martin and Good (1991) is used in the Base_Het simulation, and the TMI-O₂ oxidation pathway from Ibusuki and Takeuchi (1987) is used in the TMI_sens simulation to explore the range in ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation possible by this pathway (Table 1). For both implementations of this pathway, ionic strength impacts were added as high ionic strength has been found to limit the TMI-catalyzed oxidation of S(IV) to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (Martin and Hill, 1967, 1987).

While the TMI-catalyzed oxidation pathway can be a significant contributor to secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation, especially at low pH, some studies suggest ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation during Beijing winter haze events may be dominated by the reaction of SO₂ and NO₂ in aerosol and/or cloud and fog water under mildly acidic or neutral conditions (Cheng et al., 2016; Wang et al., 2020; Yang et al., 2019). This oxidation pathway may increase in importance with increasing ionic strength (Cheng et al., 2016). Recent chamber work by Chen et al. (2019) found that increasing ionic strength increased the rate of secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation from NO₂ oxidation. To assess the potential impact ionic strength may have on this pathway, the ionic-strength-dependent NO₂ oxidation rate of Chen et al. (2019) was included in the TMI_NO₂ sensitivity simulation (Table 1).

While the aforementioned pathways are pH dependent, SO₂ aqueous oxidation by H₂O₂ is pH independent for pH > 2 due to the opposing dependencies of the reaction rate coefficient and S(IV) solubility on pH (Clifton et al., 1988; Seinfeld and Pandis, 2016). This pathway has been studied extensively in dilute conditions representative of cloud droplets (Hoffmann and Calvert, 1985; Maaß et al., 1999; McArdle and Hoffmann, 1983). Maaß et al. (1999) found that oxidation of S(IV) by H₂O₂ increases with increased ionic strength and formulated a semi-empirical relationship between ionic strength and the reaction rate coefficient for the S(IV)–H₂O₂ oxidation pathway (the upper limit of ionic strength in this study is 5 M). Field measurements during Chinese haze events have reported ionic strengths of aerosols ranging from 14–43 M (Cheng et al., 2016; Fountoukis and Nenes, 2007). Liu et al. (2020) recently studied this pathway using aerosol flow reactors to determine an ionic-strength-based enhancement factor that encapsulates the combined ionic strength effects on Henry's law coefficients, dissociation, and condensed-phase kinetics. This study found that increasing ionic strength from 0 to 14 M resulted in an order of magnitude increase in ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ production rate (Liu et al., 2020).

O₃ is another important aqueous-phase oxidant of SO₂ with a reaction rate that increases with increasing pH (Maahs, 1983) and ionic strength (Lagrange et al., 1994). The ionic strength enhancement factor for this rate has been implemented with other ionic strength enhancement or inhibition factors in the All_Ionic sensitivity simulation (Table 1).

Table 2 Henry's law coefficients used per model simulation.

^* Aqueous-phase concentrations of S(IV) are calculated similarly to Base_Het, TMI_sens, and TMI_NO₂_sens but with ionic-strength-dependent equilibrium coefficients.

For the fine-mode aerosol, total gas and particle concentrations of key inorganic species are passed to the thermodynamic equilibrium model ISORROPIA II to calculate aerosol pH and ALW (Fountoukis and Nenes, 2007). In addition to the ALW associated with inorganic ions, ALW associated with organic aerosols is also estimated in CMAQ via hygroscopicity parameters (Pye et al., 2017).

2.3 Model base case and sensitivity simulations

Several CMAQ configurations were used here to understand the impacts of using heterogeneous sulfur chemistry, ionic strength, and alternative pseudo-first-order rate expressions. A base case CMAQ simulation (“base”) was completed using in-cloud ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation from aqueous oxidation by H₂O₂, O₃, PAA, and MHP and via TMI-catalyzed O₂ of SO₂ and gas-phase oxidation of SO₂ by OH (Fahey et al., 2017; Sarwar et al., 2013).

To account for the impacts of heterogeneous sulfur chemistry in ALW, the Base_Het model simulation was completed for all domains (see “Model configuration”) using the aforementioned heterogenous reactive uptake parameterizations (Eqs. 1–2). Parameters from CMAQ's KMT2 cloud chemistry model (Fahey and Roselle, 2019; Fahey et al., 2017) were used to calculate k_chem (Table 1). An ionic strength inhibition term was added to the TMI-catalyzed O₂ pathway in aerosol water to account for the limiting effect of ionic strength on this pathway. The ionic strength was capped at I_max= 2 M to reflect experimental constraints (Martin and Good, 1991; Martin and Hill, 1987; Seinfeld and Pandis, 2016). For this pathway in both clouds and ALW, Fe³⁺ was assumed to be 90 % of dissolved Fe at night and 10 % during daytime, with soluble fractions of Mn and Fe assumed to be 0.5 and 0.1 of total Fe and Mn respectively (Alexander et al., 2009). Other sources of less soluble Fe and Mn emissions, such as dust, are likely minimal given the snow cover for this domain and episode (Shao et al., 2019).

The k_chem for the TMI-catalyzed O₂ pathway in the Base_Het case is neither temperature nor pH dependent (Martin and Good, 1991; Martin and Hill, 1987). Ibusuki and Takeuchi (1987) found both a pH and a temperature dependence on k_chem for this pathway (Martin and Hill, 1987). To explore the effects of both a pH and a temperature dependence on the rate of a TMI-catalyzed sulfur oxidation pathway, a sensitivity simulation, TMI_sens, was run (Ibusuki and Takeuchi, 1987). Given that this k_chem is reduced in colder temperatures, the TMI_sens run likely represents a lower bound on ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation for this pathway during winter episodes. This k_chem also uses the same solubility, dissociation, and ionic strength bounds as the TMI-catalyzed O₂ oxidation pathway used in the Base_Het simulation.

In the TMI_NO₂_sens simulations, both the alternative k_chem for the TMI-catalyzed O₂ pathway and an ionic-strength-dependent k_chem for the NO₂ oxidation pathway in ALW are included (Chen et al., 2019; Ibusuki and Takeuchi, 1987). Both the TMI-O₂ and the NO₂ oxidation k_chem values favor weakly acidic pH regimes (Cheng et al., 2016; Martin and Good, 1991). This sensitivity simulation was implemented to analyze potential competition between two pathways under their favorable pH conditions and conditions that are also characteristic of wintertime haze episodes (weakly acidic) (Cheng et al., 2016).

Ionic strength impacts on k_chem for H₂O₂ and O₃ aqueous oxidation formation pathways were included on top of the previous modifications in the TMI_sens and TMI_NO₂_sens model simulations (Chen et al., 2019; Ibusuki and Takeuchi, 1987; Maaß et al., 1999). Ionic strength adjustments were also included for S(IV) dissociation constants and Henry's law coefficients for H₂O₂, O₃, and SO₂ in the All_Ionic simulations (Ali et al., 2014; Lagrange et al., 1994; Maaß et al., 1999; Millero et al., 1989; Seinfeld and Pandis, 2016) to analyze the combined effects of ionic strength on total modeled PM_2.5,sulf aerosol.

In addition to the reactions shown in Table 1 (which are treated in both aerosol and cloud water), also included is S(IV) oxidation by peroxyacetic acid (PAA) and methyl hydroperoxide (MHP) in aerosol and cloud water and in-cloud S(IV) oxidation by HNO₄, OH, and NO₃ (Lee and Schwartz, 1983a; Lind et al., 1987; Martin, 1984). With the exception of the base case, which used CMAQ's default AQCHEM cloud chemistry scheme, all other simulations used the KMT2 cloud chemistry scheme. KMT2 includes additional inorganic and organic chemistry compared to AQCHEM, including several additional S(IV) oxidation reactions as well as HMS formation and loss.

2.4 Model configuration

Results of the base and sensitivity simulations were compared for three different spatial domains: Fairbanks, Alaska; the Contiguous US (CONUS); and the Northern Hemisphere. Model simulations over Fairbanks and North Pole, Alaska, span two wintertime PM episodes (Episode 1 (E1): 25 January–11 February 2008 and Episode 2 (E2): 4–17 November 2008) with 2 d of spin-up and a horizontal resolution of 1.33 km. Model simulations over the northern hemispheric domain were performed for the winter season from December 2015–February 2016 with 2 months of spin-up and a horizontal resolution of 108 km following a standard US EPA configuration described by Mathur et al. (2017) and Appel et al. (2021). Model simulations over the CONUS domain were run for the months of January and July 2016 with 10 d of spin up and at a horizontal resolution of 12 km.

The Weather Research and Forecasting (WRF) model (Skamarock et al., 2008) was used to develop meteorology on all three domains. The Fairbanks WRF case follows the original configuration and case study of Gaudet et al. (2010, 2012) for Fairbanks, AK. This older WRF simulation was updated from WRFv3.3 (Gaudet et al., 2010, 2012) to WRFv4.1.1. with all geophysical and meteorological inputs reprocessed for compatibility with the more recent version of WRF and using the Mellor–Yamada–Nakanishi–Niino 2.5-order (MYNN2.5) closure scheme (Nakanishi and Niino, 2009). Sensitivity testing found performance improvements when updating the planetary boundary layer model change from the Mellor–Yamada–Janjić (MYJ) scheme to MYNN2.5. Evaluations also showed that the WRFv4.1.1 configuration captured the extreme temperature variations in these cases well, with 2 m temperature root mean square error (RMSE) of the order of 2–3 K, which is within historical benchmarks for complex geographical areas that are more difficult to model (Kemball-Cook et al., 2005). For the CONUS, meteorological inputs were sourced from WRFv4.1.1, and CMAQ inputs were processed with the Meteorology-Chemistry Interface Processor (MCIP; Appel et al., 2017; Otte and Pleim, 2010) version 5.0.

The emissions inputs for the two Fairbanks wintertime PM episodes were based on inputs provided by the Alaska Department of Environmental Conservation (ADEC, 2014). These emission estimates were from the base year (2008) used for the Fairbanks PM_2.5 moderate state implementation plan (SIP) and represent the best available emission estimates for the two time periods. Table 5.6-3 from Sect. 5.06 of ADEC (2014) provides a summary of the methods and inputs used to develop this emission inventory. For this model setup, we used the same inventory inputs and scripts and only updated the speciation for CMAQv5.3.2. The emission inventories for the northern hemispheric domain follow the same setup as described in Appel et al. (2021), where anthropogenic emissions were sourced from the 2010 Hemispheric Transport of Air Pollution version 2 (Janssens-Maenhout et al., 2015), biogenic emissions were sourced from the Model of Emissions of Gases and Aerosols from Nature (Guenther et al., 2012), soil and lightning NO were sourced from the Global Emissions Initiative (http://www.geiacenter.org, last access: 3 May 2021), biomass burning emissions were sourced from the Fire Inventory from the National Center for Atmospheric Research (Wiedinmyer et al., 2011), and on-road and non-road emissions were developed using the Motor Vehicle Emission Simulator v2014a. The emission inventories for the CONUS domain were sourced from the 2016v7.2 (beta and Regional Haze) modeling platform (Appel et al., 2021).

Gas-phase chemistry was simulated using the CB6r3 mechanism (Luecken et al., 2019), and aerosol dynamics were simulated using the aero7 module. Boundary conditions for the Fairbanks domain were sourced from EPA's Air QUAlity TimE Series (EQUATES) project for 2008 (USEPA, 2021). The sulfur-tracking method (STM) (Fahey and Roselle, 2019) was extended to include the new heterogeneous sulfur chemical pathways in order to track the contributions of each chemical reaction, primary emissions, and initial and boundary conditions to modeled SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (Appel et al., 2021).

2.5 Sulfate and PM_2.5,sulf observations

The model predictions were evaluated against available observations. While most monitoring networks report measurements for PM_2.5 or PM₁ ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, recent studies have indicated hydroxymethanesulfonate may be included in those observations (Dovrou et al., 2019; Ma et al., 2020; Moch et al., 2018; Song et al., 2019). Based on these findings, we compare measured ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ to modeled PM_2.5,sulf:

$$\begin{array}{}\text{(4)}& {\mathrm{PM}}_{\mathrm{2.5},\mathrm{sulf}}\left({\displaystyle \frac{\mathrm{\mu }\mathrm{g}}{{\mathrm{m}}^{\mathrm{3}}}}\right)={\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+\mathrm{HMS}\times {\displaystyle \frac{{\mathrm{MW}}_{{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}}{{\mathrm{MW}}_{\mathrm{HMS}}}}.\end{array}$$

PM${}_{\mathrm{2.5}}{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ observations in Fairbanks during 2008 were obtained from ADEC's Thermo Electron Partisol 2000 single-channel instrument w/SCC monitor 020900010 (State Office Building in Fairbanks Alaska; 64.840672, −147.722461) (USEPA, 2024a). ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ observations for the 2016 hemispheric domain were sourced from the United States Environmental Protection Agency (US EPA) Air Quality System (AQS) monitoring network; the Canadian National Air Pollution Surveillance (NAPS) monitoring network; the European Monitoring and Evaluation Programme (EMEP) monitoring network; and one monitor at Tsinghua University in Beijing, China (ECCC, 2022; Ma et al., 2020; Zhang et al., 2021a; Tørseth et al., 2012; USEPA, 2022). Modeled PM_2.5,sulf concentrations and measured ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ from the AQS and NAPS networks were cast in units of micrograms sulfur per meter cubed (µg S m⁻³) to match the measurement units used in the EMEP and Tsinghua University ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ measurements.

3 Results

3.1 Modeled particulate sulfur enhancement during dark and cold PM episodes in Fairbanks and North Pole, AK

3.1.1 Aerosol sulfur enhancements during a wintertime haze event (Episode 1 (E1))

The base simulation average E1 sulfate concentrations around Fairbanks and North Pole, AK, are ∼ 2–3.5 mg m⁻³ (Fig. 1a and c). Compared to the base, the Base_Het simulation leads to increased PM_2.5,sulf predictions concentrated around the cities of Fairbanks and North Pole as well as the region south of the Tanana River (Fig. 1b, d). The additional heterogeneous chemistry in the Base_Het simulation contributes up to an additional 11 µg m⁻³ of maximum daily PM_2.5,sulf compared to the base simulation in the region south of the Tanana River. Maximum daily differences are defined as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathrm{Maximum}\mathrm{daily}\mathrm{differences}=\\ & max({\mathrm{Avg}}_{\mathrm{daily},\mathrm{2}}-{\mathrm{Avg}}_{\mathrm{daily},\mathrm{1}}).\end{array}\end{array}$$

Figure 1 Episode average sulfate (a) and HMS (c) concentrations in the base simulation along with daily max differences in sulfate (b) and HMS (d) concentrations between the Base_Het and base CMAQ simulations over Fairbanks and North Pole, AK, for Episode 1 (from 25 January to 11 February 2008). HMS formation was not included in base CMAQ (i.e., HMS = 0 in the base simulation). The domain size is 264.67 km by 264.67 km with a grid cell resolution of 1.33 km by 1.33 km.

Enhancements in PM_2.5,sulf concentrations for the Base_Het simulation are mainly driven by increases in ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations (increasing up to 10.9 µg m⁻³ for daily maximum differences across the entire domain), with smaller impacts from HMS (increasing up to ∼ 1 µg m⁻³ for daily maximum differences across the entire domain). HMS concentrations are enhanced more in North Pole than Fairbanks, coinciding with higher HCHO emissions (Fig. S1 in the Supplement) from residential wood combustion combined with high co-located SO₂ emissions (from home-heating oil) along with lower temperatures (ADEC, 2017).

Out of all of the secondary PM_2.5,sulf formation pathways that are enhanced during dark, cold conditions (TMI-catalyzed O₂, NO₂, and the formation of HMS), the leading secondary ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation pathway in Base_Het is the TMI-catalyzed O₂ oxidation pathway in ALW (Fig. 2). The first-order condensed-phase rate constant (k_chem) of this pathway is lower than that of k_chem for NO₂ by almost 2 orders of magnitude for average modeled conditions characteristic of Fairbanks and North Pole for E1 (aerosol pH = 3.83, $\left[\mathrm{Fe}\left(\mathrm{III}\right)\right]=\mathrm{0.24}\mathrm{M},\left[\mathrm{Mn}\left(\mathrm{II}\right)\right]=$ 0.002 M, [SO₂] = 20 ppb, [NO₂] = 20 ppb, [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] = 3 µg m⁻³, [ALW] = 6 µg m⁻³, and Temp = 243 K) (Fig. S2) and is ∼ 1 order of magnitude higher than that for HMS formation in ALW. Despite the NO₂ k_chem being higher, the TMI-catalyzed O₂ heterogeneous rate of sulfate formation is dependent upon SO₂ partitioning into the particle, as Fe and Mn are both aerosol species, and simulated dark conditions reduce the conversion of Fe³⁺ to Fe²⁺ from daytime photochemical reactions (Alexander et al., 2009; Rao and Collett, 1998; Shao et al., 2019). The effective Henry's law coefficient for SO₂ increases with pH, while Henry's law coefficient for NO₂ remains low across the pH spectrum. This and a higher mass accommodation coefficient (by ∼ 2 orders of magnitude) for SO₂ compared to NO₂ contribute to the TMI-catalyzed O₂ pathway outcompeting the NO₂ pathway for this model configuration. The TMI-catalyzed O₂ heterogeneous reactive uptake pathway also outcompetes the H₂O₂ and O₃ heterogeneous reactive uptake pathways due to low photochemical activity with the dark conditions of this domain and episode.

Figure 2 Particulate sulfur process/chemistry contributions and speciation (SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS) in downtown Fairbanks (a) located at the State Office Building (64.84° N, 147.72° W; grid cell 108, 93) and North Pole (b) located at 64.76° N, 147.34° W (grid cell 122, 86) for Episode 1 (E1) speciated by source and/or formation pathway. Secondary aqueous formation of PM_2.5,sulf is divided into two categories: ALW and cloud, where ALW pathways represent the heterogeneous sulfur chemistry added in this study.

The formation of HMS is higher in North Pole, which can be colder than Fairbanks by up to ∼ 5 °C. Higher modeled HCHO emissions in North Pole (Fig. S1) along with colder temperatures increase the partitioning of HCHO into existing particles. This effect and the increase in HMS formation are more pronounced in the TMI_sens simulation (Fig. S3a and b).

Compared to k_chem used in Base_Het, the alternative k_chem in TMI_sens for the TMI-catalyzed O₂ pathway is ∼ 2 orders of magnitude lower under average conditions characteristic of this episode (Fig. S2), with the extremely cold temperatures decreasing the TMI-catalyzed O₂ oxidation k_chem. This ultimately results in a slower conversion of S(IV) species to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and subsequent competition with HMS formation, which increases with colder temperatures (Fig. S4). Maximum daily average enhancements in modeled PM_2.5,sulf concentrations in the TMI_sens simulation reach up to 25 µg m⁻³ at a grid cell in North Pole and are mostly attributed to high HMS concentrations (Fig. S3a and b). Formation rates of HMS are also dependent on aerosol pH, which can be higher in North Pole than in Fairbanks (Fig. S4). Higher pH increases the conversion of ${\mathrm{HSO}}_{\mathrm{3}}^{-}$ to ${\mathrm{SO}}_{\mathrm{3}}^{\mathrm{2}-}$, and the rate constant for HCHO reaction with ${\mathrm{SO}}_{\mathrm{3}}^{\mathrm{2}-}$ is 5 orders of magnitude higher than that of the reaction with ${\mathrm{HSO}}_{\mathrm{3}}^{-}$ (Boyce and Hoffmann, 1984). Lower ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ production rates in TMI_sens and lower modeled aerosol acidity compared to Base_Het likely contribute to the higher HMS formation (and loss) rates seen in the TMI_sens simulation. It is important to note that while aerosol acidity is modified by aqueous-phase formation of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, it is not modified by the formation of HMS in CMAQ.

When implementing an ionic-strength-dependent NO₂ rate expression (Chen et al., 2019) on top of this alternative TMI-catalyzed O₂ rate expression in the TMI_NO₂_sens model simulation, the NO₂ oxidation pathway outcompeted the formation of HMS (Fig. S3). While this model simulation compared well with the Base_Het model simulation in daily averaged PM_2.5,sulf concentrations, it predicted significantly lower PM_2.5,sulf concentrations than TMI_sens in North Pole due to a reduction in HMS predicted. Compared to downtown Fairbanks, North Pole has lower NO₂ emissions, but the switch to the ionic-strength-dependent reaction rate still leads to higher ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ production in North Pole compared to TMI_sens.

The ALL_Ionic model simulation resulted in similar PM_2.5,sulf predictions and pathway contributions to the TMI_NO₂_sens simulation for this episode (Fig. S3). In both Fairbanks and North Pole, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation attributed to heterogeneous reactive uptake via the H₂O₂ oxidation pathway increased slightly. For ionic strength factor calculations, ionic strength is capped at the maximum ionic strength considered in the experiments the parameterizations are based on (∼ 5–6 M) (Ali et al., 2014; Millero et al., 1989). The modeled ionic strength of ALW is typically at or above the maximum experimental ionic strength when deriving the effective Henry's law coefficients, dissociation coefficients, and k_chem for the H₂O₂ oxidation pathway – leading to expected higher dissolution and kinetics for this pathway. HMS concentrations for this model simulation also increased slightly in comparison to TMI_NO₂_sens. Assuming excess of dark oxidant precursors (TMI, NO₂, and HCHO) and maximum ionic strength, the leading higher HMS production rate in ALL_Ionic is likely due to a slightly higher range of aerosol pH (∼ 2–5.5) compared to the Base_Het and TMI_NO₂_sens simulations (Fig. S3d).

3.1.2 Aerosol sulfur enhancements with liquid cloud events (Episode 2 (E2))

Sulfate and HMS are known to form efficiently in cloud and fog droplets (Altwicker and Nass, 1983; Boyce and Hoffmann, 1984; Calvert et al., 1978; Clifton et al., 1988; Ibusuki and Takeuchi, 1987; Lee and Schwartz, 1983a; Martin and Good, 1991; McArdle and Hoffmann, 1983). In E1, there was minimal cloud or fog liquid water simulated; however, during E2 (4–11 November 2008), there were some periods where cloud/fog chemistry impacts on PM_25,sulf formation were evident.

Compared to E1, PM_2.5,sulf concentration enhancements were lower overall during E2. Differences between Base_Het and base simulations, however, are appreciable during this episode, with PM_2.5,sulf increasing up to 4.6 µg m⁻³ across the entire domain (daily maximum difference) (Fig. 3). Enhancements in PM_2.5,sulf are mainly driven by increased ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation in and around Fairbanks and North Pole; however, simulated HMS concentrations reached up to 4.4 µg m⁻³ south of the Tanana River (daily maximum) for this episode. Note that the base simulation included some contributions from in-cloud S(IV) oxidation (i.e., five S(IV) oxidation reactions from CMAQ's default cloud chemistry mechanism (AQCHEM)), while Base_Het included the additional in-cloud chemical reactions from the KMT2 cloud chemistry module.

Figure 3 Episode average sulfate (a) and HMS (c) concentrations in the base simulation along with daily max differences in sulfate (b) and HMS (d) concentrations between the Base_Het and base CMAQ simulations over Fairbanks and North Pole, AK, for Episode 2 (from 4 to 17 November 2008). HMS formation was not included in base CMAQ (i.e., HMS = 0 in the base simulation). The domain size is 264.67 km by 264.67 km with a grid cell resolution of 1.33 km by 1.33 km.

In the Base_Het simulation of Episode 2, PM_2.5,sulf was formed in both ALW and cloud liquid water (Fig. 4). Similar to Episode 1, the leading secondary formation pathway in downtown Fairbanks was the TMI-catalyzed O₂ pathway in both ALW and cloud water. This formation pathway split between the ALW and cloud formation pathways when modeled fog water content was around 0.025–0.05 g m⁻³ for 5 and 7 November; however, it was completely overtaken by the cloud formation pathway on 15 November when fog water content was > 0.175 g m⁻³ (at both Fairbanks and North Pole), highlighting that when surface clouds or fog water is present, reactions in cloud water can compete with those in aerosol water. In North Pole, the leading ALW PM_2.5,sulf formation pathway was also the TMI-catalyzed O₂ oxidation pathway, with contributions from this pathway in cloud water as well. Sulfate formed via S(IV) oxidation by NO₂ in cloud water was slightly higher in North Pole than downtown Fairbanks for all cloud events, even though NO₂ emissions were higher in downtown Fairbanks (Fig. S5). This could be due to higher cloud/fog pH in North Pole compared to downtown Fairbanks, as this pathway is known to favor pH > 5 (Clifton et al., 1988; Lee and Schwartz, 1983a; Littlejohn et al., 1993; Sarwar et al., 2013). HMS formation in ALW is also present in North Pole on the last day of Episode 2, corresponding with temperatures around −28 °C and aerosol pH ∼ 5.

Figure 4 Particulate sulfur process/chemistry contributions and speciation (SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS) in downtown Fairbanks (a) located at the State Office Building (64.84° N, 147.72° W; grid cell 108, 93) and North Pole (b) located at 64.76° N, 147.34° W (grid cell 122, 86) for Episode 2 speciated by source and/or formation pathway. Secondary aqueous formation of PM_2.5,sulf is divided into two categories: ALW and cloud, where ALW pathways represent the heterogeneous sulfur chemistry added in this study.

HMS contribution to PM_2.5,sulf concentrations is further increased on 16 November in North Pole in the TMI_sens model simulation (Fig. S7b). Despite this increase in HMS, total PM_2.5,sulf predicted in this simulation is lower than in the Base_Het simulation for most of the episode. The compensation of the HMS formation pathway for the TMI-catalyzed O₂ formation pathway is not as significant for this episode despite similar HCHO emissions (Fig. S5) and is likely due to higher temperatures (∼ +15 °C warmer than E1).

3.1.3 Improved model performance for PM_2.5,sulf in Fairbanks

Daily average PM_2.5,sulf concentrations for the base and Base_Het simulations were compared with 24 h ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ measurements taken every third day at the State Office Building in downtown Fairbanks (Fig. 5) (USEPA, 2024a).

Figure 5 Time series comparing modeled PM_2.5,sulf concentrations to measured ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations at the State Office Building in downtown Fairbanks (64.84° N, 147.72° W; grid cell 108, 93) for Episode 1 (a) and Episode 2 (b). The lines represent daily average modeled values, and the shading represents hourly maximum and minimum concentrations for each day. The black lines represent the sampling period of the monitor measurements, with black circles representing the mid-point of the 24 h sampling period.

With the inclusion of heterogeneous sulfur chemistry in CMAQ (Base_Het), the mean bias in PM_2.5,sulf improved by ∼ 0.62 µg m⁻³ during E1 and by ∼ 0.36 µg m⁻³ during E2, reducing the model bias by up to ∼ 1 µg m⁻³ during E1 and up to ∼ 0.85 µg m⁻³ during E2. CMAQ still underpredicts PM_2.5,sulf observations for both episodes by ∼ 11 and 7 µg m⁻³, particularly during E1 for 6 and 9 February respectively. PM_2.5 speciation measurements during both episodes were only available for a single monitor located in downtown Fairbanks, and therefore we were not able to assess model performance in North Pole, where modeled PM_2.5,sulf concentrations can be higher by up to ∼ 5–27 µg m⁻³ (daily maximum difference) for E1 and ∼ 2–2.4 µg m⁻³ (daily maximum difference) for E2. Although there were no measurements to compare with modeled North Pole results, higher PM_2.5,sulf enhancements in North Pole are consistent with higher PM_2.5 concentrations in North Pole compared to Fairbanks (ADEC, 2017).

3.2 Modeled PM_2.5,sulf formation enhancements over the Northern Hemisphere in winter

To investigate model performance with the addition of heterogeneous sulfur chemistry for other locations and time periods that experience dark and cold conditions, Hemispheric CMAQ (HCMAQ) was run over the Northern Hemisphere from October 2015 to February 2016 (with the first 2 months as model spin-up) for the same base and sensitivity simulations as above. PM_25,sulf in China was of interest in this domain, given its haze events. High secondary and heterogeneously formed ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS have been documented to coincide with the high SO₂ emissions, high PM loadings, and high ALW during Chinese haze events (Cheng et al., 2016; Elser et al., 2016; Li et al., 2017; Ma et al., 2020; Peng et al., 2021; Wang et al., 2016). Winter episode average modeled SO₂ over northeastern China and other parts of China can be > 20 ppb, episode average modeled PM_2.5 can be > 80 µg m⁻³, and episode average ALW can be > 80 µg m⁻³ in this region (Fig. S8).

Maximum enhancements in PM_2.5,sulf for the Base_ Het simulation occurred largely in the North China Plain and northeastern China, with some notable enhancement over India as well (Fig. 6), and led to a maximum daily increase in PM_2.5,sulf by up to ∼ 54 µg S m⁻³ at a grid cell in southern China (27.1651° N, 107.5234° W; grid cell 149, 67) for this winter period (Fig. 6).

Figure 6 Episode average sulfate (a) and HMS (c) concentrations in the base simulation along with maximum daily differences in sulfate (b) and HMS (d) concentrations between the Base_Het and base CMAQ simulations over the Northern Hemisphere for a wintertime episode (from 1 December 2015 to 29 February 2016). HMS formation was not included in base CMAQ. Differences are cast in micrograms of sulfur per meter cubed to be consistent with measurement units from EMEP.

This maximum daily enhancement is due almost entirely to an increase in predicted ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (max daily concentration of 53 µg S m⁻³ at this same location). On average, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ can increase by up to ∼ 9 µg S m⁻³ over a grid cell in northeastern China (45.6698° N, 127.9877° W; grid cell 122, 64). HMS contributions to PM_2.5,sulf for the Base_Het run over this domain were less significant, with a maximum daily concentration of ∼ 2.6 µg S m⁻³ in a grid cell near Tehran, Iran (36.7976° N, 51.6208° W; grid cell 137, 119).

For the TMI_sens simulation, maximum daily PM_2.5,sulf concentrations increased from the base simulation by up to ∼ 33 µg S m⁻³ and on average up to ∼ 7 µg S m⁻³ at a grid cell in Hebei, China (39.8603° N, 119.2348° W; grid cell 131, 65) (a reduced episodic enhancement in comparison to the Base_Het simulation). This maximum enhancement is also almost entirely attributed to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ increases. However, HMS concentrations in the TMI_sens simulation contributed up to ∼ 2.8 µg S m⁻³ (in maximum daily concentrations) in a grid cell in northeastern China (45.6698° N, 127.9877° W; grid cell 122, 64) (Fig. S9). Predicted PM_2.5,sulf concentrations in the TMI_NO₂_sens and All_Ionic simulations were similar to the TMI_sens simulation for both predicted ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and predicted HMS concentrations, with the exception of not reproducing as high HMS concentrations in northeastern China (Fig. S9).

The addition of heterogenous sulfur chemistry decreased the model bias during an extreme haze event on the HCMAQ domain as well. Modeled PM_2.5,sulf concentrations from a grid cell over Beijing were compared to sulfate measurements at Tsinghua University in Beijing (Zhang et al., 2021a) (Fig. 7).

Figure 7 Model–measurement comparisons of PM_2.5,sulf across all model runs for HCMAQ at a grid cell over Tsinghua University in Beijing, China, from 5 to 22 December 2015. Observations are sourced from Zhang et al. (2021b).

The Base_Het and all additional sensitivity runs predicted higher PM_2.5,sulf at this grid cell than the base model simulation and reduced modeled mean bias by 2.9 µg S m⁻³ (model mean bias with base was −4.2 µg S m⁻³, and mean bias with Base_Het was −1.3 µg S m⁻³) (Fig. 7). Despite the overall improvement in model performance in the Base_Het simulation, a substantial gap in modeled and measured PM_2.5,sulf still existed on 22 December. Daily averaged modeled SO₂, NO₂, HCHO, and TMI concentrations (from base HCMAQ, representing a lower bound for SO₂ consumption) for this time period show that peak PM_2.5,sulf concentrations coincide with the co-occurrence of heightened SO₂+ TMI + NO₂ concentrations (Fig. S10). On 22 December, while SO₂ concentrations reached a daily average of ∼ 22 ppb, NO₂ and TMI concentrations were ∼ $\mathrm{1}/\mathrm{2}$ of the concentrations they were on 7–8 December.

The enhancement in PM_2.5,sulf is mainly attributed to SO${}_{\mathrm{4}}^{\mathrm{2}-}$ increases, as HMS maximum daily concentrations over this grid cell for this time period were ∼ 0.002 µg S m⁻³ for Base_Het and the three sensitivity simulations. Simulated HMS for Base_Het and the sensitivities for the winter 2016 case were also ∼ 0.002 µg S m⁻³ (∼ 0.006 µg m⁻³). While the HCMAQ estimates may be low, they are within the range of daily filter measurements of a polluted episode (from ∼ 0 to 1.2 µg m⁻³) reported in Ma et al. (2020).

Normalized mean biases (NMBs) for PM_2.5,sulf in the base simulations over the Northern Hemisphere ranged between −90 % and 800 %, and NMB for the Base_Het simulation ranged from −90 % to 980 % (Fig. S11). The largest positive NMB occurred at a site in the western US for both simulations and is due to very low modeled and observed concentrations. The mean biases at this location were only ∼ 0.1–0.2 µg S m⁻³ for both simulations, as this region was generally not affected by our updates (Fig. 6). Improvement of negative NMB occurred in the Base_Het run in some parts of the eastern US and Canada; however it caused and/or increased NMB in the positive direction in most of the eastern US and Canada (Fig. S11). The change in NMB in Europe was not significant, despite enhancements from HMS formation. Regionally aggregated model performance metrics can be found in Table S1.

3.3 CMAQ model performance changes over the Contiguous United States with the implementation of heterogeneous sulfur chemistry

CMAQ was run over the CONUS domain for both a winter- and a summertime episode to ensure that heterogeneous sulfur chemistry updates intended to be significant during dark and cold episodes had minimal impact on domains and episodes where this chemistry is less likely to dominate. For a January 2016 simulation over the CONUS domain, PM_2.5,sulf enhancements can mainly be seen over the eastern US and western Canada (Fig. 8) and are on average enhanced by up to ∼ 1.5 µg m⁻³. Average daily enhancements in Base_Het ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ for the entire episode were up to ∼ 1.5 µg m⁻³ as well, and daily average Base_Het HMS enhancements were up to ∼ 0.7 µg m⁻³.

Figure 8 Episode average sulfate (a) and HMS (c) concentrations in the base simulation along with maximum daily differences in sulfate (b) and HMS (d) concentrations between the Base_Het and base CMAQ simulations over the Contiguous United States for a wintertime episode (January 2016). HMS formation was not included in base CMAQ, and HMS mass concentrations are multiplied by the ratio of sulfate to HMS molecular mass.

However, PM_2.5,sulf maximum daily enhancements for this episode reached 28 µg m⁻³ at a grid cell in southwestern Kansas (37.3983° N, 101.9184° W; grid cell 121, 177). It should be noted that total modeled daily averaged PM_2.5 concentrations at this grid cell were 575 µg m⁻³, signifying a major PM event. PM_2.5,sulf enhancements from ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ at this grid cell and on this day contributed ∼ 75 %, with HMS contributing ∼ 25 %. However, the leading PM_2.5,sulf formation pathway at this grid cell in the Base_Het run was the TMI-catalyzed O₂ oxidation pathway in ALW (Fig. S12). This pathway dominates at a few other locations (Fig. S12); however, gas-phase oxidation of SO₂ by OH is the leading secondary PM_2.5,sulf formation pathway spatially, followed by cloud–aqueous oxidation by H₂O₂ and HNO₄.

The maximum daily average enhancement in HMS was ∼ 13 µg m⁻³ and occurred in the Ozarks in south-central Missouri (36.6721° N, 92.9428° W; grid cell 114, 243), also coinciding with a major PM event (model daily average PM_2.5 concentrations of 301 µg m⁻³) (Fig. 8).

PM_2.5,sulf concentration enhancements in the TMI_ sens simulation are similar to those in Base_Het, with a lower contribution from ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ production and a higher contribution from HMS production (Fig. S12a and b). The highest maximum daily enhancement for the TMI_sens simulation occurred in the same grid cell in southwestern Kansas discussed previously and was slightly lower than the Base_Het simulation's maximum daily average enhancement (27.7 µg m⁻³). However, the percentage of PM_2.5,sulf that was ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS was ∼ 10 % (3.2 µg m⁻³) and ∼ 90 % (28.3 µg m⁻³ or 24.5 µg m⁻³ SO₄ eq.). Other instances of HMS concentrations higher than 5 µg m⁻³ for this run were infrequent and coincided with major PM_2.5 events (PM_2.5 concentrations > 100 µg m⁻³). The highest maximum daily enhancements in PM_2.5,sulf for the TMI_NO₂_sens and All_Ionic runs were 22 and 18 µg m⁻³ respectively and occurred at the same grid cell in southwestern Kansas. Enhancements in PM_2.5,sulf in this grid cell were mainly attributed to HMS for both the TMI_NO₂_sens and the All_Ionic simulations and were ∼ 14–15 µg m⁻³, further demonstrating the importance of HMS to total PM_2.5,sulf when a temperature- and pH-dependent TMI-catalyzed O₂ k_chem is used during a cold major PM event.

Both the base and the Base_Het CONUS simulations overestimate measurements in the western US (Fig. 9), with the highest positive NMB being ∼ 450 % in Washington state with little difference in the NMB between the two runs – as this region was generally not affected by the implementation of heterogeneous sulfur chemistry (Fig. 8). It should also be noted that mean bias in the western US is fairly low (Fig. S13), and therefore the high NMB in this region is a result of overall low PM_2.5,sulf concentrations (both modeled and measured). Negative NMB of PM_2.5,sulf in the eastern US in the base run are ameliorated in the Base_Het runs, particularly in the Ohio River valley and southeastern US; however, areas in the east with good model PM_2.5,sulf performance in the base run were mostly overpredicted in the Base_Het run. Model performance differences across the additional sensitivity runs were minimal; however, all sensitivity runs had a lower positive NMB at the AQS sites around Chicago compared to Base_Het (Fig. S14), and the TMI_sens simulation had the lowest positive NMB and mean bias at the AQS site in south-central Missouri out of all the newly implemented model runs.

Figure 9 Normalized mean bias in PM_2.5,sulf concentrations by monitor over the CONUS domain for a wintertime episode (January 2016) for base CMAQ (a) and Base_Het CMAQ (b).

For a July 2016 simulated episode over the CONUS domain, overall daily enhancements in Base_Het PM_2.5,sulf were on average 0.05 µg m⁻³ across the entire domain. Daily ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ enhancements were more prevalent spatially over the eastern part of the US, particularly in the south and Ohio River valley (Fig. 10), and HMS formation was more prevalent in the western part of the domain.

Figure 10 Episode average sulfate (a) and HMS (c) concentrations in the base simulation along with maximum daily differences in sulfate (b) and HMS (d) concentrations between the Base_Het and base CMAQ simulations over the Contiguous United States for a summertime episode (July 2016). HMS formation was not included in base CMAQ, and HMS mass concentrations are multiplied by the ratio of sulfate to HMS molecular mass.

Maximum daily PM_2.5,sulf enhancements in the Base_Het model run reached up to 789 µg m⁻³, at a grid cell in Monterey County, California (36.2839° N, 121.9208° W; grid cell 135, 30), on 26 July. This was largely due to high daily averaged HMS concentrations (∼ 533 µg m⁻³) (Fig. 10) and coincided with the Soberanes fire at this location in late July (Queally, 2016). ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration enhancements were also significant at this grid cell (∼ 255 µg m⁻³). It should be noted that modeled daily averaged PM_2.5 and HCHO concentrations for this day and grid cell were also extremely high (9389 µg m⁻³ and 577 ppbv respectively).

In the same grid cell and at the same time as the occurrence of maximum daily enhancements in PM_2.5,sulf in the Base_Het run, TMI_sens-simulated PM_2.5,sulf increased by ∼ 791 µg m⁻³, with HMS contributing ∼ 732 µg m⁻³. Similar to results in other domains and episodes, in the TMI_sens simulation, HMS formation is higher, and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ enhancement decreases in comparison to Base_Het. Average daily ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ enhancements in TMI_sens reached up to 4 µg m⁻³ and in Base_Het reached up to 11 µg m⁻³, and both occurred in the same grid cell in Monterey, California, mentioned before. For the TMI_NO₂_sens simulation, maximum daily enhancements in PM_2.5,sulf concentrations are similar to TMI_sens (∼ 791 µg m⁻³), with a slightly lower contribution from HMS (∼ 720 µg m⁻³). Maximum daily PM_2.5,sulf enhancements were the lowest in the All_Ionic model run but were still a substantial enhancement (∼ 788 µg m⁻³; occurring in the same grid cell and at the same time as all other model run maximum daily enhancements), with HMS contributing ∼ 709 µg m⁻³ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ contributing ∼ 79 µg m⁻³ to this enhancement. The spatial distributions of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ enhancements for all sensitivity runs were similar to those of the Base_Het ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ enhancements (Fig. S15).

Both the base and the Base_het simulations mainly overestimate PM_2.5,sulf concentrations in the northwestern US and have similar performance (Figs. 11 and S16) due to this region generally not being affected by heterogeneous sulfur chemistry updates (with the exception of a few locations that are not in the same grid cells as an AQS monitor). Although there were significant enhancements in HMS and off the Californian coast, HMS concentrations greater than 5 µg m⁻³ were limited to a few grid cells on 26 July (when maximum concentrations of HMS were predicted) and did not reach the nearest AQS monitors. Higher daily HMS concentrations were predicted south of Monterey down the Californian coast the following 2 d (which did not have a corresponding AQS measurement). Increased ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation in the Base_Het simulation improved modeled underestimates in some parts of the eastern US (Figs. 11 and S16), particularly in parts of the southeast (Fig. S17).

Figure 11 Normalized mean bias in PM_2.5,sulf concentrations by monitor over the CONUS domain for a summertime episode (July 2016) for base CMAQ (a) and Base_Het CMAQ (b).

Although the implementation of heterogeneous sulfur chemistry generally led to minimal enhancements in both ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS over the entire CONUS domain for both episodes, it did slightly increase ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in the eastern US and substantially increased HMS concentrations during major PM events, including but not limited to wildfires. Enhancements in ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ could potentially impact the formation of secondary organic aerosol (SOA; Fan et al., 2022). ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is an important nucleophile in the formation of isoprene and epoxydiol (IEPOX)-derived organosulfates (Surratt et al., 2010), and enhancements in ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration may help reduce model underpredictions for IEPOX organosulfates (Budisulistiorini et al., 2017). Aqueous ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in the presence of Fe³⁺ and other isoprene-derived SOA precursors has been shown to enhance the production of C2–C4 organosulfur compounds as well (Huang et al., 2019). Modeled enhancements in both ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and HMS during wildfires are both particularly high due to extraordinarily high emissions during these events in general. TMI_sens predicted the highest HMS concentrations across all model simulations; however, max daily ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in TMI_NO₂_sens were higher in comparison to TMI_sens, indicating the importance of the heterogeneous NO₂ oxidation pathway during wildfire events. Wildfires release high concentrations of black carbon, which can inhibit modeled photolysis of NO_2,, and high in-plume oxidant concentrations can facilitate its regeneration (Buysse et al., 2019; Xing et al., 2017).

4 Discussion

4.1 Improved model performance with heterogeneous sulfur chemistry in dark, cold episodes

Traditional mechanisms of secondary sulfate formation have not been able to reproduce the high observed sulfate concentrations experienced during Fairbanks and North Pole, AK, winters. The added heterogeneous sulfur chemistry in ALW in this study enhanced modeled wintertime PM_2.5,sulf concentrations in this domain and in China, where high wintertime PM_2.5,sulf concentrations have also been observed (ADEC, 2017; Cheng et al., 2016; Ma et al., 2020; Moch et al., 2018; Wang et al., 2016; Zhang et al., 2021a). Inclusion of these reactions improved model performance in these regions (Figs. 5 and 7). PM_2.5,sulf enhancements between the base and Base_Het simulations during a PM episode (E1) in the Fairbanks domain showed that heterogeneous sulfur chemistry in ALW can increase daily PM_2.5,sulf concentrations by up to 11 µg m⁻³ across the entire domain for this dark and cold modeling episode (Fig. 2). At grid cells in downtown Fairbanks and North Pole, maximum daily enhancements in PM_2.5,sulf during E1 can range between ∼ 1–1.5 and ∼ 4–25 µg m⁻³ respectively (Figs. 2 and S3) across all of the heterogeneous sulfur chemistry configurations tested in this work. In this same domain during E2, modeled secondary PM_2.5,sulf formation takes place primarily via in-cloud formation pathways (Figs. 4 and S7), highlighting the enduring importance of cloud–aqueous chemistry in modeling PM_2.5,sulf formation when substantial liquid cloud water is present.

PM_2.5,sulf enhancements from heterogeneous sulfur chemistry for the wintertime HCMAQ simulations occurred mostly in China, South Asia, and Europe (Fig. 6). Maximum daily enhancements in PM_2.5,sulf ranged between 33–54 µg S m⁻³ across heterogeneous sulfur chemistry configurations tested in this work and were all located in China. Heterogeneous sulfur chemistry updates also improved negative model bias for PM_2.5,sulf in Canada (Table S1), where wintertime temperatures can be extremely cold and residential home heating along with industrial emissions may be high (Cho et al., 2009; Liggio et al., 2016; Stroud, 2009).

In general, there were minimal changes in predicted PM_2.5,sulf in regions and during episodes that do not experience extremely cold and dark conditions. Positive model bias in estimating PM_2.5,sulf in Europe for the HCMAQ simulation increases slightly in the Base_Het, TMI_sens, and TMI_NO₂_sens simulations but remains unchanged in the All_Ionic simulation. Similarly, model bias in estimating PM_2.5,sulf in the US for HCMAQ also increases slightly from the base in all the newly implemented model runs by 0.05 µg m⁻³ (Table S1).

The inclusion of heterogeneous sulfur chemistry had a smaller impact on modeled PM_2.5,sulf over the CONUS domain (in comparison to the Fairbanks and HCMAQ domains), with the exception of high-PM events. Interestingly, the inclusion of these missing pathways reduced some of the negative bias in PM_2.5,sulf concentrations in the southeast during warm and light summertime conditions (Figs. 11, S16, and S17), and this bias reduction is mostly due to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ enhancements (Figs. 10 and S15). ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ can participate as a reactant and a modulator of pH in the heterogeneous formation of IEPOX-derived organosulfates (Marais et al., 2016; Pye et al., 2017; Pye et al., 2020), and improved ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ negative bias may also improve model performance of IEPOX SOAs.

4.2 Contribution of HMS to PM_2.5,sulf

In this study we found that HMS, a previously untracked aerosol species in CMAQ, can contribute substantially to total PM_2.5,sulf, depending on the HCHO and SO₂ concentrations, temperature, and heterogeneous sulfur chemical kinetics chosen. During E1, HMS concentrations in North Pole are higher in comparison to Fairbanks. Given the stagnant conditions for this domain and episode, emissions tend to stay local (ADEC, 2017; Gaudet et al., 2012; Tran and Molders, 2011). Thus, HCHO from residential wood combustion and SO₂ from home-heating oil in North Pole, a largely residential area, compared to downtown Fairbanks are the likely dominating reasons for higher HMS concentrations when comparing the two areas.

Modeled HMS generally had a smaller impact on both the northern hemispheric and the CONUS domains (generally contributing < 1 µg m⁻³), with a few anomalous exceptions. Modeled HMS in HCMAQ runs primarily appeared in Europe and China, where wintertime HMS has previously been predicted (Moch et al., 2020). In the CONUS domain, HMS concentrations were predicted to be high during major PM events (including but not limited to wildfires), which illuminates the importance of HMS during high atmospheric loadings of both SO₂ and HCHO.

Across all modeled domains and episodes, the TMI_sens model run predicted the highest HMS concentrations and may represent an upper bound on modeled HMS concentrations. Ultimately, more resolved measurements of speciated PM_2.5 that can separate HMS and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (Campbell et al., 2022) can help discern their relative contributions to PM_2.5,sulf mass and help constrain future modeled heterogeneous sulfur chemical kinetics.

4.3 PM_2.5,sulf formation pathways of interest during cold and dark episodes

In addition to the inclusion of both heterogeneous SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and both in-cloud and heterogeneous HMS formation in CMAQ, we determined which PM_2.5,sulf formation pathways are the most important given ionic strength, pH, and temperature regimes characteristic of dark and cold conditions. Across both the Fairbanks and the CONUS domains in the Base_Het simulation during wintertime, the most prevailing heterogeneous SO${}_{\mathrm{4}}^{\mathrm{2}-}$ formation pathway was the TMI-catalyzed O₂ pathway (Figs. 2, 4, S12). In the TMI_sens E1 in Fairbanks, however, this formation pathway was the third-most important, behind HMS formation and the NO₂ pathway (Fig. S3). Although the modeled pH for TMI_sens ranged between 3–6 for Fairbanks and North Pole and for both episodes (Fig. S4), which included the optimal pH for this pathway (pH = 4.2; Ibusuki and Takeuchi, 1987), the dampening of this pathway can also be attributed to the extremely cold temperatures (modeled average −30 °C or 243 °K).

TMI_sens-modeled aerosol pH was seen to be the least acidic in comparison to all of the other model simulations, especially in North Pole (Fig. S4). As noted before, HMS was the largest contributor to secondary PM_2.5,sulf formation in North Pole, the formation (and loss) rates of which increased with increasing pH (Ervens et al., 2003; Kok et al., 1986) (Fig. 2). Aerosol pH and ALW calculations in ISORROPIA II only consider inorganic species. Organic species (e.g., organic acids) may also increase aerosol acidity (Zuend and Seinfeld, 2012; Zuend et al., 2011), and therefore the predicted aerosol pH in TMI_sens might represent an overprediction. Aerosol pH for the Base_Het, TMI_NO₂_sens, and All_Ionic model simulations was similar at both North Pole and Fairbanks, with both sensitivity simulations predicting slightly higher pH than the Base_Het simulation during E1 and slightly lower pH during E2.

The impacts of increased ionic strength were explored with respect to the NO₂, O₃, and H₂O₂ oxidation pathways (note that ionic strength inhibition of the TMI-O₂ pathway is included in the Base_Het simulations as well). Ionic strength impacts added to the NO₂ pathway (in TMI_NO₂_sens) had the highest impact on the formation of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ via this pathway (Figs. S3 and S7c and d) in the Fairbanks domain for both episodes. Although the ionic strength for this pathway was bounded at 1.14 M, an increase in aerosol ionic strength from 0.1 M to the upper bound of 1.14 M increased k_chem for this pathway by ∼ 2 orders of magnitude (Chen et al., 2019). Ionic strength impacts on the H₂O₂ and O₃ heterogeneous sulfur oxidation pathways had minimal impact during the wintertime PM episodes in Fairbanks and North Pole, only accounting for ≤ ∼ 0.2 µg m⁻³ in the All_Ionic model simulation (Figs. S3 and S7). These particular pathways were not assumed to be prolific given the lack of photochemistry during this episode; however with an ionic strength change from 0 to 5 M, the third-order aqueous-phase rate coefficient for the H₂O₂ heterogeneous sulfur pathway can increase by more than 40 % regardless of pH or temperature (Maaß et al., 1999; Millero et al., 1989). The ionic strength used to calculate the ionic strength effect factor for this k_chem was limited to a maximum of 5 M; however recent studies have observed significant ionic strength enhancement up to 14 M (Liu et al., 2020). For ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ formation via heterogeneous oxidation by O₃, k_chem for this pathway can increase by ∼ 80 %, with an ionic strength increase from 0.1 to 0.8 M (Lagrange et al., 1994; Song et al., 2021a) at temperatures characteristic of Fairbanks winters; however these effects were not seen due to the lack of ozone modeled given the dark and cold conditions.

5 Conclusion

Air quality modeling of secondary sulfate has traditionally only included in-cloud aqueous- and gas-phase SO₂ oxidation pathways, often resulting in underpredictions of observed PM_2.5,sulf, especially during the cold and polluted conditions characteristic of wintertime PM and haze events (ADEC, 2017; Gao et al., 2016). In this study, we implemented heterogeneous sulfur chemistry in aerosol liquid water in CMAQ to resolve model–measurement gaps in PM_2.5,sulf concentrations during extreme wintertime PM episodes in and around Fairbanks, Alaska. We compared modeled PM_2.5,sulf (SO${}_{\mathrm{4}}^{\mathrm{2}-}+$ HMS) concentrations to sulfate measurements at several measurement sites (under the assumption that HMS may be included in the sulfate observations, Dovrou et al., 2019). Negative model bias improved in Fairbanks during winter and decreased with these updates and also improved in Beijing, another location known to experience wintertime haze events. When applied more broadly to larger domains and other seasons, the update also resolved underestimations of PM_2.5,sulf concentrations in both the United States and globally; however, it did not have a huge impact when applied over domains that were not as dark and cold. HMS was found to be an important contributor to PM_2.5,sulf mass during dark and cold episodes; however, to better understand the ratio of sulfate to HMS, more observations of HMS are necessary. Recently, the Alaskan Layered Pollution and Chemical Analysis (ALPACA) field campaign (Simpson et al., 2024) was conducted in and around Fairbanks during January–February 2022 and offers observations to elucidate important sulfate and HMS formation pathways in the area and better characterize source apportionment of PM_2.5,sulf. Ultimately, HMS and sulfate formed via the TMI-catalyzed O₂ and NO₂ pathways proved to be the most important to PM_2.5,sulf formation pathways in Fairbanks and North Pole and require further investigation in the context of PM_2.5,sulf control strategies. Finally, while this study aims to include ionic strength, pH, and temperature impacts on PM_2.5,sulf formation in CMAQ, the ionic strength, pH, and temperature ranges under which Henry's law, reaction rates, and other coefficients and parameters were derived experimentally may not be representative of concentrated aerosol water or that of extremely cold and dark wintertime conditions. Laboratory studies are needed to extend the bounds of these parameters and determine rate expressions appropriate for the concentrated conditions characteristic of aerosol water.