Interface-dominated hydroxymethanesulfonate and its isomer formation provides key mechanisms for reconciling the atmospheric sulfur budget gap in polluted and cold environments

Hydroxymethanesulfonate (HMS) is a critical source of particulate sulfur, formed by formaldehyde (HCHO) and sulfur dioxide (SO₂) in droplets. Current models relying on bulk aqueous-phase HMS formation only explain ∼ one-third of unexplained sulfate concentrations, leaving gaps in atmospheric sulfur budget, especially in polluted and cold environments. Using Born–Oppenheimer molecular dynamics simulations, we explored HMS and its isomer hydroxymethyl sulfite (HMSi) formation mechanisms across aqueous phase and air–water/ice interfaces. Air–water interfaces enable nearly barrierless HMS formation (0.6 kcal mol⁻¹) via unique stepwise water-mediated proton transfer, preferring HMS over HMSi (0.6 vs. 6.1 kcal mol⁻¹), which contrasts sharply with the competitive pathways observed in the bulk aqueous phase (7.7 vs. 7.6 kcal mol⁻¹). In contrast, protonation of formaldehyde under strongly acidic conditions reverses reaction selectivity, favoring HMSi formation over HMS. Importantly, these reaction mechanisms remain viable at air–ice interfaces in cold environments including polar areas and the upper troposphere, revealing ice surfaces as previously overlooked yet significant sites for atmospheric organosulfate formation. Our findings suggest that interfacial mechanisms may provide efficient pathways for HMS and HMSi formation in both polluted and cold environments, helping to reconcile model-observation discrepancies in the atmospheric sulfur budget.

1 Introduction

Atmospheric fine particulate matter (PM) impacts human health, visibility, ecosystems, and global climate (Zhang et al., 2015). These fine particles can exist not only as liquid aerosol droplets but also transform into ice crystals under cold atmospheric conditions (Knopf et al., 2018). Among these constituents, sulfates represent critical species driving atmospheric particle formation (Sipila et al., 2010; Zhang, 2010; Zhong et al., 2019). Despite decades of research, a persistent discrepancy between modeled and observed sulfate concentrations remains, particularly in polluted, strongly acidic aerosol environments such as winter haze episodes in northern China and remote cold regions (Wang et al., 2014, 2016, 2025; Campbell et al., 2024). This gap reveals critical deficiencies in our understanding of atmospheric sulfur cycling, particularly regarding organosulfur intermediates that may contribute significantly to sulfate formation.

