Global modeling of heterogeneous hydroxymethanesulfonate chemistry

Hydroxymethanesulfonate (HMS) has recently been identified as an abundant organosulfur compound in aerosols during winter haze episodes in northern China. It has also been detected in other regions although the concentrations are low. Because of the sparse field measurements, the global significance of HMS and its spatial and seasonal patterns remain unclear. Here, we modify and add to the implementation of HMS chemistry in the GEOS-Chem chemical transport model and conduct multiple global simulations. The model accounts for cloud entrainment and gas–aqueous mass transfer within the rate expressions for heterogeneous sulfur chemistry. Our simulations can generally reproduce quantitative HMS observations from Beijing and show that East Asia has the highest HMS concentration, followed by Europe and North America. The simulated HMS shows a seasonal pattern with higher values in the colder period. Photochemical oxidizing capacity affects the competition of formaldehyde with oxidants (such as ozone and hydrogen peroxide) for sulfur dioxide and is a key factor influencing the seasonality of HMS. The highest average HMS concentration (1–3 μgm−3) and HMS / sulfate molar ratio (0.1–0.2) are found in northern China in winter. The simulations suggest that aqueous clouds act as the major medium for HMS chemistry while aerosol liquid water may play a role if its rate constant for HMS formation is greatly enhanced compared to cloud water.


458
S. Song et al.: Global modeling of heterogeneous hydroxymethanesulfonate chemistry et al., 2017;Brüggemann et al., 2020). Organosulfates are primarily formed by the reactive uptake of gas-phase epoxides on acidic sulfate particles (Froyd et al., 2010;Surratt et al., 2010;Xu et al., 2015). The most abundant organosulfate observed in ambient fine particulate matter (PM 2.5 ) is the isoprene-derived methyltetrol sulfate (C 5 H 11 SO − 7 ), with an average concentration of 1.8 µg m −3 found during August 2015 in Atlanta, Georgia, USA (Hettiyadura et al., 2019). MSA is produced primarily by the oxidation of biogenic dimethyl sulfide (DMS, mainly from marine phytoplankton) and is likely the major organosulfur species in many regions over the oceans (Chen et al., 2018;Hodshire et al., 2019). The concentrations of aerosol-phase MSA in marine environments range from tens to a few hundred ng m −3 (Phinney et al., 2006;Sciare et al., 2009;Huang et al., 2017).
Recently, high mass concentrations of hydroxymethanesulfonate (HMS, CH 2 (OH) SO − 3 ), the most abundant hydroxyalkylsulfonate species commonly found in the atmosphere, have been detected in winter in Beijing, China using an aerosol mass spectrometer by Song et al. (2019a), using an improved ion chromatography method by Ma et al. (2020), and using a UHPLC-LTQ-Orbitrap mass spectrometry by Wei et al. (2020). The mass spectrometry quantification of HMS in ambient aerosols may be subject to the interference of other inorganic and organic sulfur compounds, as suggested by Dovrou et al. (2019). Another difficulty with interpreting HMS in the observational record is the potential for HMS to decompose during sample storage and analysis Moch et al., 2020). The average HMS concentration in the 2015/16 and 2016/17 winters in Beijing was observed to be 1.9 µg m −3 . The highest daily average HMS concentration reached 15 µg m −3 , accounting for 6 % of PM 2.5 concentration . Song et al. (2019a) argued that HMS was likely the major organosulfur compound during winter haze events in northern China. Prior to the two studies, only low levels of HMS, with averages on the order of 0.01 µg m −3 , had been observed in the United States, Japan, and Germany (Dixon and Aasen, 1999;Suzuki et al., 2001;Scheinhardt et al., 2014). More recently, Moch et al. (2020) reported observational evidence for a ubiquitous presence of HMS across 160 locations in North America, Europe, and Asia. Most of the observations of HMS reported by Moch et al. (2020) were unquantified, but mass concentrations of up to 7.6 µg m −3 were reported for Shijiazhuang, China and up to 0.6 µg m −3 for Singapore.
Our knowledge of the chemical mechanism for HMS stems largely from studies in the 1980s when it was recognized as part of the aqueous sulfur chemistry (Pandis and Seinfeld, 1989). Field measurements of cloud water in the Los Angeles Basin showed the coexistence of H 2 O 2 and S(IV) that was much larger than expected based on the phase equilibrium with gaseous SO 2 (Richards et al., 1983). The formation of HMS by the reaction of dissolved SO 2 and HCHO was postulated, and then proved, to explain the observed excess of S(IV) (Munger et al., 1986). The laboratory experi-ments from several groups determined the kinetics and thermodynamics of HMS reactions in aqueous solutions (Boyce and Hoffmann, 1984;Deister et al., 1986;Dong and Dasgupta, 1986;Kok et al., 1986;Olson and Fessenden, 1992). Briefly, both formation and decomposition of HMS depend strongly on pH, i.e., the hydrogen ion activity expressed on a logarithmic scale. HMS is resistant to oxidation by H 2 O 2 and O 3 but reacts with hydroxyl radicals (OH) in the aqueous phase. These studies suggested that the atmospheric conditions favorable for the formation and stability of HMS involved abundant gas-phase SO 2 and HCHO, high aqueous water content, low temperature, intermediate pH, and low photochemical activity.
The integration and reconciliation of data from field observations, laboratory experiments, and chemical modeling are crucial for obtaining a better understanding of how HMS is processed in the atmosphere. Moch et al. (2020) recently implemented in-cloud HMS chemistry into the GEOS-Chem chemical transport model to explore its large-scale spatiotemporal distribution. The HMS chemistry is heterogeneous in nature since the reactions occur in the aqueous phase with reactants transported from the gas phase (Jacob, 2000). Sometimes heterogeneous chemistry is referred to as multiphase chemistry (Ravishankara, 1997). As shown in Fig. 1, the overall heterogeneous reaction rates are controlled not only by rate constants in the aqueous phase but also by mass transfer limitations between the gas and aqueous phases (Jacob, 1986;Ravishankara, 1997;Seinfeld and Pandis, 2016). In partly cloudy conditions, heterogeneous reactions may also be influenced by the entrainment and detrainment of air into and out from clouds (Holmes et al., 2019). Building upon Moch et al. (2020), the current study designed and conducted multiple GEOS-Chem global model simulations, in order to explore the controlling factors and processes of the spatiotemporal distribution of HMS, including cloud entrainment and gas-aqueous mass transfer limitations, kinetics and thermodynamics of HMS chemistry, and the formation media of HMS. The model is driven by the kinetic and thermodynamic data obtained from available laboratory experiments. The simulated results are compared with observations made in four field campaigns. Compared with the control simulation that follows the parameterization in the standard GEOS-Chem model, the default simulation improves treatments of entrainment and mass transfer processes for heterogeneous cloud sulfur chemistry. Both aqueous cloud droplets (Jacob, 1986;Olson and Hoffmann, 1989;Moch et al., 2018;Moch et al., 2020) and aqueous aerosols (Song et al., 2019a;Ma et al., 2020) have been suggested to provide the media for HMS reactions. However, kinetic and thermodynamic data have been determined only in dilute solutions, which are suitable for application in clouds. The lack of corresponding data in concentrated solutions poses a key challenge to modeling HMS chemistry for aerosol water. Therefore, following Moch et al. (2020) we assume that cloud water serves as the only medium in the control and default simulations. The role . Entrainment and detrainment of air into and out from clouds. The volume occupied by aqueous clouds in the grid cell is represented by the cloud fraction (f c ), which is provided by the MERRA-2 meteorological reanalysis in this study. The cloud-free fraction is thus 1−f c . Aqueous aerosols are assumed to be evenly distributed in the grid cell. For aqueous cloud droplets and aqueous aerosols, the same mass transport processes are considered and are shown in the right panel. (Top Right). Gas-phase, interfacial, and aqueous-phase mass transport limitations for the molecules A and B. (Bottom Right). Concentration (C) profiles of A and B are a function of radial distance (r) from the surface of a spherical particle. The subscripts g and aq refer to gas and aqueous phases, respectively. The concentrations are in arbitrary units and their scales are different for gas and aqueous phases. The entrainment/detrainment processes for clouds have been described in detail by Holmes et al. (2019). The right panel is adapted from Fig. 4 in Ravishankara (1997). of aerosol water is explored through sensitivity simulations. Aerosol water chemistry also considers the physiochemical processes in Fig. 1, allowing an evaluation of the importance of the two aqueous media.
This article is organized as follows. In the Methods section, we first provide an overview of the aqueous chemical reactions for HMS, including its formation, decomposition, and oxidation (Sect. 2.1). From existing laboratory studies, we critically estimate the best values and uncertainties of their rate constants. The general configuration of the GEOS-Chem model is described in Sect. 2.2, including its version, simulation period, spatial and temporal resolutions, meteorological field, chemical mechanisms, and underlying emissions. A brief introduction of sulfur simulation in the standard model is given in Sect. 2.3. The two major simulations in this study, control and default, are described in Sect. 2.4 and 2.5, respectively. Based on settings in the standard model, the control simulation implements hetero-geneous HMS chemistry using cloud as the only aqueous medium (Moch et al., 2020). We find that the in-cloud SO 2 titration by various reactants is inappropriately represented in the control simulation, very likely leading to an overestimation of HMS formation. The default simulation fixes this issue. Section 2.6 describes the sensitivity simulations designed to investigate the key factors leading to uncertainty in the modeled HMS levels. In the Results and discussion section, we first show in Sect. 3.1 the spatial and seasonal distributions of HMS from the default simulation and discuss the underlying factors. Differences in the modeled HMS between the default and control simulations are presented and discussed in Sect. 3.2. Section 3.3 demonstrates the key uncertain parameters and processes in the HMS model identified from sensitivity simulations. Sect. 3.4 compares the observations of HMS in four different regions with these model results. The knowledge gained in this study and the remaining gaps are summarized in Sect. 3.5. Finally, the conclusions are given in Sect. 4.