Hydroxymethane sulfonate (HMS), formed from the reaction of formaldehyde (HCHO) and sulfur dioxide (SO₂) in cloud droplets (Reactions R1–R3), serves as an important tracer for cloud and fog processing (Wei et al., 2020). However, HMS has often been overlooked as a potential contributor to atmospheric sulfur budgets, primarily because common aerosol measurement techniques (e.g., commercial aerosol mass spectrometry) frequently misidentify HMS as sulfate, thereby compromising the accuracy of sulfate simulation models (Moch et al., 2018). With advances in analytical techniques, recent field studies have detected HMS in diverse environments worldwide, including urban areas such as Beijing, tropical regions such as Singapore, remote cold regions such as Alaska, and marine settings (Campbell et al., 2022; Moch et al., 2020; Song et al., 2019; Zhao et al., 2024). HMS impacts atmospheric chemistry not only through its widespread presence but also via its subsequent chemical transformations to sulfate (Lai et al., 2023, 2024). Multiple studies have demonstrated that HMS formation can substantially contribute to ambient sulfate concentrations (Huang et al., 2023; Zhang et al., 2024b). These findings highlight the potentially significant contribution of HMS to atmospheric sulfur chemistry. Song et al. (2019) estimated that HMS could account for approximately one-third of the discrepancy between model-simulated and observed sulfate in haze events.

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\mathrm{SO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons {\mathrm{SO}}_{\mathrm{2}}\cdot {\mathrm{H}}_{\mathrm{2}}\mathrm{O}\text{(R2)}& {\displaystyle }{\mathrm{SO}}_{\mathrm{2}}\cdot {\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons {\mathrm{HOSO}}_{\mathrm{2}}^{-}+{\mathrm{H}}^{+}\text{(R3)}& {\displaystyle }\mathrm{HCHO}+{\mathrm{HOSO}}_{\mathrm{2}}^{-}\to {\mathrm{CH}}_{\mathrm{2}}\left(\mathrm{OH}\right){\mathrm{SO}}_{\mathrm{3}}^{-}\left(\mathrm{HMS}\right)\text{(R4)}& {\displaystyle }\mathrm{HCHO}+{\mathrm{HOSO}}_{\mathrm{2}}^{-}\to {\mathrm{CH}}_{\mathrm{2}}\left(\mathrm{OH}\right){\mathrm{OSO}}_{\mathrm{2}}^{-}\left(\mathrm{HMSi}\right)\end{array}$$

Nevertheless, current studies on the HCHO and SO₂ reaction mechanisms suffer from critical limitations that hinder the closure of the sulfur budget gap, especially in polluted and cold atmospheric environments. Chen and Zhao (2020) employed density functional theory calculations with implicit solvent models to study aqueous-phase reactions between HCHO and bisulfite (HOSO${}_{\mathrm{2}}^{-}$). Their study revealed that this reaction can produce hydroxymethyl sulfite (HMSi, Reaction R4), a structural isomer of HMS with distinct sulfur-oxygen binding, as a competitive product in aqueous-phase (Chen and Zhao, 2020). However, implicit solvent models, which treat the solvent as a structureless continuum, inherently lack the molecular-level detail necessary to capture subtle mechanistic variations and dynamic solvent-solute interactions (Zhang et al., 2017). In contrast, explicit solvent models allow reactive intermediates and transition states to interact directly with individual solvent molecules, fully capturing the dynamic influence of solvation on reaction pathways (Norjmaa et al., 2021). This gap in molecular-level understanding of HMS/HMSi aqueous formation mechanisms remains a nonnegligible barrier to accurately modeling atmospheric sulfur cycling.

In addition to bulk aqueous phase, air-water interfaces represent ubiquitous yet fundamentally distinct reactive environments in atmospheric systems, including cloud droplets, aerosols, and sea-salt particles. The molecular environment at these interfaces differs significantly from bulk aqueous phases, characterized by unique hydrogen bonding networks that can enhance the adsorption and condensation of atmospheric species onto aerosols (Li et al., 2021; Ruiz-Lopez et al., 2020). Notably, both HCHO and bisulfite (HOSO${}_{\mathrm{2}}^{-}$) demonstrate preferential accumulation at air-water interfaces (Martins-Costa et al., 2012; Yang et al., 2019), creating localized concentration enhancements and reactive microenvironments. These distinct interfacial properties may substantially enhance reaction kinetics or facilitate entirely novel reaction mechanisms compared to bulk solution (Ning et al., 2023, 2024). Despite their potential significance, a systematic understanding of how HMS and HMSi formation mechanisms differ between bulk and interfacial environments and whether reaction selectivity toward HMS versus HMSi varies at air-water interfaces, remains absent. This knowledge gap directly limits the accuracy of organosulfur formation representations in atmospheric models.

Additionally, the chemical complexity introduced by extreme aerosol acidity (pH) presents another critical knowledge gap. Based on global surface-layer aerosol pH estimated by Li et al. (2022) using annual-mean (2016) GEOS-Chem simulations coupled with E-AIM thermodynamic calculations, aerosol pH in the range of −1 to 1 occurs over at least ∼ 20 % of the global surface area (Li et al., 2022) where the model successfully converges (RH ≥ 60 % and T ≥ 263.15 K). These highly acidic regions span diverse environments including marine and continental systems (Angle et al., 2021; Ding et al., 2019; Jia et al., 2020; Kakavas et al., 2021; Pye et al., 2018; Zheng et al., 2020; Zhou et al., 2022). These highly acidic environments trigger fundamental changes in formaldehyde speciation, leading to the formation of protonated formaldehyde (HCHOH⁺). Experimental evidence (Jayne et al., 1996) has demonstrated that HCHOH⁺ can exist and accumulate at droplet surfaces under such acidic conditions, establishing its potential significance in atmospheric reaction. Despite this potential significance, the interfacial chemistry of protonated formaldehyde with sulfur species remains uncharacterized. Hence, the contributions of HMS and HMSi to atmospheric sulfur budgets in these acidic environments also remain poorly understood.

In addition to air-water interfaces, ice crystals exhibit high surface-to-volume ratios, providing extensive reactive interfaces that play crucial roles in atmospheric chemistry (George et al., 2015; Kerbrat et al., 2008; Zhong et al., 2020). Notably, HCHO and SO₂, the precursors to HMS and HMSi, are ubiquitous in cold regions (Sumner and Shepson, 1999; Anderson et al., 2017; Höpfner et al., 2015; Joppe et al., 2024). However, it remains unclear whether air–ice interfaces favor HMS or HMSi formation. If ice surfaces do facilitate HMS or HMSi formation, current models, which overlook the organosulfur chemistry on ice surface, are likely to substantially underestimate their contributions to atmospheric sulfur budgets.

To address these critical knowledge gaps, we employed Born–Oppenheimer molecular dynamics (BOMD) simulations with metadynamics enhancement and high-level quantum chemical calculations to elucidate HMS and HMSi formation mechanisms in polluted and cold environments. Our specific objectives are: (i) to reveal aqueous-phase reaction mechanisms with explicit solvation, focusing on the competitive HMS/HMSi formation pathways; (ii) to clarify the formation mechanisms and selectivity rules of HMS/HMSi at the air-water interface; (iii) to determine how protonated formaldehyde alters reaction pathways and selectivity in strongly acidic environments; (iv) to assess the feasibility of HMS and HMSi formation at the air–ice interface at temperatures characteristic of polar regions and the upper troposphere; (v) to systematically characterize the interfacial dynamics of HMS and HMSi products, including their hydrogen bonding networks, orientational preferences, and implications for atmospheric oxidation. By integrating these investigations, we aim to reveal interfacial sulfur cycling mechanisms, reconcile model-observation discrepancies, and improve organosulfur chemistry representations in global atmospheric models.

2 Methods

2.1 Quantum Chemistry Calculations

Geometry optimizations of stationary points (reactants, products, and transition states) were performed at the B3LYP/6-311+G(d,p) level (Krishnan et al., 1980; McLean and Chandler, 1980; Stephens et al., 1994) with the DFT-D3 dispersion correction of Grimme et al. (2010) using Gaussian 16 package (Frisch et al., 2016). Harmonic frequency analyses were conducted at the same level. The transition states are verified to be connected with the corresponding reactants or products by means of intrinsic reaction coordinate (IRC) calculations. The solvation model based on density implicit solvent model (SMD) was applied to represent the aqueous phase. To obtain more accurate electronic energies, single-point energy calculations were subsequently performed at the DLPNO-CCSD(T) level of theory employing the aug-cc-pVTZ basis set (Guo et al., 2018; Riplinger and Neese, 2013; Riplinger et al., 2013). These calculations were carried out using the ORCA 5.0.4 quantum chemistry program package (Neese, 2011). Fukui function analysis evaluated the nucleophilic reactivity of atoms within the HOSO${}_{\mathrm{2}}^{-}$ anion. Electron density difference maps were analyzed to reveal charge transfer between reactant fragments. All wave function analyses were performed by Multiwfn 3.8 (dev) code (Lu and Chen, 2012, 2023; Zhang and Lu, 2021).

2.2 Born–Oppenheimer Molecular Dynamics Simulations

Born–Oppenheimer molecular dynamics (BOMD) simulations were performed using the CP2K package (Kuhne et al., 2020). All BOMD simulations employed the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional (Perdew et al., 1996) with the DFT-D3 dispersion corrections of Grimme et al. (2010). This D3 correction has been shown to improve the description of water at the water–air interface at the GGA level (Dodia et al., 2019). A double-ζ Gaussian basis set (DZVP-MOLOPT-SR-GTH) was used for valence electrons in conjunction with Goedecker–Teter–Hutter (GTH) pseudopotentials for core electrons (Goedecker et al., 1996; VandeVondele and Hutter, 2007). The plane wave and Gaussian basis set cutoff energies were set to 300 and 40 Ry, respectively. The simulations were conducted in the canonical (NVT) ensemble using the canonical sampling through velocity rescaling (CSVR) thermostat (Bussi et al., 2007) with a time step of 1 fs. For aqueous-phase simulations, a cubic box of 1.7 nm × 1.7 nm × 1.7 nm containing 128 water molecules was employed. The air–water interface model was constructed by expanding the z-axis of this cubic box to 4.7 nm while keeping the xy-plane dimensions (1.7 nm × 1.7 nm) and water molecules (128) unchanged, which generated vacuum regions above and below the aqueous slab to establish two air–water interfaces. Free energy profiles were computed via well-tempered metadynamics (see Sect. S1 in the Supplement for details) as implemented in PLUMED (Bonomi et al., 2009). Collective variables for HMS and HMSi formation are defined in Fig. S1 in the Supplement. Convergence was assessed by monitoring Gaussian hill height evolution (Fig. S2), collective variable sampling behavior (Fig. S3), and the time evolution of free energy profiles (Fig. S4). Detailed configurations for air–ice interface simulations are provided in the Supplement (Fig. S5).

3 Results and Discussion

3.1 Bulk-Phase Mechanisms with Explicit Solvation

Bulk aqueous reactions represent an important pathway for HCHO–HOSO${}_{\mathrm{2}}^{-}$ interactions in atmospheric chemistry. Although previous DFT studies employing implicit solvation models have provided valuable understanding of this reaction (Chen and Zhao, 2020), explicit solvation approaches offer complementary insights by capturing discrete water-reactant interactions and dynamic hydrogen bonding that may play important roles in determining reaction mechanisms and kinetics. To provide a more detailed molecular-level picture of bulk aqueous HMS formation, we employed BOMD simulations with metadynamics to investigate the HCHO–HOSO${}_{\mathrm{2}}^{-}$ reaction with explicit water molecules.

Figure 1 Snapshot structures captured from the metadynamics-biased BOMD simulations of HCHO and HOSO${}_{\mathrm{2}}^{-}$ reaction in the bulk phase. (a) Free energy profile for Path I. (b) Snapshot structures along Path I. (c) Free energy profile for Path II. (d) Snapshot structures along Path II. RC, TS, and P denote the reactant complexes, transition states, and products, respectively.

Given that HOSO${}_{\mathrm{2}}^{-}$ contains multiple potential nucleophilic sites (sulfur, non-hydroxyl oxygens, and hydroxyl oxygen), we first assessed their relative reactivities through preliminary quantum chemical calculations (Fig. S6). These analyses confirmed that the hydroxyl oxygen O(H) in HOSO${}_{\mathrm{2}}^{-}$ exhibits significantly lower nucleophilic reactivity compared to the sulfur atom and non-hydroxyl oxygen atoms. Furthermore, reactions involving the O(H) did not yield viable products. Based on these findings, reactions involving the O(H) were excluded from subsequent metadynamics simulations. As illustrated in Fig. 1, two distinct reaction pathways dominate the bulk-phase HCHO + HOSO${}_{\mathrm{2}}^{-}$ reaction. In Path I (Fig. 1b), the sulfur atom of the HOSO${}_{\mathrm{2}}^{-}$ attacks the carbon atom of formaldehyde, forming a seven-membered ring transition state with the participation of a water molecule. The calculated free energy barrier for the formation of HMS is 7.7 kcal mol⁻¹. Path II (Fig. 1d) proceeds via nucleophilic attack of a non-hydroxyl oxygen atom of HOSO${}_{\mathrm{2}}^{-}$ on the carbonyl carbon atom of HCHO. This pathway proceeds through a ring transition state mediated by water molecules (Fig. 1d). The free energy barrier for HMSi formation via this route is 7.6 kcal mol⁻¹. These values are in close agreement with barriers derived from high-level quantum chemical calculations (Fig. S7), validating the reliability of the simulation results. The comparable barriers for HMS and HMSi formation in bulk water (7.7 vs. 7.6 kcal mol⁻¹) indicate that both isomers form competitively in bulk water. Despite these moderate free energy barriers, bulk phase HMS formation alone is insufficient to reconcile the atmospheric sulfur budget gap. Atmospheric modeling studies incorporating bulk-phase HMS chemistry can account for only approximately one-third of the discrepancy between modeled and observed sulfate concentrations during severe haze episodes (Song et al., 2019). This shortfall may be partly attributed to the high surface-to-volume ratios of atmospheric aerosols and cloud droplets, where interfacial processes may play a more significant role than bulk-phase reactions.

3.2 Rapid and Selective HMS Formation at the Air-Water Interface

Given the limitations of bulk-phase chemistry identified above, understanding HCHO-HOSO${}_{\mathrm{2}}^{-}$ reactivity at air–water interfaces becomes critical. Recent investigations reveal significant accumulation of both HCHO and HOSO${}_{\mathrm{2}}^{-}$ at the air-water interface (Martins-Costa et al., 2012; Yang et al., 2019). This interfacial enrichment suggests potentially significant atmospheric reaction mechanisms previously overlooked. Herein, we employed BOMD simulations to investigate these interfacial heterogeneous reactions at the air-water interface.

Similar to bulk-phase reactions, interfacial processes proceed via two distinct pathways. Remarkably, HMS formation via Path I proceeds rapidly at the air–water interface. Starting from reactant complex configurations optimized at the B3LYP-D3/6-311+G(d,p) level, unbiased BOMD trajectories reveal that the reaction completes within ∼ 3 ps (Figs. 2a, S8), in contrast to the bulk-phase reaction. During the reaction, the distance between the C atom of HCHO and the S atom of HOSO${}_{\mathrm{2}}^{-}$ gradually decreases to approximately 1.95 Å (Fig. 2d, e), with a corresponding Mayer bond order of 0.57, indicating that an S–C bonding interaction has been established at this intermediate stage. The resulting anionic intermediate structurally resembles HMS, differing only in proton location, which remains on the O atom of the HOSO${}_{\mathrm{2}}^{-}$. Electron density difference analysis demonstrates decreased electron density surrounding the S atom with concurrent increase at the formaldehyde O atom during intermediate formation (Fig. 2f). Mulliken population analysis quantifies this charge transfer: the charge on S atom in HOSO${}_{\mathrm{2}}^{-}$ increases from approximately 0.6 to 0.81 e, while the charge on O atom in HCHO decreases from approximately −0.3 to −0.57 e (Fig. S9). These analyses can provide clear evidence for electron transfer from HOSO${}_{\mathrm{2}}^{-}$ to HCHO in the initial nucleophilic addition step. Subsequently, at approximately 2.42 ps, a transition state structure forms wherein the O2 atom of proximal H₂O approaches the H2 atom of the reaction intermediate, while the H1–O2 bond in H₂O lengthens. Interfacial water molecules facilitate the final step of the reaction, completing HMS formation while shortening the S–C bond from approximately 1.95 to 1.81 Å (Fig. 2g). Unlike the bulk-phase reaction, where S–C bond formation and proton transfer occur synchronously (Fig. 1b), the interfacial HMS formation proceeds via a stepwise mechanism: the S–C bond forms initially, followed by water molecule-mediated proton transfer from HOSO${}_{\mathrm{2}}^{-}$.

Figure 2 Mechanistic investigation of HCHO and HOSO${}_{\mathrm{2}}^{-}$ at the air-water interface. (a) Time evolution of key bond distances during HMS formation, with corresponding structural snapshots shown in (d). (b) Free energy profile for the HMS formation via Path I. (c) Free energy profile for the HMSi formation via Path II, with representative structural snapshots presented in (h). (e) Color-mapped electron localization function analysis illustrating bond reorganization in key intermediates during HMS formation. (f) Electron density difference (Δρ) of the reaction intermediates (gray region: Δρ>0; cyan region: Δρ<0). (g) Combined distribution functions involving the bond distance (r, Å) and angular distribution function (α, degree) for HCHO and HOSO${}_{\mathrm{2}}^{-}$.

To obtain quantitative free energy profiles and systematically explore the interfacial reaction pathways, we performed computationally expensive metadynamics-biased BOMD simulations. Metadynamics simulations reveal an exceptionally low free energy barrier of 0.6 kcal mol⁻¹ for interfacial HMS formation (Fig. 2b), corroborating the fast reaction observed in unbiased trajectories. We validated our BOMD findings using high-level quantum chemical calculations using a ten-water-molecule cluster model (W₁₀) at the DLPNO-CCSD(T)/aug-cc-pVTZ//B3LYP-D3/6-311+G(d,p) level. These high-level calculations further validated the near-barrierless HMS formation, providing complementary evidence to our metadynamics simulations (Fig. S10).

Our results differ from Li et al. (2025)'s recent report of a ∼ 7.6 kcal mol⁻¹ free energy barrier for this process. This discrepancy stems from their use of a single distance collective variable d(C⋯S) in metadynamics simulations, which overlooked the critical water-mediated cyclic proton transfer mechanism. In contrast to the rapid formation of HMS, no fast formation of HMSi was observed across multiple simulation trajectories. Metadynamics simulations revealed a higher free energy barrier of 6.1 kcal mol⁻¹ for HMSi formation at the interface (Fig. 2c), with the reaction mechanism paralleling that observed in the bulk phase. Although this barrier represents a 1.5 kcal mol⁻¹ reduction from the bulk-phase value (7.6 kcal mol⁻¹), it remains substantially higher than the negligible barrier for HMS formation.

To elucidate the potential contributor to the enhanced interfacial reactivity, we analyzed the solvation environments of HCHO and HOSO${}_{\mathrm{2}}^{-}$ at the interface versus bulk water (Fig. S11). The partial solvation characteristic of the air-water interface results in reduced water coordination around both reactants relative to their bulk-phase counterparts, potentially facilitating reactant approach. Furthermore, the HOMO–LUMO gap of the reactive complex is smaller at the interface than in bulk water, suggesting that the partially solvated environment may contribute to the reaction acceleration (Fig. S12) (Ishiyama et al., 2022; Kusaka et al., 2021). It is important to note that the interfacial acceleration mechanism is multifaceted. Beyond partial solvation effects, other factors including the interfacial electric field and synergistic enthalpy-entropy effects play crucial roles in promoting HMS formation at the air-water interface (Li et al., 2025). These computational findings demonstrate a clear mechanistic preference for HMS formation over HMSi at the air-water interface, with a markedly lower free energy barrier (ΔG=0.6 kcal mol⁻¹ for HMS versus 6.1 kcal mol⁻¹ for HMSi). The exceptionally low barrier indicates that HMS formation proceeds very rapidly, especially in the polluted regions with abundant HCHO and SO₂ in the atmosphere, which can provide an important theoretical foundation for more accurate simulation of atmospheric sulfur budgets.

3.3 Product Selectivity Reversal in Strongly Acidic Environments

The strongly acidic nature of many atmospheric environments, as discussed earlier, leads to protonation of HCHO, forming HCHOH⁺. Experimental studies (Jayne et al., 1996) have demonstrated that when pH drops below 2, formaldehyde undergoes efficient uptake by acidic droplets, leading to its protonation and formation of HCHOH⁺. However, the detailed mechanistic consequences of this protonation, especially as they manifest in heterogeneous environments, remain poorly understood. To capture the complex reactivity landscape of HCHOH⁺ with HOSO${}_{\mathrm{2}}^{-}$ under strongly acidic conditions, we conducted 25 independent BOMD simulations (detailed configurations in Fig. S13). Figure 3a displays the temporal evolution of product formation across all simulations, with each horizontal line representing an individual BOMD trajectory. The colored markers indicate specific chemical transformations. Hydroxymethanesulfonic acid (HMSA) and hydroxymethyl hydrogen sulfite (HMHSi) formed within 2 ps in most trajectories. The subsequent deprotonation processes led to the formation of HMS and HMSi within approximately 10 ps. The mechanistic diversity enabled by formaldehyde protonation under strongly acidic conditions represents a fundamental departure from reactions involving neutral formaldehyde molecules. Figure 3b illustrates the three distinct mechanistic pathways identified through our simulations. Path A proceeds via nucleophilic attack by sulfur atom of HOSO${}_{\mathrm{2}}^{-}$ at the carbonyl carbon of HCHOH⁺, leading to HMSA formation before transforming to HMS. Path B involves nucleophilic attack by either O1 or O2 of HOSO${}_{\mathrm{2}}^{-}$ on the carbonyl carbon, forming HMSi.

Figure 3 Mechanistic investigation of the reaction between HCHOH⁺ and HOSO${}_{\mathrm{2}}^{-}$ in acidic environment. (a) Time evolution of product formation across 25 independent BOMD simulations. Each horizontal line represents a single trajectory, with the starting point indicating the formation time of intermediate species (HMSA or HMHSi) and the end point indicating the formation time of final products (HMS or HMSi) after deprotonation. Colored markers represent different reaction products: HMSA (red), HMS (purple), HMHSi (yellow), and HMSi (green). Simulation time is shown in picoseconds (ps). (b) The three distinct mechanistic pathways (Paths A, B, and C).

Our trajectory analysis shows no significant preference between O1 and O2 attack, suggesting that both oxygen atoms experience similar activation in the interfacial acidic environment. The subsequent deprotonation of HMHSi to HMSi proceeds through water-mediated proton transfer, with interfacial water molecules playing a crucial catalytic role. Perhaps most remarkably, Path C demonstrates that the typically unreactive hydroxyl oxygen O(H) of HOSO${}_{\mathrm{2}}^{-}$ becomes a viable nucleophile under acidic conditions.

Although the O(H) site of HOSO${}_{\mathrm{2}}^{-}$ is generally weakly nucleophilic (Fig. S6), protonation fundamentally alters the electronic structure of HCHO by introducing a new δ-hole on the oxygen atom while substantially intensifying the δ-hole at the carbonyl carbon (Fig. S14). These multiple electron-deficient sites work synergistically to attract even weak nucleophiles such as O(H), thereby enabling direct HMSi formation. This reactivity highlights the unique ability of highly acidic interfaces to modulate electronic structure, thereby activating HCHOH⁺ toward unconventional nucleophilic pathways. Statistical analysis of the product distribution reveals a notable preference for HMSi formation under acidic conditions: among the 25 trajectories, 6 proceeded through Path A to form HMS, while 19 yielded HMSi via either Path B or Path C, establishing an HMS : HMSi ratio of approximately 1 : 3. This product distribution contrasts sharply with the preferential HMS formation observed at neutral pH interfaces, highlighting the profound influence of acidity on reaction selectivity. The pH-dependent reaction selectivity becomes even more pronounced under extreme acidic conditions. Under extremely acidic conditions (pH ≈ 0), the speciation of sulfurous acids at interfaces undergoes significant shifts. As observed by Buttersack et al. (2024) HSO${}_{\mathrm{3}}^{-}$ becomes the predominant species at the air-water interface under extremely acidic conditions. With the sulfur atom protonated and therefore unavailable for nucleophilic attack, S–C bond formation leading to HMS becomes inhibited, yielding only HMSi (Fig. S15). These results establish HMSi as a potentially significant contributor to atmospheric sulfur budgets in strongly acidic environments. From an analytical perspective, HMS and HMSi differ in charge distribution, molecular geometry, and bonding arrangements at the sulfur center, which may in principle enable their separation by targeted ion chromatography and their differentiation by high-resolution tandem mass spectrometry. Both techniques have been successfully applied to HMS detection (Wei et al., 2020; Campbell et al., 2022), although dedicated method development and experimental verification will be required to confirm the resolution of these two isomers in ambient samples.

Figure 4 Interfacial reaction of HCHO + HOSO${}_{\mathrm{2}}^{-}$ on the ice surface. (a) Snapshot structures of BOMD simulation for the air–quasi-liquid-layer (QLL) interface at T=243 K. (b) Energy decomposition analysis (sobEDAw) plotted against S–C atomic distance (Å) during the reaction between HCHO and HOSO${}_{\mathrm{2}}^{-}$ fragments, where ΔE_els (blue triangles) represents the electrostatic interaction energy, ΔE_xrep (red stars) denotes the exchange repulsion term, ΔE_orb (green circles) indicates the orbital interaction energy and ΔE_disp (yellow triangles) shows the dispersion term. Additional energy decomposition details provided in Sect. S2.

Figure 5 (a) Schematics showing on the hydration behavior of HMS/HMSi on the water surface, where x and y represent the number of interfacial water molecules bound by HMS and HMSi, respectively, via hydrogen bonding. (b) Histograms of probabilities of HMS (upper panel) and HMSi (lower panel) conformations combining different numbers of water molecules. (c) Schematic illustration of the spatial orientation of oxygen atoms (Oa, Ob, Oc, and Od) in HMS and HMSi molecules with respect to the Gibbs dividing surface (GDS), defined as the plane where water density equals half its bulk value. (d) Probability distributions of distinct oxygen atoms across the interface for HMS (upper panel) and HMSi (lower panel), with the yellow zone indicating the interfacial region determined by GDS and water density profiles. Additional methodological details are provided in Supplement.

3.4 Previously Unrecognized Reactions at the Air-Ice Interface

Although air–water interfaces have been widely recognized as critical sites for heterogeneous reactions, the cold conditions prevalent in the upper troposphere and polar regions give rise to an additional interfacial environment in the form of atmospheric ice. These ice surfaces provide a distinct medium for heterogeneous chemistry beyond conventional water vapor and aqueous surfaces, particularly relevant given that HCHO and SO₂ are ubiquitous throughout the atmosphere, including regions where temperatures favor ice formation (Anderson et al., 2017; Bai et al., 2023; Höpfner et al., 2015; Joppe et al., 2024; Sumner and Shepson, 1999). Despite the atmospheric abundance of ice crystals and the widespread distribution of these reactive species, the heterogeneous chemistry occurring at ice surfaces remains poorly characterized compared to reactions in other atmospheric environments. Considering the potential significance of ice-surface heterogeneous reactions for atmospheric sulfur cycling, we investigated the reaction between HCHO and HOSO${}_{\mathrm{2}}^{-}$ at the air-ice interface to elucidate the associated reaction pathway. Given that the freezing process can lead to ion redistribution, we employed classical molecular dynamics (MD) simulations to track HOSO${}_{\mathrm{2}}^{-}$ spatial distribution during ice formation (see Sect. S3 for details). The results showed that HOSO${}_{\mathrm{2}}^{-}$ ions preferentially accumulate at the ice-air interface (Fig. S16), which is likely attributed to their exclusion from the growing ice crystal during crystallization (Ning et al., 2024; Tsironi et al., 2020). Next, we explored the reaction between HCHO and HOSO${}_{\mathrm{2}}^{-}$ at the air-ice interface, with a quasi-liquid layer (QLL) on the ice crystal (Nagata et al., 2019) at T=243 K as shown in Fig. 4a. This temperature represents a typical polar temperature (Zhang et al., 2024a).

Energy decomposition analysis (Lu and Chen, 2023) (Fig. 4b; see Sect. S2 for methodology) revealed that the interaction between HCHO and HOSO${}_{\mathrm{2}}^{-}$ fragments is driven by electrostatic (ΔE_els) and orbital interaction (ΔE_orb) components, with modest contributions from dispersion forces (ΔE_disp). Electronic structure analysis (Fig. S17) confirmed HOSO${}_{\mathrm{2}}^{-}$ as the active electron donor throughout the reaction pathway. These results show that strong electrostatic attraction and orbital interactions jointly facilitate the approach of the two reactants and promote S–C bond formation on the ice surface. Subsequent water-mediated proton transfer yielded HMS. In addition, we investigated HMS formation at temperatures characteristic of the upper troposphere and lower stratosphere. The reaction remains feasible at 223 K, which represents a typical temperature of the high-altitude free troposphere (Li et al., 2024), and HMS formation proceeded rapidly in unbiased BOMD trajectories, completing within ∼ 3 ps through similar mechanistic pathways (Fig. S18). In contrast, HMSi formation was not observed in unbiased trajectories. Metadynamics simulations revealed free energy barriers of 5.4 kcal mol⁻¹ at 243 K and 5.7 kcal mol⁻¹ at 223 K for HMSi formation (Figs. S19 and S20). Additionally, the protonated formaldehyde pathway remained viable on ice surfaces (Figs. S21 and S22). Our simulations revealed three distinct reaction pathways on ice surfaces (Fig. S21): nucleophilic attack at the sulfur atom (Path A, forming HMS), at oxygen atoms O1/O2 (Path B, forming HMSi), and at the hydroxyl group O(H) (Path C, also forming HMSi). Furthermore, when HSO${}_{\mathrm{3}}^{-}$ serves as the sulfur species, the reaction exclusively yields HMSi (Fig. S22).

Our mechanistic findings provide a potential explanation for field observations reporting enhanced HMS concentrations during cold atmospheric periods (Song et al., 2019; Campbell et al., 2024). While HMS formation involves multiple complex atmospheric processes, our results demonstrate that ice-mediated heterogeneous chemistry may significantly enhance this process through some key contributing factors. First, the substantially increased Henry law constants for both HCHO and SO₂ at reduced temperatures facilitate enhanced partitioning of these precursors to condensed phases (Campbell et al., 2022). Second, efficient HMS formation on ice surfaces establishes a previously unrecognized heterogeneous pathway operating across cold atmospheric environments from polar regions to the upper troposphere. This ice-mediated chemistry represents an important yet underappreciated contributor to global organosulfur compound formation and atmospheric sulfur cycling. It should be noted that real atmospheric aerosol and cloud-droplet interfaces are chemically more complex than the simplified system examined here, and they may contain dissolved oxidants, organic species, transition metal ions, and other S(IV) components that could influence interfacial reactivity by competing for reactive sites or modifying the local environment. In future work, we intend to confirm the impacts of other atmospheric components.

Scheme 1 Conceptual overview of four atmospheric scenarios for HMS and HMSi formation (bulk aqueous phase, air–water interface, strongly acidic aerosol, and air–ice interface) and their respective contributions toward reconciling the atmospheric sulfur budget gap. The four labeled panels represent distinct reaction regimes: (1) Bulk-Hindered, (2) Interface-Accelerated, (3) Acid-Catalyzed, and (4) Ice-Mediated.

3.5 Interfacial Dynamics of Products

To further analyze the atmospheric implications of these products, the dynamic interfacial behavior of HMS and HMSi was characterized through quantitative analysis of hydrogen-bonded water molecules at the air-water interface. As shown in Fig. 5, the reaction products are stabilized at the interface through hydrogen bonds (HBs) with surrounding water molecules. Specifically, HMS and HMSi exhibit the highest probability of forming hydrogen bonds with 4–7 interfacial water molecules, exceeding 80 % probability. These products remain bound to the interface through multiple hydrogen-bonding interactions with water molecules, which consequently impedes their transfer into the gas phase. This interfacial retention reflects their strong hydrophilicity, as evidenced by the extensive hydration distribution of HMS and HMSi under different relative humidities (RH = 25 %–100 %) and temperatures (T = 240–298 K; Fig. S23), which further suggests their potential to enhance aerosol hygroscopic growth.

To further elucidate the interfacial behavior of HMS and HMSi, we investigated the positional correlations between their active sites (Oa, Ob, Oc, and Od) and the air–water interface. Figure 5c depicts the distances between these active sites and the Gibbs dividing surface (GDS) to determine their positions at the interface, while Fig. 5d quantifies the probability distributions of each oxygen atom (Oa, Ob, Oc, and Od) across the air-water interface. For HMS, the Ob, Oc, and Od atoms predominantly reside within the interfacial region (yellow shaded area), while the hydroxyl oxygen, designated as Oa, exhibits maximum probability at positive distances from the GDS plane. This indicates preferential orientation of the hydroxyl group toward the gas phase. For HMSi, while the distribution pattern is comparable, the probability distribution curve for Oa shows a significant reduction in extension toward the gas phase compared to HMS. Notably, the spatial distributions of HMS and HMSi at the air-ice interface exhibit similar orientational preferences to those observed at the air-water interface (Fig. S24). This differential orientation of the hydroxyl oxygen may have atmospheric implications, as the greater gas-phase exposure of the hydroxyl group in HMS compared to HMSi suggests a potential enhancement of interfacial oxidation due to orientation preferences, which could facilitate its interaction with atmospheric oxidants such as hydroxyl radicals (OH) and oxygen (O₂). These molecular-level insights into the interfacial behavior of HMS and HMSi provide mechanistic understanding of their atmospheric processing and environmental fate, enabling a more accurate assessment of the contributions of organosulfur compounds to aerosol formation.

4 Conclusions

The overall framework of this study is summarized in Scheme 1. This study clarifies the formation mechanisms of HMS and its isomer HMSi via Born–Oppenheimer molecular dynamics simulations combined with high-level quantum chemical calculations. Our findings provide critical theoretical insights for resolving the persistent atmospheric sulfur budget gap in polluted winter haze episodes and cold environments including polar regions and the upper troposphere.

In the bulk aqueous phase, HMS and HMSi form competitively with near-identical free energy barriers of 7.7 and 7.6 kcal mol⁻¹, respectively. However, bulk aqueous phase contributions to atmospheric sulfur cycling remain limited because aerosols and ice crystals exhibit high surface-to-volume ratios that favor interfacial chemistry. At air–water interfaces in polluted aerosols, HMS formation proceeds with a substantially reduced barrier of 0.6 kcal mol⁻¹ through stepwise water-mediated proton transfer, while HMSi formation faces a higher barrier of 6.1 kcal mol⁻¹, suggesting that interfacial chemistry may substantially enhance HMS formation relative to the bulk aqueous phase. This interfacial selectivity reverses under strongly acidic conditions where pH drops below 2. In these acidic environments, protonated formaldehyde activates unconventional nucleophilic sites on bisulfite, shifting product distribution to favor HMSi with an HMS to HMSi ratio of approximately 1 to 3. At extreme acidity approaching pH 0, HMSi forms exclusively. Based on global surface-layer aerosol pH estimated using annual-mean GEOS-Chem simulations coupled with E-AIM thermodynamic calculations, aerosol pH in the range of −1 to 1 occurs over at least ∼ 20 % of the global surface area where the model successfully converges (Li et al., 2022).

HMS and HMSi are stably retained at air-water interfaces through four to seven hydrogen bonds with surrounding water molecules, occurring with greater than 80 % probability. Their strong hydrophilicity prevents gas-phase evaporation while enhancing aerosol hygroscopic growth. Despite comparable interfacial positioning of oxygen atoms, HMS exhibits a crucial structural distinction: its hydroxyl group orients preferentially toward the gas phase relative to HMSi at both air-water and air-ice interfaces. This enhanced gas-phase exposure may render HMS more accessible to atmospheric oxidants including hydroxyl radicals, potentially facilitating its interfacial oxidation due to orientation preferences. These molecular insights provide a critical basis for assessing organosulfur contributions to aerosol formation and their atmospheric fate.

Collectively, these mechanistic insights provide essential molecular-level understanding for reconciling the persistent discrepancies between modeled and observed sulfate levels. This work establishes a foundation for improving atmospheric sulfur budget models across polluted and cold regions, ultimately enhancing our ability to predict air quality and assess climate impacts. Extending these mechanistic findings to atmospheric models requires further effort, including additional simulations across a wider temperature range and improved characterization of aerosol acidity.