Kinetics and thermodynamics of HMS chemistry
Hydroxymethanesulfonic acid (HMSA, CH 2 (OH) SO 3 H) is a diacid with pK a1 < 0 (Reaction R1) and pK a2 ∼ 12 (Reaction R2). Thus, it primarily exists as HMS (CH 2 (OH) SO − 3 ) in tropospheric clouds and aerosols. In the aqueous phase, HMS is produced by the nucleophilic addition of HSO − 3 and SO 2− 3 to the carbonyl C atom of HCHO (Reactions R3-R6). As SO 2− 3 is a much stronger nucleophile than HSO − 3 , the rate constant of HCHO (aq) + SO 2− 3 , k 2 , is a few orders of magnitude higher than that of HCHO (aq) + HSO − 3 , k 1 , as shown in Table 1. HCHO (aq) refers to the free, unhydrated formaldehyde dissolved in the aqueous phase, and maintains an equilibrium with its hydrated form, CH 2 (OH) 2 (methylene glycol). The equilibrium constant of Reaction (R7), K h , represents the extent of hydration (Eq. 1). Reactions (R1)-(R6) are all reversible and can be summarized by Reaction (R8). SO T 2(aq) is the sum of SO 2 ·H 2 O, HSO − 3 , and SO 2− 3 (Eq. 2). k f (M −1 s −1 ) and k d (s −1 ) represent the forward and backward reaction (HMS formation and decomposition) rate constants of Reaction (R8) and K eq (M −1 ) is its equilibrium constant (Eq. 3). k f is a combination of k 1 and k 2 weighted by the fractions of HSO − 3 and SO 2− 3 in SO T 2(aq) (Eqs. 4-6). K s1 and K s2 denote the first and second dissociation constants for dissolved SO 2 (Table 2). Figure 2 shows the values of k f and k d obtained from the available laboratory experiments as a function of pH (Blackadder and Hinshelwood, 1958;Sørensen and Andersen, 1970;Boyce and Hoffmann, 1984;Deister et al., 1986;Dong and Dasgupta, 1986;Kok et al., 1986;Lagrange et al., 1999). In general, we find a large discrepancy 460 S. Song et al.: Global modeling of heterogeneous hydroxymethanesulfonate chemistry for k f and good agreement for k d .
Therefore, two sets of HMS formation kinetic data can be obtained from Boyce and Hoffmann (1984) and are designated here as the high and low rate constants, as shown in Table 1 and Fig. 2. The high rate constants are the same as those used in Moch et al. (2020). The calculated high and low k f differ by a factor of about 3 at pH < 2 and by a factor of about 6 at pH > 4. The low k f agrees very well (within a factor of 1.1) with the results determined by Kok et al. (1986) and Deister et al. (1986) at higher pH 4, 5, and 5.6 (Fig. 2). The low kinetic data are also closer to the rate constants from the recent quantum chemical calculations by Zhang et al. (2019) (k 1 = 0.9 M −1 s −1 , k 2 = 2 × 10 6 M −1 s −1 , at 298 K). Consequently, the low formation rate constants from Boyce and Hoffmann (1984) are adopted for the default model simulation, while the high ones are used in a sensitivity simulation. Lagrange et al. (1999) proposed another value of k f which was about 1-4 orders of magnitude smaller than the low k f from Boyce and Hoffmann (1984) at pH > 4 (Fig. 2). The simulated HMS concentration is negligible everywhere when applying the k f from Lagrange et al. (1999) in the model, and thus, will not be discussed further.

HMS decomposition
The most complete analysis of K eq was done by Deister et al. (1986). We calculate the expression of k d using the low k f from Boyce and Hoffmann (1984) and K eq from Deister et al. (1986) (Eq. 3 and Table 1). As shown in Fig. 2, k d estimated in this way agrees within a factor of about 2 with results from the other laboratory studies (Blackadder and Hinshelwood, 1958;Sørensen and Andersen, 1970;Dong and Dasgupta, 1986;Kok et al., 1986;Lagrange et al., 1999). Therefore, this expression of k d is adopted in the default simulation, and its value is doubled in a sensitivity simulation. If we use the high k f from Boyce and Hoffmann (1984) and the K eq from Deister et al. (1986), we will obtain a k d that is several times higher than estimates from the other studies. This may serve as circumstantial evidence in favor of the low k f . Moch et al. (2020) uses a value of 3.6 × 10 3 s −1 for k d following Diester et al. (1986).

HMS oxidation
HMS is resistant to oxidation by H 2 O 2 and O 3 but can be oxidized by OH in the aqueous phase (Martin et al., 1989;Olson and Fessenden, 1992). Reaction (R9) produces HCHO and peroxysulfate radical (SO − 5 ) with a rate constant of 2.7× 10 8 M −1 s −1 (Olson and Fessenden, 1992) (Table 1). This value is lower by a factor of about 4 than the results reported in two earlier laboratory studies (Martin et al., 1989;Deister et al., 1990). Olson and Fessenden (1992) argued that these two studies were subject to artifacts and interferences from secondary reactions.
x SO 2 ·H 2 O = SO 2 · H 2 O / SO T

aq
= H + 2 / H + 2 + K s1 H + + K s1 K s2 a EF is obtained by fitting the experimental data shown in Fig. 2C in Liu et al. (2020). µ b is the molality-based ionic strength (mol kg −1 ). b The relationship between k and µ is: Shao et al., 2019). c k 7 is believed to be the lower limit .
is the total dissolved HONO and NO − 2 . e k 9 is determined at 25 • C, µ = 0.5 M. We consider the reaction of HOBr and SO 2− 3 but not the one between HOBr and HSO − 3 , which is included in the standard GEOS-Chem model. The original lab experiments (Liu, 2002) seemed to be interfered by Br 2 , a stronger oxidizing reagent which also reacts with HSO − 3 . A recent study by Liu and Abbatt (2020) suggested that the rate constant of HOBr and HSO − 3 was much lower that of HOBr and SO 2− 3 .
sources and sinks, which are linked to photochemical processes (e.g., photolysis of nitrate and peroxides), transition metal ions (Fenton reactions), and/or reactions with halogen anions and organic matter, are not yet fully understood (Tilgner and Herrmann, 2018 in real aerosols and clouds seem to be 1 order of magnitude lower than modeled concentrations and 1 order of magnitude higher than measured lev- Sander (2015) * The Henry's law constant of HOBr is very uncertain, ranging from 90 to 6000 M atm −1 . HOBr can undergo acid dissociation and has a pKa of 8.65 at 25 • C. We do not consider its acid dissociation because it is only partially dissociated in the interested pH range and because of the high uncertainty of its intrinsic Henry's law constant.  Boyce and Hoffmann (1984) at µ = 1 M. (b) the low k f is also from Boyce and Hoffmann (1984) and corrected for µ and K h . (c) Lagrange et al. (1999): Kok et al. (1986): the reported k f is limited by the dehydration rate of CH 2 (OH) 2 , k dh , and is thus corrected here. (e) is calculated using the k d and K eq determined by Deister et al. (1986) and is also corrected for k dh . The calculated k f values are 2.6 × 10 3 , 2.2 × 10 4 , and 9.1 × 10 4 M −1 s −1 , respectively, at pH = 4, 5, and 5.6 in (d) and (e). For comparison, the extrapolation of the low k f data (b) are 2.7 × 10 3 , 2.4 × 10 4 , and 9.3 × 10 4 M −1 s −1 , respectively, at pH = 4, 5, and 5.6. (g) k d is calculated using the K eq from Deister et al. (1986) and the low k f from Boyce and Hoffmann (1984). (h) Lagrange et al. (1999) Kok et al. (1986) measured k d of 4.8 × 10 −7 and 3.5 × 10 −6 s −1 , respectively, at pH 4 and 5. (m) Dong and Dasgupta (1986) measured K eq at pH 4 and µ = 0.05 M, which translated to a k d of 4 × 10 −7 s −1 . (n) Blackadder and Hinshelwood (1958): k d = 1 × 10 −5 s −1 at pH 5 and µ ≈ 0.1 M. (p) Sørensen and Andersen (1970): k d = 8.5 × 10 −2 s −1 at pH 9 and µ = 0.1 M. For comparison, values of k d calculated by (g) are 5.4 × 10 −7 , 4.9 × 10 −6 , and 4.8 × 10 −2 s −1 at pH 4, 5, and 9, respectively. els (Tilgner and Herrmann, 2018). Since GEOS-Chem does not have a detailed representation of aqueous OH chemistry, we follow Moch et al. (2020) and simply estimate [OH] aq using the modeled [OH] g and a pseudo Henry's law constant H * OH (Eq. 7). In the default simulation, H * OH is set to 4 × 10 −20 M cm 3 molecules −1 . H * OH is more than 1 order of magnitude lower than its intrinsic Henry's law constant, H OH (Table 2), reflecting our presumption that the various organic and inorganic compounds in the aqueous phase act as a net sink for OH radicals. A global mean [OH] g of about 1 × 10 6 molecules cm −3 implies a mean [OH] aq of 4 × 10 −14 M, 1 order of magnitude higher than the mean of the above-mentioned measured [OH] aq . Moch et al. (2020) use a value for H * OH of 1 × 10 −19 M cm 3 molecules −1 based on Jacob et al. (2005).
The products of Reaction (R9) are HCHO (aq) and SO − 5 . Interestingly, the net effect of HMS formation (Reaction R8) and its subsequent oxidation (Reaction R9) is the oxidation of SO T 2(aq) by OH (aq) , which thus represents an indirect oxidation pathway for SO 2 . The sinks for SO − 5 are mainly the reactions with O − 2 , HCOO − , and itself (Reactions R10-R12). The reaction of SO − 5 and HSO − 3 is slow (Jacob et al., 1989). The peroxymonosulfate radical (HSO − 5 ) produced by Reactions (R10) and (R11) can oxidize HSO − 3 to sulfate (Reaction R13) with a similar rate constant to H 2 O 2 + HSO − 3 (Betterton and Hoffmann, 1988). The sulfate radical (SO − 4 ) produced by Reaction (R12) is a very strong oxidant and can react rapidly with HSO − 3 and SO 2− 3 (Reactions R14 and R15) as well as with many other species such as Cl − , NO − 2 , O − 2 , HCOO − , and HO 2 (Jacob, 1986). The rate constants for Reactions (R10)-(R15) can be found in Jacob et al. (1989). It is convenient to define the sulfate yield as the number of SO 2− 4 ions produced due to each attack of OH (aq) on HMS. If SO − 5 reacts with O − 2 /HCOO − (Reactions R10 and R11) and the product HSO − 5 oxidizes HSO − 3 (Reaction R13), the yield is 2. If SO − 5 undergoes self-reaction (Reaction R12) and the produced SO − 4 reacts with HSO − 3 /SO 2− 3 (Reactions R14 and R15), a reaction chain is triggered as the products include SO − 5 . In certain conditions, the sulfate yield can reach 20 or more (Jacob et al., 1989). However, as mentioned above, other oxidizable species also compete for SO − 4 , thereby terminating this chain and leading to a sulfate yield of 1. In remote environments where SO 2 is very low, HSO − 5 may be a stable species, resulting in a sulfate yield < 1. Our low [OH] aq assumption implies the existence of important oxidizable species, and therefore, the chain propagation is limited. As in Moch et al. (2020) the sulfate yield is assumed to be 2 in our simulations.

Phase equilibrium
The gas/aqueous phase equilibriums of HCHO (Reaction R16) and SO 2 (Reaction R17) are described by intrinsic Henry's law constants, H HCHO and H SO 2 , respectively (Table 2). HCHO (aq) is subject to hydration and the apparent Henry's law constant, H * HCHO , is much larger than H HCHO (Eq. 8). SO 2 ·H 2 O dissociates twice in the aqueous phase and thus H * SO 2 depends on pH (Eq. 9). The rates for the hydration of HCHO (aq) (k h in Table 2) and the acid dissociations of SO 2 ·H 2 O (Schwartz and Freiberg, 1981) are fast enough and we assume that these reactions are always in equilibrium.

General model description
We perform global simulations of heterogeneous HMS chemistry using the three-dimensional GEOS-Chem chemical transport model ( The first 6 months are used for initialization and we focus on the 1-year simulation results from September 2015 to August 2016. These months are selected to obtain a continuous boreal winter. We use the tropospheric chemistry mechanism with detailed reactions for O 3 -NO x -VOC (volatile organic compound)-aerosolhalogen interactions. The time step for species advection, vertical mixing, and convection is set to 10 min. The time step is 20 min for emissions, dry deposition, photolysis, and chemistry, as recommended by Philip et al. (2016). The simulated aerosol species include secondary inorganic (sulfate, nitrate, and ammonium) and organic aerosols, primary organic aerosols, black carbon, dust, and sea salt.
Because of the importance of acidity for heterogeneous HMS chemistry, more details are provided for the calculation of cloud water and aerosol pH. The standard model calculates cloud water pH iteratively with an initial estimate of 4.5, as described in Alexander et al. (2012)

Sulfur simulation in the standard model
The sulfur simulation in GEOS-Chem has been developed and improved based on multiple studies (Chin et al., 2000;Park et al., 2004;Alexander et al., 2005Alexander et al., , 2009Chen et al., 2017;Shao et al., 2019). The simulated sulfur species include DMS, SO 2 , MSA, and sulfate. It includes primary emissions of DMS, SO 2 , and sulfate (Sect. 2.2). SO 2 , MSA, and sulfate can also be formed by chemical reactions. The model contains three gas-phase reactions of DMS oxidation, producing SO 2 and MSA (Reactions R18-R20). An expanded chemistry mechanism for DMS can be found in Chen et al. (2018). The oxidation of SO 2 to sulfate occurs in the gas phase by OH (Reaction R21) and in the aqueous clouds. The aqueousphase oxidants are O 3 , H 2 O 2 , O 2 (catalyzed by transition metal ions Mn 2+ and Fe 3+ ), and HOBr (Reactions R22-R25). The effect of the heterogeneity in cloud droplet pH on sulfate production rates is accounted for using the parameterization by Yuen et al. (1996) and Fahey and Pandis (2001). This parameterization is restricted over the ocean since the heterogeneity in pH is believed to be caused by alkaline seasalt aerosols (Alexander et al., 2012). The model also includes the oxidation of SO 2 by O 3 on sea-salt aerosol surface (Reaction R26) (Alexander et al., 2005).

Control simulation
Based on the standard model v12.1.0, we implement heterogeneous HMS chemistry and assume that cloud water provides the only aqueous medium, following Moch et al. (2020). As described in Sect. 2.1, HMS is produced by dissolved SO 2 and HCHO, undergoes decomposition, and is oxidized to sulfate by aqueous OH. Two other cloud sulfate formation pathways are also incorporated , in which SO 2 is oxidized by NO 2 and HONO (Reac-tions R27 and R28).
Tables 1 and 2 show all the rate constants of aqueousphase reactions and the Henry's law constants of the reactants. The solubilities of transition metals Fe and Mn are reduced following Shao et al. (2019). Ten advected tracers are added: one is the aerosol HMS species and the others represent different sulfate formation pathways. Transport and deposition of these tracers are treated in the same way as the sulfate tracer as in Moch et al. (2020). In addition, several other changes are made in the control simulation to the standard model. First, we update the dry deposition scheme and the reactive uptake coefficients of NO 2 , NO 3 , and N 2 O 5 on aerosols, following Jaeglé et al. (2018) and Shah et al. (2018). Second, this simulation includes some updates developed by Luo et al. (2019Luo et al. ( , 2020 in the treatments of wet processes, including spatially and temporally varying in-cloud condensation water contents, empirical washout rates for water-soluble aerosols and nitric acid, the cloud fraction available for aqueous chemistry, and rainout efficiencies for water-soluble aerosols and gases. Third, more ions are included in the cloud water pH calculation. We consider Ca 2+ , Mg 2+ , NH + 4 , Na The Newton-Raphson method is used to find the solution to the cubic electroneutrality equation, following Luo et al. (2020) and Moch et al. (2020). Ca 2+ and Mg 2+ are assumed to constitute 3 % and 0.6 %, respectively, of the dust by mass (Claquin et al., 1999;Fairlie et al., 2010;Nickovic et al., 2012;Shao et al., 2019;Moch et al., 2020). Only Na + and Cl − from sea-salt aerosols are considered. HMS and CH 3 SO − 3 are from the cloud scavenging of aerosols. NO − 2 , HCOO − , and CH 3 COO − are from the scavenging of HONO, HCOOH, and CH 3 COOH, respectively. Fourth, HMS, CH 3 SO − 3 , and Ca 2+ and Mg 2+ in fine dust are included in the ISORROPIA calculations. We assume the same hygroscopicity of HMS and MSA as sulfate .
We evaluate the performance of the control simulation by comparing it with the standard GEOS-Chem v12.1.0 (GC12.1.0). Figure S1 in the Supplement shows the horizontal distributions of surface SO 2− 4 and SO 2 concentrations. The global average SO 2− 4 in the control simulation is reduced by 24 % compared to GC12.1.0. The updates in the treatments of wet processes by Luo et al. (2019Luo et al. ( , 2020 are primarily responsible for this difference. The SO 2− 4 concentrations modeled in the control simulation are consistent with the improved model results in Luo et al. (2020), which have been found to agree well with SO 2− 4 observed in the United States, Europe, and Asia (Luo et al., 2020). Moreover, since GC12.1.0 was released in late 2018, it is necessary to compare it with a more recent model ver-sion. Accordingly, we conduct a simulation using the standard GEOS-Chem v12.7.0 (GC12.7.0, released in February 2020, http://wiki.seas.harvard.edu/geos-chem/index.php/ GEOS-Chem_12, last access: 10 June 2020). We find that the global average SO 2− 4 in GC12.7.0 only differs little (3 %) compared with that in GC12.1.0 (Fig. S2 in the Supplement).
Below, we provide details on the calculation of cloud sulfur chemistry and highlight the need for more accurate representations of in-cloud SO 2 titration by various reactants, which include O 3 (Reaction R22), H 2 O 2 (Reaction R23), O 2 (Reaction R24), HOBr (Reaction R25), NO 2 (Reaction R27), HONO (Reaction R28), and HCHO (Reaction R8). Cloud sulfur chemistry is calculated locally in the model grid cells where aqueous clouds are present. f c (dimensionless, ≤ f c ≤ 1) denotes the fraction of aqueous cloud in a grid cell, and L (m 3 H 2 O m −3 air) denotes the in-cloud liquid water content. In each chemistry time step ( t = 20 min), the losses of SO 2 in the above reactions (R8, R22-R25, R27, and R28) are calculated. Reaction (R24) is treated as a first-order reaction of SO 2 (O 2 is in large excess), while the other reactions are second order. The first-and second-order rate constants for the aqueous reaction of SO T 2(aq) and X i (aq) , k 1,aq,i (s −1 ) and k 2,aq,i (M −1 s −1 ), are obtained by Eq. (10) from the kinetic data in Table 1. X i (i = 1 : 7) represents the ith reactant with SO 2 . R aq,i is the reaction rate (M s −1 ). k 1,aq,i and k 2,aq,i are used to derive the first-and second-order rate constants for the heterogeneous reaction of SO 2(g) and X i (g) , k 1,g,i (s −1 ) and k 2,g,i (mol mol −1 s −1 ) (Eq. 11). H * SO 2 and H X i indicate their Henry's law constants. f g,SO 2 and f g,Xi are the gas-phase partitioning fractions of SO 2 and X i , respectively (Eq. 12). R is the gas constant. T (K) is the temperature. P (atm) is atmospheric pressure. The loss of SO 2 over time t, SO 2g,i , is solved analytically (Eq. 13). SO 2,t=0 g and [X i t=0 ] g are the mixing ratios (mol mol −1 ) for SO 2(g) and X i (g) at the beginning of this time step. The grid-average losses of SO 2(g) from all seven reactions are limited by the availability of SO 2(g) within the cloud fraction f c . k 1,aq,i =R aq,i / SO T 2 aq and k 2,aq,i k 1,g,i =k 1,aq,i H * SO 2 f g,SO 2 LRT and k 2,g,i = k 2,aq,i H * Since multiple in-cloud reactions consume SO 2 simultaneously, it is important to allow them to compete effectively and fairly. As shown in Eq. (13), the contribution of the ith reaction to the total SO 2(g) loss depends on its rate constant (k 1,g,i or k 2,g,i ), its relative abundance ([X i t=0 ] g / SO 2,t=0 g ), and the choice of t. Ideally, t should be smaller than the lifetime (τ i ) of SO 2(g) for any ith reaction. τ i is the inverse of the pseudo-first-order rate constant, k ∼ 1,g,i , which equals to k 1,g,i for a first-order reaction and to k 2,g,i [X i t=0 ] g for a secondorder reaction. Figure 3 shows the probability density distributions of the calculated k ∼ 1,g,i and the total rate constant for the seven reactions, 7 i=1 k ∼ 1,g,i , in the lower troposphere for a randomly selected week in boreal summer. k ∼ 1,g,i (and thus τ i ) can vary by several orders of magnitude in different model grid cells. Notably, there is a > 50 % possibility that the lifetime of SO 2(g) is smaller than 20 min, the t used in this simulation. The rapid consumption of SO 2(g) is mainly via O 3 and H 2 O 2 , as shown in Fig. 3 and Table S1 in the Supplement (statistics of probability distributions). The other five reactions consuming SO 2 (HCHO, TMI+O 2 , NO 2 , HONO, and HOBr) can be considered as minor pathways. This means that using t = 20 min for the sulfur chemistry will in general lead to an underestimation of the contribution of O 3 and H 2 O 2 and an overestimation of the importance of the other reactants such as HCHO. An example is provided in Text S1 in the Supplement to conceptually explain the effect of t on the competition of different reactions. We conduct a sensitivity simulation in which t is set to 10 min and, as we expect, the SO 2− 4 concentrations through the cloud O 3 chemistry increase significantly (Fig. S3 in the Supplement). A simple way to solve this problem is to reduce t. The possibility of τ < t decreases to only 4 % when t = 1 min ( Fig. 3 and Table S1). Also, most (> 80 %) of the cases of τ < 1 min arise from the rapid reaction of SO 2− 3 with O 3(aq) when cloud water pH is high. The remaining cases are from the reactions of SO 2− 3 with HOBr (aq) and HCHO (aq) . The other four reactions can hardly lead to τ < 1 min. We change the time step to 1 min when calculating in-cloud SO 2 titration in the default simulation (Sect. 2.5).
Another issue in the control simulation is, in a partly cloudy (< f c < 1) model grid, that the mixing of air between the cloudy fraction (f c ) and the cloud-free fraction (1−f c ) occurs in the same timescale as the chemistry time step of the model (Holmes et al., 2019). In each time step, the gridaverage loss of SO 2(g) from all in-cloud reactions cannot exceed the amount of SO 2(g) available within the cloudy fraction and at the beginning of this time step (Eq. 14). This socalled "cloud-partitioning method" is unphysical as the entrainment/detrainment rates are affected by the setting of the chemistry time step (Holmes et al., 2019). Since many chemical transport models such as GEOS-Chem do not resolve individual clouds, Holmes et al. (2019) developed a more realistic and stable "entrainment-limited uptake" method, which . Probability density distributions of the pseudo-first-order rate constants with respect to SO 2(g) for cloud reactions in the control simulation. The shaded area shows the sum of rate constants for the seven reactions consuming SO 2 . The red, green, blue, and orange curves indicate the distributions for reactions with O 3 , H 2 O 2 , HCHO, and the sum of rate constants for the other four reactions consuming SO 2 (TMI+O 2 , NO 2 , HONO, and HOBr), respectively. Data shown are for the first week of July and in the lower troposphere (13 vertical layers above surface up to about 800 hPa). Since the chemistry time step ( t) of this simulation is 20 min, there are 504 steps in this week. The total number of data points is 72 × 46 (number of 5 • × 4 • grids) × 13 (vertical layers) × 504 ≈ 2.2× 10 7 . About 1.4 × 10 7 data points have aqueous clouds (cloud fraction f c > 0), accounting for about 2/3. The probability density distributions are plotted based on these data points. The dashed and solid vertical black lines indicate the rate constants corresponding to t of 20 and 1 min, respectively. accounts for cloud entrainment/detrainment within the chemical rate expression. We apply this method to the default simulation (Sect. 2.5).

Default simulation
Three major changes are made in this simulation based on the control simulation, as mentioned in Sect. 2.4. The first is applying the entrainment-limited uptake method developed by Holmes et al. (2019) to more realistically model the entrainments and detrainments of air in cloudy grid cells. The second is reducing the time step to 1 min when calculating cloud sulfur reactions to better quantify the competition of different chemical pathways consuming SO 2 . The third is adding the reaction of H 2 O 2 and SO 2 in aerosol water using the kinetic data reported recently by Liu et al. (2020). The control simulation only includes the reaction of H 2 O 2 and SO 2 in cloud water. Figure S4 in the Supplement shows the horizontal distributions of surface SO 2− 4 and SO 2 concentrations in the control and default simulations, and only very small differences (4 % for SO 2− 4 and 1 % for SO 2 ) are found for their global average values.
In the entrainment-limited uptake method (Holmes et al., 2019), the first-order loss rate of SO 2(g) in a model grid cell due to heterogeneous cloud chemistry, k 1 (s −1 ), depends on the cloud fraction (f c ), the detrainment rate (k c , s −1 ), and the in-cloud total pseudo-first-order rate constant, k * 1,g = 7 i=1 k * 1,g,i (s −1 ) (Eq. 15). As shown in Holmes et al. (2019), the entrainment/detrainment (k c term) limits its reactive uptake. In a completely cloudy condition (f c = 1), Eq. (15) reduces to k 1 = 7 i=1 k * 1,g,i . k c is the reverse of the in-cloud residence time of air (τ c ), which varies with cloud types and ranges from 15 to 120 min for stratus and cumulus clouds (Holmes et al., 2019). We use τ c = 30 min in this work since MERRA-2 does not provide this information. A sensitivity simulation shows that assuming a τ c of 60 min decreases the global average surface SO 2− 4 concentration by 10 %. Holmes et al. (2019) have pointed out that future studies are needed to specify the spatiotemporal variability of τ c in the global reanalysis data sets. Within the cloudy fraction of a model grid cell, as shown in Fig. 1 and Eq. (16), the heterogeneous reaction rates are limited by a series of resistances associated with the mass transfer processes from the gas phase to the aqueous phase, including gas-phase diffusion, transfer of the reactants across the air-water interface, and aqueous-phase diffusion (Ravishankara, 1997;Jacob, 2000). In Eq. (16), the α SO 2 term represents the limitation due to mass accommodation at the air-water interface and the D g,SO 2 term represents that due to gas-phase diffusion. A dimensionless parameter Q, whose expression is given by Eq. (17) (0 < Q < 1), is used to account for aqueous-phase mass transport limitations when calculating k 1,aq,i and k 2,aq,i (Eq. 10).
Here, r is the radius of cloud droplets and is assumed to be 10 −5 m. A (m 2 m −3 air) is the surface area density of cloud droplets and is derived using L and r. v SO 2 (m s −1 ) is the molecular mean speed of SO 2 (Eq. 18). α SO 2 (dimensionless) is the mass accommodation coefficient of SO 2 ( Table 3). D g,SO 2 (m 2 s −1 ) is the gas-phase diffusion coef-ficient of SO 2 (Eq. 19). q is a dimensionless parameter determined by r, D aq , and k ∼ 1,aq,i (the pseudo-first-order rate constant with respect to SO T 2(aq) for the ith reaction). For a first-order and second-order reaction, k ∼ 1,aq,i is equal to k 1,aq,i and k 2,aq,i [X i ] aq , respectively (Eq. 10). D aq is the aqueous-phase diffusion coefficient (10 −9 m 2 s −1 ) (Song et al., 2019a). M SO 2 (g mol −1 ) represents the molar mass of SO 2 . ρ n,air (molecule cm −3 ) is the number density of air. k ∼ 1,g,i (s −1 ) is the pseudo-first-order rate constant with respect to SO 2(g) for the ith aqueous-phase reaction, and is equal to k 1,g,i and k 2,g,i [X i ] g for the first-order and second-order reactions, respectively. In addition, as illustrated in Fig. 1, the second-order reaction rate may also be limited by the mass transfer of X i . Thus, the in-cloud pseudo-second-order rate constant, k * 2,g,i , is given by Eq. (20). v X i , α X i , and D g,X i are the molecular mean speed, the mass accommodation coefficient, and the gas-phase diffusion coefficient of X i , respectively. v X i and D g,X i are calculated similarly to Eqs. (18) and (19). α X i can be found in Table 3.
As mentioned in Sect. 2.4, we do not change the chemistry time step of the model ( t = 20 min) but only the time step (to 1 min) when identifying cloud SO 2 reactions. For each 1 min time step, the loss of SO 2(g) for the ith reaction, SO 2g,i , is solved analytically using Eq. (13), in which k 1,g,i and k 2,g,i are replaced by k * 1,g,i from Eq. (16) and k * 2,g,i from Eq. (20), respectively. This change reflects the mass transport limitations. The grid-average first-order loss rate of SO 2(g) , k 1 , is calculated using Eq. (15), in which is replaced by the in-cloud total pseudo-first-order rate con- SO 2g,i . The grid-average loss of SO 2(g) and the contributions of different reactions are then calculated using k 1 and SO 2g,i . The mixing ratios of the relevant chemical species are updated at the end of this 1min time step and used as initial condition for the next time step. The calculations are repeated 20 times in a chemistry time step. The cloud water pH and rate constants for the heterogeneous reaction of SO 2(g) and X i (g) are calculated only at the beginning of each chemistry time step. We conduct a sensitivity simulation that redoes cloud SO 2 calculations using  the cloud water pH and rate constants estimated at the end of each chemistry time step. The resulting change is insignificant (global mean SO 2− 4 concentration decreases by < 2 %). The aqueous-phase sulfur reactions are hard-coded into the model. Ideally, further model development of cloud chemistry should apply the advanced numerical solvers generated by the Kinetic PreProcessor (KPP), which may not only allow a full coupling of gas-phase and cloud chemistry but also make it easier for the model to incorporate additional aqueous reactions (Fahey et al., 2017;Viral Shah, personal communication, 18 December 2019).
The implementation of sulfur chemistry in aerosol water is similar to that for cloud sulfur chemistry. As shown in Fig. 1, the heterogeneous reaction rates are also controlled by the mass transfer of reactants from the gas to the aqueous phase. The difference is that aerosols (and aerosol water) can be considered evenly distributed in a model grid cell, and it is unnecessary to include the entrainment/detrainment processes. The major difficulty in parameterizing the aerosol water sulfur chemistry is the lack of suitable reaction rate constants. Liu et al. (2020) have recently found that the high ionic strength of deliquesced aerosols significantly enhances the rate constants for the reaction of H 2 O 2 and SO 2 . The enhancement factor (EF) relative to its rate constant in dilute solutions is derived by fitting the data in Liu et al. (2020) as a function of the molality-based ionic strength, µ b (mol kg −1 ) ( Table 1). The water content, pH, µ b , and the absorbed water volume fraction of inorganic aerosols are calculated by the ISORROPIA II model (Sect. 2.2). The aerosol water volume fraction of 0.25 is used as a threshold for the occurrence of aqueous reactions as it governs the transition of aerosols to a liquid state (Bateman et al., 2016).

Sensitivity simulations
In addition to the control and default simulations, we conduct ten sensitivity simulations to investigate the key factors leading to uncertainty in the modeled HMS concentrations. As shown in Table 4, all these sensitivity simulations are based on the default simulation. HiKf, HiKd, and HiOH make changes to HMS formation, decomposition, and oxidation, respectively, in heterogeneous cloud chemistry. HiKf uses the high k f instead of the low k f in the default simulation (Sect. 2.1.1). HiKd increases k d by a factor of 2, the upper limit of its estimate (Sect. 2.1.2). HiOH increases [OH] aq by a factor of 10, matching its average value in current multiphase models (Sect. 2.1.3). CWpH considers less ions, i.e., NH + 4 , H + , OH − , SO 2− 4 , NO − 3 , HSO − 3 , SO 2− 3 , HCO − 3 , CO 2− 3 , HMS, and CH 3 SO − 3 , in cloud water pH calculations. AWOH, AWK0, and AWKE examine the potential role of aerosol water in heterogeneous HMS chemistry (Table 4). Since the rate constants of HMS chemical reactions in concentrated solutions have not been determined experimentally, we have to make assumptions about these data. AWOH implements the oxidation of HMS by OH in aerosol water and assumes the same rate constant as those for cloud water. AWK0 adds the formation and decomposition of HMS in aerosol water also using the rate constants for cloud water. Theoretically, we anticipate that the rate constant of HMS formation, k f , in concentrated solutions should be enhanced relative to dilute solutions (Song et al., 2019a), similar to the situation found for the reaction of H 2 O 2 and SO 2 . AWKE arbitrarily increases k f by the same EF for the H 2 O 2 and SO 2 reaction ( Table 1). The implementation of the above chemical reactions of HMS in aerosol water follows the approach described in Sect. 2.5.
Three sensitivity simulations, HiNH3, HiFA, and AppHet, focus on East Asia (Table 4). SO 2 , HCHO, and NH 3 emissions may influence the modeled HMS. Although SO 2 emissions are well understood, recent studies have shown that there may be large uncertainties in emissions of NH 3 and HCHO in China (Pan et al., 2018;Kong et al., 2019;Song et al., 2019a). An inverse study found that the MEIC inventory underestimated NH 3 emissions by 30 % nationally and by > 40 % in eastern and central regions using observations over the same time period as our study (Kong et al., 2019). HiNH3 increases the anthropogenic emissions of NH 3 in MEIC by 50 %. Less information is available regarding the emissions of HCHO due to its sparse observations and complex chemistry. Model-observation comparisons in Beijing suggested a strong underestimation of HCHO emis-sions during winter (Song et al., 2019a). Mobile and residential emission sources may be responsible for its underestimation Song et al., 2019a). HiFA increases HCHO emissions from the transportation and residential sectors by a factor of 5. Chemical transport models commonly underestimate SO 2− 4 during winter haze episodes in China (Chu et al., 2020), and thus some studies have adopted an apparent heterogeneous parameterization for SO 2 reactive uptake in order to compensate for the missing SO 2− 4 (Zheng et al., 2015;Cheng et al., 2016;Chu et al., 2016;Wang et al., 2016;Liu et al., 2019;Li et al., 2020). This parameterization is applied in AppHet during the cold season, in which the reactive uptake coefficient of SO 2 increases from 2 × 10 −5 to 5 × 10 −5 with 50 % < RH ≤ 100 % (Zheng et al., 2015). The surface HMS concentrations and HMS / sulfate molar ratios exhibit distinct seasonal patterns with maxima in DJF (boreal winter) and minima in JJA (boreal summer), and thus our analyses focus on these two seasons. It is noted that there are hotspots of HMS in JJA in Siberia that are linked to massive forest fires in that region in July 2016 (Sitnov et al., 2017). As highlighted in Fig. 4, three regions with relatively high HMS levels, East Asia (EA), Europe (EU), and North America (NA), are selected for quantitative analysis.  3 µg sm −3 ) and HMS / sulfate ratios (0.1-0.2) are found during the winter season in northern China (Fig. S7).
As mentioned in Sect. 1, previous studies have suggested that the formation and existence of HMS in the condensed phase are generally favored under the following conditions: high precursor (SO 2 and HCHO) concentrations, low photochemical oxidant levels, low temperature, abundant aqueous water, and moderate acidity (Munger et al., 1984(Munger et al., , 1986Moch et al., 2018;Song et al., 2019a;Ma et al., 2020). The seasonal variability of HMS does not follow that of the precursor levels ( Fig. S6 in the Supplement). The seasonal variation of the geometric mean of the two precursors, √ SO 2 × HCHO, is weak because of their opposite seasonality (more SO 2 but less HCHO in winter). The cloud liquid water content (L) in the lower troposphere shows a spatial distribution with higher values over the ocean and lower values over land (Fig. S8). There is no consistent seasonal pattern of L between DJF and JJA over the three regions (EA, EU, and NA). The modeled cloud water pH exhibits a seasonal difference. The average pH in DJF (JJA) is 4.3 (5.8), 4.7 (5.6), and 4.7 (5.7) for EA, EU, and NA, respectively. The higher pH in JJA is related to more abundant gaseous NH 3 (Fig. S8 in the Supplement), given the buffer capacity of NH 3 in moderating the acidity of atmospheric condensed water (Song et al., 2019b).
One of the above-mentioned factors favoring HMS is the moderate acidity. This term is somewhat ambiguous but is used to represent the pH range that allows for relatively rapid formation and slow decomposition of HMS. We show in Fig. 2 that both k f and k d increase with pH. The lifetime of HMS with respect to decomposition is about 60, 6, and 0.6 hours at pH 5, 6, and 7, respectively, at 298 K, and is even larger at lower T . For the range of pH in the three regions (its average from 4.3 to 5.8), the decomposition of HMS is so slow that its chemical equilibrium is difficult to achieve. Accordingly, the modeled HMS levels are predominantly controlled by formation kinetics. This is supported by the results from the HiKd simulation, in which k d × 2 makes little difference in the modeled HMS compared to the default simulation (Fig. 5). The higher cloud pH in JJA should lead to faster HMS formation rates than those in DJF. However, in the default simulation, the modeled HMS levels show an opposite pattern. This is believed to be linked to the different photochemical oxidizing abilities in the two seasons. Globally, the two main aqueous oxidants for SO 2 are O 3 and H 2 O 2 , which compete with HCHO. The competition of different pathways can be influenced by the levels of these gases, T (changing gas solubilities and rate constants), and pH (changing rate constants). O 3 , H 2 O 2 , and HCHO all have higher concentrations in JJA (Fig. S6). The lower T in DJF favors the H 2 O 2 reaction most and the O 3 reaction least. Notably, the response of the O 3 reaction to pH is essentially the same as that for HCHO since both react rapidly with SO 2− 3 . The HCHO + SO 2 reaction is significant only when the two photochemical oxidants are inefficient.

Difference between the default and control simulations
Compared with the default simulation, the control simulation realizes very different spatial and seasonal distributions of HMS concentrations and HMS / sulfate molar ratios (Fig. S9 in the Supplement). Figure 6 shows the differences in surface HMS concentrations for DJF and JJA. The corresponding information for MAM and SON is presented in Fig. S10 in the Supplement. Two features are evident. First, a weak seasonality is found for the control simulation, but for the default simulation, HMS is much more abundant in DJF. Second, the control simulation predicts significantly higher HMS concentrations almost everywhere except in parts of East Asia and Europe in DJF. Interestingly, the only region where the default simulation gives higher HMS concentrations is wintertime in northern China, the focus of several studies (Moch et al., 2018;Song et al., 2019a;Ma et al., 2020). Specifically, as shown in Fig. 5, the average HMS concentrations modeled by the control simulation in DJF (JJA) are 0.60 (0.70), 0.22 (0.12), 0.14 (0.11) µg sm −3 , respectively, for East Asia (EA), Europe (EU), and North America (NA). The average HMS / sulfate ratios inferred from the control simulation in DJF (JJA) are 0.09 (0.14), 0.08 (0.06), 0.09 (0.07), respectively, for EA, EU, and NA. As described in Sect. 2.5, based on the control, the default simulation improves the representations of heterogeneous cloud sulfur chemistry in the model by applying the entrainment-limited uptake method from Holmes et al. (2019) and by reducing the time step when calculating aqueous sulfur reactions. These changes allow for a more realistic simulation of entrainments and detrainments of air in partly cloudy grid cells and for an effective competition of different aqueous reactions consuming SO 2 . In the control simulation, the time step for calculating in-cloud sulfur reactions is the same as the chemistry time step of the model, t = 20 min. But there may be a > 50 % possibility that the lifetime of in-cloud SO 2 is less than this t, as shown in Fig. 3 and described in Sect. 2.4. Given that the main reactants with in-cloud SO 2 are O 3 and H 2 O 2 , this setting leads to a general underestimation of the contribution of O 3 and H 2 O 2 and an overestimation of importance of the minor reactants such as HCHO. The bias is larger in JJA than in DJF, as suggested by the probability distribution statistics for the in-cloud lifetime of SO 2 (Table S1 in the Supplement).

Key uncertain factors
Section 2.6 and Table 4 describe the ten sensitivity simulations we conduct with an aim to find out the key parameters and processes leading to HMS modeling uncertainties. All of these simulations are modified based on the default simulation and can be classified into three groups: heterogeneous cloud chemistry (HiKf, HiKd, HiOH, and CWpH); heterogeneous aerosol water chemistry (AWOH, AWK0, and AWKE); and East Asia only (HiNH3, HiFA, and AppHet). A comparison of the surface HMS concentrations and HMS / sulfate molar ratios from these sensitivity simulations is provided in Fig. 5, focusing on three regions (EA, EU, and NA) and two seasons (DJF and JJA).
First, we examine the sensitivity simulations in terms of the formulation of heterogeneous cloud chemistry. HiKf, HiKd, HiOH, and CWpH make changes in HMS formation, decomposition, oxidation, and cloud water pH calculations, respectively. The surface HMS concentrations and HMS / sulfate ratios in the latter three indicate relative differences of less than ±20 % compared to the default simulation. However, HiKf shows a very large increase, by a factor of 2 to 6, in modeled HMS. This is expected since the high and low k f differ by a factor of about 3 at pH < 2 and by a factor of about 6 at pH > 4 (Sect. 2.1.1). As described in Sect. 2.1.3, the formation of HMS and its oxidation by OH represent an indirect oxidation pathway for SO 2 . The sulfate yield, defined as the number of SO 2− 4 ions produced due to each attack of OH on HMS, is set to 2 in our simulations. The small difference between the HiOH and default simulations suggests that this indirect pathway should be insignificant.
Second, three sensitivity simulations are conducted for East Asia, as it is found most suitable for the existence of HMS. HiNH3 and HiFA increase the concentrations for modeled HMS in DJF by 60 % and 20 %, respectively, whereas the concentrations by AppHet are decreased by about 30 %. The changes due to HiNH3 and HiFA are much smaller in JJA (Fig. 5). The increase of HMS in HiNH3 can be attributed to higher cloud water pH, and its decrease in AppHet should be related to a decrease in SO 2 available for cloud chemistry. Interestingly, HiNH3 increases the HMS / sulfate ratios in DJF by only 20 %. The higher cloud water pH enhances the formation of sulfate through the pH-sensitive pathways such as the reaction of SO 2 with O 3 .
Third, AWOH, AWK0, and AWKE explore the potential role of heterogeneous aerosol water HMS chemistry. The challenge of modeling aerosol water HMS chemistry is the lack of its reaction rate constants in concentrated aqueous solutions. We use the rate constants from dilute solutions in AWOH and AWK0. The oxidation of HMS by OH in aerosol water leads to losses of 10 %-20 % (DJF) and 40 %-60 % (JJA) (Fig. 5). The formation and decomposition of HMS in aerosol water result in negligible changes in the modeled HMS concentrations, as shown in Fig. 7 (DJF) and Fig. S11 (the other seasons). Results from AWKE suggest that aerosol water might play a role in the formation of HMS only when the k f is strongly enhanced in concentrated solutions like the rate constant of the SO 2 reaction with H 2 O 2 .
Overall, our sensitivity simulations suggest that the key uncertain parameter in the model is k f . Based on existing experimental results, the low value for k f seems more reason- able (Sect. 2.1.1). The key uncertain process is modeling the aerosol water chemistry of HMS in the absence of reliably defined rate constants.

Comparison with prior observations
The quantitative observations of HMS in ambient aerosols remain sparse (Moch et al., 2020) and we provide here a comparison between the observations at four sites and two model simulations (control and default) in Table 5. Since these observations have been collected over the past three decades while our simulations cover only one year, it is more appropriate to use the molar ratios of HMS / SO 2− 4 or HMS/MSA rather than absolute HMS concentrations. The measurement techniques of HMS in different studies were different (Table 5) and their effects on the reported results are unclear. Moreover, observations of HMS may be subject to measurement artifacts that may make quantitative comparisons between model results and observations difficult Moch et al., 2020). Among the observations shown in Table 5, the highest HMS / SO 2− 4 ratio of 11 % has been found in winter in Beijing by . Model results from both default and control simulations agree well with this observed ratio. Less HMS was observed in other regions, including New Mexico (USA), Germany, and Osaka Bay (Japan). The default simulation overestimates the HMS / SO 2− 4 or HMS/MSA ratios by a factor of 2-3, whereas the control simulation overestimates these ratios by an order of magnitude.
A more detailed comparison of the model with observations in Beijing is provided below . These samples from Beijing were stored at −20 • C between sample collection and analysis, had HCHO added to solution during sample extraction, and were examined by ion chromatography using a more acidic eluent than normal all in order to limit HMS decomposition and misidentification. The observations in Ma et al. (2020) cover 73 d in winter and 11 polluted days in other seasons. The data for the other seasons is presented only in their discussion paper. Because of the coarse resolution of global model, we do not expect our simulations to capture the day-to-day variability that is ob-served at a single site. Accordingly, we examine the ability of our simulations to reproduce the observed relationships between HMS and its influencing factors. Figure 8 provides scatter plots of HMS concentrations (and HMS / SO 2− 4 ratios) versus two variables (O 3 and RH) and compares the data from observations and model simulations (control and default). The level of O 3 represents photochemical oxidizing capacity and RH may indicate the abundance of aqueous water in the lower troposphere.
We find a similar relationship between HMS and O 3 from the observations and default simulation (Fig. 8). Significant HMS levels are observed and modeled only under low O 3 conditions (< 20 ppb). However, the control simulation obtains another cluster of days with high HMS levels when O 3 is abundant (> 40 ppb). This cluster is linked to the inappropriate representation of heterogeneous cloud sulfur chemistry in the control simulation. The large time step for SO 2 titration excessively favors the reaction of SO 2 with HCHO, as described in Sect. 2.4 and Fig. 3. It should be noted in Ma et al. (2020) that only 11 daily samples had O 3 levels larger than 20 ppb. There might be a possibility that the days with high levels of both HMS and O 3 were missed in their sampling coverage, but we think it is unlikely given the rapid oxidation of SO 2 by the photochemical oxidants (gaseous OH and aqueous O 3 and H 2 O 2 ) under such conditions.
The scatter plots of HMS and RH show a similar exponential-like relationship in the observations and model simulations (Fig. 8). Such an exponential-like relationship has been interpreted to support the hypothesis that HMS is produced through aerosol water (Song et al., 2019a;Ma et al., 2020). This is because the amount of aerosol water also exhibits an exponential relationship with RH (Song et al., 2018(Song et al., , 2019b. Interestingly, our model simulations using aqueous clouds as the only medium can obtain a similar relationship between HMS and RH, which reduces the credibility of this aerosol water hypothesis and lends support to the hypothesis that cloud water is a dominant medium for HMS formation (Moch et al., 2018(Moch et al., , 2020. Global atmospheric models, including the numerical weather prediction model employed in the MERRA-2 reanalysis, are usually not capable of resolving subgrid cloud processes, and cloud properties are parameterized using an RH-related statistical scheme (Molod, 2012;Molod et al., 2015). Thus, it is not surprising to find a relationship between RH and HMS in the simulations.

Knowledge gained and remaining gaps
The different spatiotemporal patterns of the HMS levels modeled by the control and default simulations indicate the importance of an appropriate representation of heterogeneous cloud sulfur chemistry. Our modeling suggests that photochemical oxidizing capacity is a key influencing factor for HMS formation because it affects the competition of HCHO with oxidants (e.g., O 3 and H 2 O 2 ) for SO 2 . This factor is partly responsible for the distinct seasonality in HMS  modeled by the default simulation. On a regional scale, the most suitable place for the formation and existence of HMS is parts of East Asia in the lower troposphere during the cold season. Aqueous clouds are the major medium for HMS chemistry since the model simulations can reasonably reproduce both the observed HMS levels and the relationship between in situ HMS and RH when assuming this as the only medium. Aerosol water may play a role if the rate constant of HMS formation is greatly enhanced in concentrated solutions. This finding is consistent with several studies (Jacob, 1986;Olson and Hoffmann, 1989;Whiteaker and Prather, 2003;Moch et al., 2018Moch et al., , 2020. Quantitative observations of HMS are sparse, and more data are required to validate the model. The quantification of HMS in different seasons and over different photochemical conditions, with samples analyzed shortly after collection to minimize potential HMS decomposition, is particularly valuable. Two knowledge gaps are identified from our sensitivity simulations. First, the key uncertain factor in the model is k f , the rate constant for HMS formation. Large discrepancies exist among existing laboratory experiments (Fig. 2), and future laboratory studies are required to narrow its uncertainty. Second, the lack of kinetic and thermodynamic data for HMS chemistry in concentrated solutions poses a key challenge to modeling HMS processing in aerosol water, and new laboratory studies are needed. Also, we did not consider the uncertainty in the meteorological reanalysis. Although MERRA-2 assimilates extensive observations and represents the atmo-spheric states accurately, cloud properties are modeled exclusively. Studies have shown biases in seasonal and spatial variations of cloudiness when comparing the reanalysis data with lidar and satellite observations (Kennedy et al., 2011;Stengel et al., 2018;Miao et al., 2019). Previous work has found that simulated HMS is extremely sensitive to cloud distributions and properties, and these biases in MERRA-2 cloud properties have been identified as a possible source of error for simulation of HMS driven by MERRA-2 (Moch et al., 2018. Recently, the quantum chemical calculations by Chen and Zhao (2020) suggested that hydroxymethyl sulfite (HMSi), an isomer of HMS, might also be produced by an aqueous reaction of HCHO and SO 2 . The laboratory experiments of De Haan et al. (2020) demonstrated that HMS was one of the major products from the aqueous processing of glyoxal monobisulfite (CH(OH) 2 CH(OH)SO − 3 ), the adduct of glyoxal and SO 2 . These new mechanisms need to be considered in future model studies.

Conclusions
Based on appropriate implementation of heterogeneous HMS chemistry and assuming aqueous clouds as the only medium, the global GEOS-Chem model can reasonably reproduce the observations of HMS in Beijing and three other sites. The modeled HMS concentrations and HMS / sulfate ratios show a clear seasonal pattern with higher values in the cold period. The spatial distributions of HMS in descending order are East Asia, Europe, and North America. Our model simulations find the highest average HMS concentrations (1-3 µg m −3 ) and HMS / sulfate molar ratios (0.1-0.2) in northern China during the winter season. Photochemical oxidizing capacity affects the competition of HCHO with oxidants (e.g., O 3 and H 2 O 2 ) for SO 2 , and is a key factor influencing HMS formation. Aqueous clouds act as the primary medium for HMS chemistry, agreeing with prior studies (Moch et al., 2018(Moch et al., , 2020, while aerosol liquid water could play a role if the rate constant for HMS formation is greatly enhanced. This study identifies future research needs. Laboratory experiments should reduce the uncertainty in the formation rate constant of HMS and determine the kinetics for HMS chemistry in concentrated solutions. More field observations of HMS, especially its quantification in different seasons and photochemical conditions, are helpful to validate the model. It is important for future observations of HMS to limit possible HMS decomposition so that the observations can be compared directly with model results. Limiting HMS decomposition can be achieved by analyzing samples shortly after collection and by using analysis methods that are designed to limit HMS decomposition Ma et al., 2020;Moch et al., 2020). The coarse resolution of the global model does not allow it to capture day-to-day observations at a single site, and we are preparing another paper to demon-strate the capacities of regional model with a finer resolution to reproduce individual haze events in northern China.
Code and data availability. The standard GEOS-Chem model is available at: https://doi.org/10.5281/zenodo.3860693 (The International GEOS-Chem User Community, 2020). The code changes made in this study are available at: https://github.com/shaojiesong/ GC1210_sulfchem_Song2020 (Song, 2021). The laboratory and observational data used in this study were all obtained from published papers and books.
Author contributions. ShS initiated the study, carried out analysis, and wrote the initial draft. All authors helped interpret the data, provided feedback, and commented on the manuscript.