Revisiting the critical role of stabilized Criegee intermediates (sCIs) in sulfuric acid formation: coupling mechanistic updates with interpretable machine learning

Sulfuric acid, formed via gas-phase SO₂ oxidation, drives new particle formation and secondary sulfate aerosol under low relative humidity (RH). Beyond the well-established OH-mediated pathway, stabilized Criegee intermediates (sCIs) from alkene ozonolysis provide an additional SO₂ oxidation route whose rate ($L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$) and fractional contribution (μ_sCIs) vary markedly across environments. Although prior studies have documented this spatiotemporal heterogeneity, quantifying the individual factor contributions governing these differences remains lacking. Furthermore, because sCI and OH formation share common precursors yet differ mechanistically, μ_sCIs can alter the sensitivity of total SO₂ oxidation ($L_{{\mathrm{SO}}_{\mathrm{2}}}$) to precursor variations, but the magnitude of this influence remains unclear. Here, we updated Criegee-related reactions in MCM v3.3.1 following recent IUPAC recommendations, conducted extensive AtChem box-model simulations, and trained interpretable XGBoost surrogate models to attribute $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ and μ_sCIs to their controlling factors and examine how μ_sCIs regimes modulate $L_{{\mathrm{SO}}_{\mathrm{2}}}$ sensitivity. $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ was primarily controlled by RH and isoprene. The dominant drivers of μ_sCIs differed between day and night: daytime factors included alkene speciation, O₃, and the alkene fraction in VOCs; nighttime factors included alkene speciation, RH, and NO_x. Under high-μ_sCIs regimes, positive $L_{{\mathrm{SO}}_{\mathrm{2}}}$responses to O₃, VOCs, and alkene fraction were consistently amplified relative to low-μ_sCIs regimes. Field observation-based modeling of a summer episode in Wuhai supported these regime-dependent patterns. These findings indicate that reducing VOC emissions provides a more effective mechanism for regulating sulfuric acid formation under conditions with active alkene ozonolysis.

1 Introduction

The environmental and climatic impacts of fine particulate matter (PM_2.5) are strongly influenced by sulfate aerosol (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), which plays a key role in regulating aerosol acidity (Yao et al., 2018), hygroscopicity (Gunthe et al., 2011), and direct radiative forcing (Quaas et al., 2022). Although a minor fraction of particulate sulfate is directly emitted from anthropogenic sources such as coal and biomass combustion (Dai et al., 2019), industrial activities (Song et al., 2025), and transportation (Cao et al., 2025), it is formed predominantly from the atmospheric oxidation of sulfur dioxide (SO₂) via gas-phase, aqueous-phase, and heterogeneous pathways (Ye et al., 2023). During severe haze events, elevated relative humidity and aerosol loadings significantly enhance aqueous-phase and heterogeneous SO₂ oxidation. These multiphase processes dominate sulfate formation and have therefore been the focus of extensive mechanistic investigations (Cao et al., 2021). Nevertheless, the gas-phase oxidation of SO₂ remains critically important because its product, sulfuric acid (H₂SO₄), is a key precursor for atmospheric new particle formation (NPF). Atmospheric nucleation and the subsequent growth of newly formed particles represent a major source of aerosol particles in terms of number concentration, and NPF has been frequently observed even in highly polluted urban environments characterized by strong condensation sinks (Yao et al., 2018). Moreover, over the past decade, stringent emission-control policies in East Asia, particularly in China, have substantially reduced PM_2.5 concentrations (Geng et al., 2021), thereby weakening conditions favorable for aqueous-phase and heterogeneous sulfate formation. In this progressively cleaner atmosphere, the relative importance of gas-phase pathways is expected to increase (Wang et al., 2025). Therefore, under increasingly stringent air pollution control targets (WHO, 2021), the role of SO₂ gas-phase oxidation pathways in secondary aerosol production is expected to become more important and warrants further investigation.

Hydroxyl radicals (OH) have long been recognized as the dominant oxidant in gas-phase SO₂ oxidation. Since the pioneering work of Cox and Penkett (1971) first demonstrated that Criegee intermediates (CIs, zwitterionic carbonyl oxides) react with SO₂ in the gas phase, extensive laboratory and theoretical studies (Taatjes et al., 2013; Vereecken et al., 2012) have advanced understanding of the formation, stabilization, and fate of CIs. Stabilized Criegee intermediates (sCIs) are therefore increasingly considered an important complementary pathway for gas-phase SO₂ oxidation and a potentially significant contributor to H₂SO₄ formation. Multiple field observational studies conducted in the boreal forests in Finland (Mauldin et al., 2012), the SMEAR II station in Finland and the Hohenpeissenberg station in Germany (Boy et al., 2013), in urban Beijing (Guo et al., 2021), and at a remote Mediterranean site strongly influenced by biogenic VOC emissions (Kukui et al., 2021) suggest that sCIs can provide an additional, environment-dependent pathway for SO₂ oxidation. Notably, as the direct measurement of atmospheric sCIs concentrations remains a formidable challenge, observational studies have primarily relied on the analysis of precursors and OH radical concentrations to provide qualitative evidence. Consequently, numerical simulations grounded in CIs chemical kinetics are indispensable for quantifying the contribution of sCIs to the gas-phase SO₂ oxidation rate.

To date, numerical modeling studies have extensively evaluated the contribution of sCIs to SO₂ oxidation across a broad range of spatial scales and environmental conditions. For instance, utilizing a steady-state approximation, Khan et al. (2017) estimated that the OH-driven SO₂ oxidation rate in urban London (0.64 Tg yr⁻¹) was approximately 17 times that of the sCI pathway (37.6 Gg yr⁻¹). However, proxy calculations of H₂SO₄ concentrations have revealed pronounced regional disparities: at the semi-pristine boreal forest site in Hyytiälä, the photochemical pathway (SO₂ + OH) and the alkene ozonolysis pathway (SO₂ + sCIs) accounted for approximately 34 % and 60 % of H₂SO₄ production, respectively, whereas in heavily polluted Beijing, the corresponding proportions shifted to 28 % and 22 % (Dada et al., 2020). Vereecken et al. (2017) predicted that in equatorial regions, the contribution of sCIs to the gas-phase conversion of SO₂ to sulfuric acid could reach up to 75 %. Furthermore, simulations using the STOCHEM-CRI global atmospheric chemistry and transport model demonstrated (Khan et al., 2018) that the SO₂ oxidation contribution from sCIs significantly exceeds that of OH radicals in terrestrial rainforests and high-latitude boreal forests. This pronounced variability is primarily driven by two intertwined factors. First, considerable uncertainty in numerical simulations arises from inherent limitations in our understanding of CIs chemical kinetics. With continuous breakthroughs in monitoring technologies, the scientific understanding of CI evolution processes is dynamically evolving, resulting in discrepancies in the reaction rate coefficients incorporated into models across different periods. Notably, the unimolecular decomposition rates of CIs/sCIs and the bimolecular reaction rates of sCIs with water monomers and dimers (H₂O/(H₂O)₂) exert a decisive influence on numerical calculations (Sarwar et al., 2014). For example, a modeling study under summer conditions in the eastern United States demonstrated that once the competitive reaction between sCIs and H₂O/(H₂O)₂ was introduced, the sCI-driven enhancement in surface sulfate plummeted from 18 % to less than 0.5 % (Li et al., 2013). Similarly, in simulations of nocturnal power plant plumes in the southeastern US, whether CI thermal decomposition was considered led to a seven-fold difference in the estimated contribution to secondary sulfate aerosol (SSA) formation (Meidan et al., 2019). Second, pronounced regional heterogeneity in precursor composition, such as biogenic isoprene and monoterpenes in forested regions and complex anthropogenic alkene mixtures in urban areas, together with differences in relative humidity and NO_x levels, leads to highly variable sCI production and scavenging fluxes. While existing studies have made progress in quantifying the contribution of sCIs to gas-phase SO₂ oxidation in specific environments and identifying their precursors, the relationship between sCIs and their precursors remains unclear. Consequently, a systematic understanding of the critical driving factors dictating the relative contribution of sCIs is still lacking. Moreover, the production and loss processes of OH and sCIs are not mutually independent; rather, they are deeply coupled within the underlying chemical mechanisms (Lu et al., 2019). This coupled nature dictates that when evaluating the drivers controlling the contribution of the sCIs + SO₂ pathway, their concurrent impacts on the OH + SO₂ pathway and the overall SO₂ oxidation rate cannot be overlooked. Fundamentally, the interplay between precursor abundances and competitive sinks determines the fractional contribution of the sCIs pathway. Yet, how this interplay shapes the sensitivity of the total gas-phase SO₂ oxidation to environmental drivers has received remarkably little attention. In particular, it remains elusive how the dominant controlling factors and their directional impacts shift as the atmospheric gas-phase SO₂oxidation regime transitions from being strictly OH-dominated to being co-driven by both OH and sCIs.

To this end, this study presents a systematic assessment of the role of sCIs in atmospheric SO₂ gas-phase oxidation. As a prerequisite for all subsequent analyses, we first updated the gas-phase kinetics of Criegee intermediates in the Master Chemical Mechanism (MCM v3.3.1, via website: https://mcm.york.ac.uk, last access: 28 August 2026) using the latest evaluated rate coefficients recommended by the International Union of Pure and Applied Chemistry (IUPAC), thereby minimizing the propagation of mechanistic uncertainties into our conclusions. Utilizing an atmospheric box model coupled with this revised mechanism, we integrated interpretable machine learning techniques to quantify the importance of key controlling factors on sCI-mediated SO₂oxidation. Building on this insight, we employed a robust experimental design to construct a machine learning surrogate model capable of predicting sCI contributions across diverse and highly variable environmental conditions. Furthermore, analytical methods, including SHAP, Sobol sensitivity analysis, Partial Dependence Plots (PDPs), and individual conditional expectation (ICE) plots, were utilized to quantify how these factors dictate the relative contribution of sCIs, and to elucidate the differential sensitivities of total SO₂ oxidation under high- and low-sCI contribution regimes. Finally, leveraging online hourly observational data of atmospheric pollutants, we assessed the disparities in their sensitivity relationships with SO${}_{\mathrm{4}}^{\mathrm{2}-}$ across varying sCI contribution levels, ultimately validating the reliability of our box model-derived surrogate model.

2 Data and methodology

2.1 Updates to alkene ozonolysis and Criegee intermediate chemistry in MCM v3.3.1

The Master Chemical Mechanism (MCM v3.3.1; available at http://mcm.york.ac.uk, last access: 28 August 2026) is a near-explicit chemical mechanism that describes the gas-phase oxidation of 143 volatile organic compounds (VOCs) to carbon dioxide (CO₂) and water (H₂O). In MCM v3.3.1, the rate coefficients and associated reaction mechanisms for alkene ozonolysis reactions follow the recommendations of Jenkin et al. (1997, 2015) and Saunders et al. (2003). Although this mechanism provides a detailed representation of CI formation, loss, and subsequent chemical evolution, several key rate coefficients and parameters differ substantially from the latest IUPAC recommendations (Cox et al., 2020). In addition, several important processes are not explicitly represented, including the specific fates of CIs, the unimolecular decomposition of sCIs, the reactions of sCIs with water dimers ((H₂O)₂), and the distinct reactivities of sCI stereoisomers. To address these limitations, we updated the gas-phase alkene + O₃ chemistry in MCM v3.3.1 primarily in accordance with current IUPAC recommendations. The updates include the alkene ozonolysis rate coefficients (k_O3+alkene), OH yields (Y_i,OH), sCI yields (Y_i,sCI), and sCI bimolecular reaction rate coefficients (k_b), as well as accounting for sCI unimolecular decomposition pathways and adding sCI + (H₂O)₂ reactions. Where IUPAC provides temperature-dependent kinetic expressions, these were adopted; in addition, a temperature-dependent rate expression was incorporated for the CH₂OO + SO₂ reaction, based on Onel et al. (2021). Specifically, the branching ratios for the unimolecular decomposition or isomerization pathways of CIs and sCIs were constrained using reported OH yields and the yields of other relevant products. For sCIs exhibiting stereoisomerism, including CH₃CHOO, C₂H₅CHOO, and the C₄ intermediates formed during isoprene ozonolysis, the Z-stereoisomers were assumed to undergo predominantly rapid unimolecular decomposition under tropospheric conditions, such that their contribution to bimolecular reactions was considered negligible (Newland et al., 2015). Accordingly, in the updated mechanism, hereafter denoted MCM v3.3.1g, stereoisomeric sCIs participating in bimolecular reactions are represented explicitly by their E-stereoisomers. The unimolecular decomposition pathways of sCIs were considered together with the decomposition or isomerization channels of their corresponding chemically activated CIs. The kinetic parameters and product yields for the alkene ozonolysis reactions before and after the updates are summarized in Tables S1 and S2 in the Supplement. The complete set of updated reaction equations, together with the corresponding rate coefficients and explicit branching ratios, is provided in Table S3. These updates were applied to C₁–C₄ Criegee intermediates, which are widely recognized to play a critical role in tropospheric chemistry.

The concentration of (H₂O)₂ is determined by assuming rapid thermal equilibrium between water monomers and dimers. Given this fast atmospheric exchange, the dimer concentration is considered to be in thermodynamic equilibrium with the monomer, as expressed in Eq. (1):

$$\begin{array}{}\text{(1)}& K_{\mathrm{eq}}^{\mathrm{D}}={\displaystyle \frac{\left[{\left({\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right)}_{\mathrm{2}}\right]}{{\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]}^{\mathrm{2}}}},\end{array}$$

where $K_{\mathrm{eq}}^{\mathrm{D}}$ is the temperature-dependent equilibrium constant for dimer formation (Scribano et al., 2006). Following the methodology of Lade et al. (2024), the reaction of sCIs with water dimers is parameterized using an effective third-order rate coefficient, k_eff (cm⁶ molec.⁻² s⁻¹), and the rate of sCI loss via reaction with water dimers, $r_{\mathrm{sCI}+({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{2}}}$, is given by Eq. (2):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{r}_{\mathrm{sCI}+({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{2}}}& =k_{\mathrm{b},\mathrm{dimer}}\left[\mathrm{sCI}\right]\left[{\left({\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right)}_{\mathrm{2}}\right]\\ & =k_{\mathrm{b},\mathrm{dimer}}{K}_{\mathrm{eq}}^{\mathrm{D}}\left[\mathrm{sCI}\right]{\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]}^{\mathrm{2}}\\ & =k_{\mathrm{eff}}\left[\mathrm{sCI}\right]{\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]}^{\mathrm{2}},\end{array}\end{array}$$

where k_eff is defined as the product of the bimolecular rate coefficient for the sCI + (H₂O)₂ reaction (k_b,dimer, cm³ molec.⁻¹ s⁻¹) and $K_{\mathrm{eq}}^{\mathrm{D}}$. To calculate $K_{\mathrm{eq}}^{\mathrm{D}}$ accurately, thermochemical data were retrieved from the Active Thermochemical Tables (ATcT). Standard Gibbs free energies of formation (${\mathrm{\Delta }}_f{G}_T^{\mathit{°}}$) for both the water monomer and dimer were extracted from Tables 1 and 3 in Ruscic (2013). These discrete data points were then used to calculate the reaction Gibbs free energy (${\mathrm{\Delta }}_r{G}_T^{\mathit{°}}$) over the temperature range of 200–360 K, which spans typical tropospheric conditions. The temperature dependence of $K_{\mathrm{eq}}^{\mathrm{D}}$ was subsequently obtained via a linear regression of $\ln \left(K_{\mathrm{eq}}^{\mathrm{D}}\right)$ as a function of $\frac{\mathrm{1}}{T}$, yielding the following empirical van't Hoff–type parameterization in Eq. (3):

$$\begin{array}{}\text{(3)}& K_{\mathrm{eq}}^{\mathrm{D}}\left(T\right)=A\exp \left({\displaystyle \frac{B}{T}}\right),\end{array}$$

where the coefficients A and B represent the intercept and slope, respectively, derived from the fitting procedure, with $A=\mathrm{1.15}\times {\mathrm{10}}^{-\mathrm{23}}$ cm³ molec.⁻¹ and B=1549.32 K.

2.2 Box model configuration and observation data

2.2.1 AtChem model setup and chemical diagnostics

AtChem (Sommariva et al., 2020), an open-source zero-dimensional atmospheric box model, was used in this study to simulate the gas-phase chemical evolution of SO₂, CIs, and their precursors. The AtChem model is coupled with MCM v3.3.1 and MCM v3.3.1g (see Sect. 2.1). H₂O concentrations for sCI loss reactions were calculated from relative humidity, temperature, and atmospheric pressure.

Two types of AtChem simulations were performed in this study. The first used an unconstrained configuration, in which the initial concentrations of selected chemical species were prescribed and allowed to evolve according to the chemical mechanism without time-varying constraints. This configuration was used in Sect. 3.1 to quantify the oxidation capacity of sCIs derived from the ozonolysis of specific alkenes, including ethene, propene, but-1-ene, trans-but-2-ene, cis-but-2-ene, 2-methylpropene, and isoprene. The second used a constrained configuration, in which selected trace gases, meteorological parameters, and photolysis frequencies were prescribed as hourly time-varying inputs from observational data or designed perturbation scenarios. This configuration was used in Sect. 3.2 and 3.3 to evaluate the role of sCIs in gas-phase SO₂ oxidation under atmospherically relevant conditions. A spin-up period of 2–3 d was implemented to initialize reactive free radicals and intermediates, such as OH, HO₂, RO₂, and sCIs, allowing them to reach realistic steady-state concentrations.

Chemical production and loss rates were quantified using the Rate of Production/Destruction Analysis (ROPA/RODA) within AtChem. The gas-phase SO₂ oxidation rate by sCIs was defined according to Eq. (4):

$$\begin{array}{}\text{(4)}& L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}={\sum }_i{k}_{{\mathrm{SO}}_{\mathrm{2}}+\mathrm{sCI}}\left[{\mathrm{sCI}}_i\right]\left[{\mathrm{SO}}_{\mathrm{2}}\right],\end{array}$$

where i denotes an individual sCI species. The OH-driven SO₂ oxidation rate was defined as in Eq. (5):

$$\begin{array}{}\text{(5)}& L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{OH}}=k_{{\mathrm{SO}}_{\mathrm{2}}+\mathrm{OH}}\left[\mathrm{OH}\right]\left[{\mathrm{SO}}_{\mathrm{2}}\right],\end{array}$$

The total gas-phase SO₂ oxidation rate and the fractional sCI contribution were then calculated as in Eqs. (6) and (7):

$$\begin{array}{}\text{(6)}& {\displaystyle }{L}_{{\mathrm{SO}}_{\mathrm{2}}}=L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}+L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{OH}},\text{(7)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{sCIs}}={\displaystyle \frac{L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}}{L_{{\mathrm{SO}}_{\mathrm{2}}}}},\end{array}$$

here, $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ and μ_sCIs provide complementary perspectives on the absolute rate and fractional contribution, respectively, of the sCI pathway.

2.2.2 Observation data

The observational data used in this study were obtained from the automatic monitoring system at the Wuhai City Atmospheric Environment Super Monitoring Station from 9 October 2019 to 30 June 2022. Wuhai is a semi-arid coal-chemical industrial city in Inner Mongolia, China, characterized by relatively high emissions of SO₂ and anthropogenic VOCs (including alkenes), frequent O₃ pollution, and comparatively dry atmospheric conditions. The dataset includes air pollutants (PM_2.5, carbon monoxide (CO), SO₂, nitric oxide (NO), nitrogen dioxide (NO₂), and ozone (O₃)), meteorological variables (wind speed (WS), temperature (T), pressure (P), and relative humidity (RH)), water-soluble inorganic ions (including NO${}_{\mathrm{3}}^{-}$ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$), elemental components (including Fe), VOCs, and photolysis frequencies. PM_2.5, NO, NO₂, SO₂, CO, and O₃ were measured using a Metone BAM-1020 PM_2.5 Monitor, an API T201 NH₃–NO_x Analyzer, an API T100 SO₂ Analyzer, an API T300 CO Analyzer, and an API T400 O₃ Analyzer, respectively. VOCs were measured using an ENTECH BCT-7800 VOCs Analyzer. Water-soluble inorganic ions were measured using a Metrohm MARGA ADI 2080 online ion chromatograph, and elemental components were measured using a CES Xact-625 atmospheric heavy metal analyzer. Meteorological variables (T, RH, P, WS) were measured using a Lufft WS601-UMB six-parameter meteorological sensor. Photolysis frequencies were measured using a Metcon UF-CCD photolysis spectroradiometer. The observational data had an hourly time resolution and were screened for invalid values. Short missing intervals were linearly interpolated when appropriate, whereas longer gaps were excluded from the analysis.

2.3 Machine learning model and interpretation framework

An interpretable machine learning (ML) framework was developed in this study to systematically evaluate the role of sCIs in SO₂ oxidation, encompassing model development, validation, and interpretation. The learning algorithm used throughout was Extreme Gradient Boosting (XGBoost) (Chen and Guestrin, 2016), which constructs an additive ensemble of regression trees by iteratively minimizing a regularized loss function, with each successive tree fitted to the negative gradients (i.e., pseudo-residuals) of the current ensemble's predictions. This formulation enables nonlinear responses and predictor interactions to be represented far more flexibly than a single regression tree. Regularization, subsampling, cross-validation, and early stopping were applied to suppress overfitting, and hyperparameters were tuned via random search or Bayesian optimization. The same sampling–training–testing workflow was applied to every model in this study. Two classes of XGBoost models were constructed (using the XGBoost Python package, version 2.1.3; https://github.com/dmlc/xgboost, last access: 28 August 2026): surrogate models trained on AtChem-simulated chemical evolution, and an observation-based model trained on ambient measurements. Four complementary interpretation methods were used to explain the XGBoost model outputs: SHAP (SHapley Additive exPlanations), Sobol sensitivity analysis, partial dependence plots (PDPs), and individual conditional expectation (ICE) plots.

2.3.1 Model development details

Three surrogate models were developed based on SO₂ gas-phase oxidation simulations from AtChem to successively quantify the key controls on sCI-mediated SO₂ oxidation, identify the atmospheric conditions under which the sCI pathway becomes important, and examine how this importance affects the sensitivity of total SO₂ oxidation. An observation-based model was developed to assess whether the sensitivities inferred from the surrogate models are also reflected in ambient measurements. These XGBoost models are introduced in turn below.

First, a surrogate model was developed to emulate the sCI-mediated SO₂ oxidation rate ($L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$), with SHAP analysis applied to diagnose the contributions of its key controlling factors, namely the principal sCI precursors and competing sinks (most notably NO₂, H₂O, and (H₂O)₂). The feature set therefore comprises O₃, RH, NO₂, and six alkene precursors: ethene, propene/but-1-ene, trans-but-2-ene, cis-but-2-ene, 2-methylpropene, and isoprene, where propene and but-1-ene are merged into a single feature owing to the close similarity of their ozonolysis mechanisms. To capture individual perturbations and interactions across all features, 3898 simulation scenarios were designed spanning the prescribed feature ranges. In each scenario, the feature values prescribed the initial input concentrations for a 1 h unconstrained AtChem simulation, with SO₂ held constant; the mean $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ over the simulated hour, obtained under both MCM v3.3.1 and v3.3.1g, served as the target output. The resulting feature matrix and corresponding simulation outputs together constitute the training dataset for the surrogate model. The prescribed feature ranges, full scenario design, training details, and model performance are provided in Sect. S1 in the Supplement. Note that in all simulation scenarios, temperature and pressure were held constant at T=290 K and P=888 hPa. Temperature-dependent results for the full scenario set are provided in Sect. S4.

Second, to identify the atmospheric conditions that favor sCI contribution to total gas-phase SO₂ oxidation, a surrogate model was developed to predict the fractional importance of the sCI pathway (μ_sCIs), combined with Sobol sensitivity analysis, PDPs, and ICE plots. Unlike $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$, which depends solely on sCI production and loss, μ_sCIs is also governed by the parallel OH-mediated oxidation pathway and therefore requires a feature set that captures both sCI and OH chemistry. Accordingly, the feature set comprises O₃, NO_x, NO₂ fraction (NO₂ % = NO₂/NO_x), total VOC concentration (VOCs), alkene fraction (alkene % = ∑alkenes/VOC), RH, and the Aggregated Reactivity Index (ARI). The NO_x–NO₂ % and VOC–alkene % pairs jointly characterize the oxidant state and reactive VOC composition while minimizing multicollinearity among features. ARI is an ordinal variable encoding the sCI-oxidation potential of the prevailing alkene mixture, derived directly from the importance ranking of alkene features established by the $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ surrogate model. The training samples generated for this surrogate differ from those of the $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ model. Rather than prescribing initial concentrations, the feature values here represent normalized scaling factors applied to time-varying constraints, as the AtChem model was run in a constrained configuration. The baseline constraints are the three-year mean diurnal cycles of hourly observations at the Wuhai supersite (Sect. 2.2.2), and the perturbation range of each feature spans the 5th–95th percentile range of the observed concentrations (normalized range is [0, 1]); ARI, as an ordinal discrete variable, is the sole exception and takes integer values of 0, 1, or 2, corresponding to the three tiers of sCI-oxidation potential. Because photolysis frequencies play a decisive role in OH production and radical cycling, daytime and nighttime conditions were treated as entirely separate modeling problems, with target time periods of 08:00–18:00 and 22:00–05:00 LT, respectively. For each scenario, the constraints of the target period were adjusted according to the prescribed feature values, while all remaining constraints (non-feature variables) were held at baseline (Fig. S2 in the Supplement); a constrained AtChem simulation under MCM v3.3.1g then yielded the mean μ_sCIs and $L_{{\mathrm{SO}}_{\mathrm{2}}}$ over the target period as the target outputs. The resulting feature matrix and simulated μ_sCIs values together constitute the training datasets for the daytime and nighttime μ_sCIs surrogate models. These daytime and nighttime datasets were then each divided into low- and high-μ_sCIs subsets based on the median μ_sCIs value, and a further pair of surrogate models was developed to predict $L_{{\mathrm{SO}}_{\mathrm{2}}}$ separately in each regime, combined with partial dependence analysis to assess whether the sensitivity of $L_{{\mathrm{SO}}_{\mathrm{2}}}$ to each feature depends on the relative importance of the sCI pathway. The baseline constraints, the mapping from normalized feature values to constraint magnitudes, full scenario design, training details, and model performance are provided in Sect. S2.

Finally, to assess whether the regime-dependent sensitivities inferred from the surrogate models are also expressed in the real atmosphere, we developed an observation-based model trained on ambient observations from the Wuhai supersite (9 October 2019 to 30 June 2022, hourly resolution) to predict SO${}_{\mathrm{4}}^{\mathrm{2}-}$. The features, including PM_2.5, NO${}_{\mathrm{3}}^{-}$, SO₂, O₃, total VOCs, alkene %, NO_x, NO₂ %, RH, wind speed, and Fe, were selected to represent sulfate aerosol precursor abundance, photochemical oxidation capacity, aerosol loading, secondary inorganic aerosol formation, meteorological conditions, and gas-phase and aqueous-phase oxidation pathways. SHAP analysis was applied to evaluate the differential contributions of key features to SO${}_{\mathrm{4}}^{\mathrm{2}-}$ under contrasting μ_sCIs regimes, enabling comparison with the regime-dependent sensitivities diagnosed from the AtChem-based surrogate models. Training details and model performance are provided in Sect. S3.

2.3.2 Model interpretation methods

In this study, XGBoost was not used merely as a black-box surrogate for predicting chemical reaction rates or target pollutant concentrations. Instead, it served as a computationally efficient emulator of the box-model simulations, enabling systematic sensitivity analysis and interpretation of the chemical response within the prescribed parameter space. To interpret the trained XGBoost models, we combined SHAP analysis, Sobol sensitivity analysis, partial dependence plots, and individual conditional expectation plots. These four approaches were selected to provide complementary perspectives spanning global versus local attribution, main effects versus interaction effects, and average versus sample-level response behavior. Specifically, SHAP values were used to quantify the relative importance and direction of predictor effects at both individual-sample and aggregate levels, while SHAP interaction values were used to diagnose key two-factor interactions. Sobol sensitivity analysis was used to provide a variance-based decomposition of the model output, separating the independent contribution of each input factor from interaction-driven contributions. PDPs were used to visualize the average nonlinear response of the model output to selected predictors and to identify potential thresholds, whereas ICE plots were used to examine whether these average responses were consistent across individual samples or instead varied among different chemical regimes, revealing heterogeneity that PDPs alone cannot capture. For models trained on controlled numerical experiments, the feature values were prescribed through independent, design-of-experiments-based sampling, and the target outputs were generated by a deterministic, mechanism-based chemical model (AtChem). Because the training data are thus not subject to uncontrolled confounding, and because XGBoost was trained to approximate the fixed relationships embedded in the chemical mechanism, the resulting SHAP, Sobol, PDP, and ICE analyses were interpreted as model-estimated chemical sensitivities that reasonably reflect the causal structure inherent in the chemical mechanism, within the sampled parameter space. For the model trained on ambient observations, by contrast, the input features are inherently correlated, and the observed feature–outcome relationships may be subject to unknown confounding. Consequently, SHAP-based interpretations for this model were treated as statistical associations used to evaluate consistency with the chemical sensitivities inferred from the AtChem-based surrogate models, rather than as evidence of direct causal effects.

1. SHAP (SHapley Additive exPlanations) is a game-theoretic approach for interpreting machine learning model predictions (Lundberg and Lee, 2017; Lundberg et al., 2018). Each input feature is treated as a “player” in a cooperative game, and the model prediction for any individual sample is expressed as an additive decomposition of per-feature contributions. Specifically, for a sample x_i with N features, the predicted value f (x_i) is decomposed as Eq. (8):

   $$\begin{array}{}\text{(8)}& f\left(x_i\right)={\mathit{\varphi }}_{\mathrm{0}}+{\sum }_{j=\mathrm{1}}^N{\mathit{\varphi }}_j\left(f,x_i\right),\end{array}$$

   where ϕ₀=E[f(x)] is the base value, i.e., the expected model output averaged over the entire dataset, and ϕ_j(f,x_i) is the SHAP value of feature j for sample x_i, quantifying the contribution of feature j to the deviation of the prediction from the base value ϕ₀. The SHAP value ϕ_j(f,x_i) is computed as the weighted average marginal contribution of feature j across all possible subsets of the remaining features, as shown in Eq. (9):

   $$\begin{array}{}\text{(9)}& {\mathit{\varphi }}_j\left(f,x_i\right)={\sum }_{S\subseteq F\backslash \left\{j\right\}}{\displaystyle \frac{\left|S\right|\mathrm{!}\left(\left|F\right|-\left|S\right|-\mathrm{1}\right)\mathrm{!}}{\left|F\right|\mathrm{!}}}\left[f\left(S\cup \left\{j\right\}\right)-f\left(S\right)\right],\end{array}$$

   where F is the full set of features, S denotes a subset of features that excludes feature j, and f (S) denotes the conditional expectation of the model output given the features in S. The sign of ϕ_j(f,x_i) indicates the direction of the feature effect for a specific sample (positive: the feature j increases the predicted value relative to the baseline; negative: it decreases it), and the magnitude indicates the strength of that effect. For XGBoost models, SHAP values are computed exactly and efficiently using the TreeExplainer algorithm (Lundberg et al., 2020).

   The global importance of feature j is summarized by its mean absolute SHAP value (Eq. 10):

   $$\begin{array}{}\text{(10)}& \stackrel{\mathrm{‾}}{\left|{\mathit{\varphi }}_j\right|}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\left|{\mathit{\varphi }}_j\left(f,x_i\right)\right|,\end{array}$$

   where n is the total number of samples. $\stackrel{\mathrm{‾}}{\left|{\mathit{\varphi }}_j\right|}$ represents the average magnitude of feature j's contribution to model predictions across the dataset and is used to rank features by their global importance.

   Pairwise feature interactions are further quantified via SHAP interaction values (Lundberg et al., 2019). The mean absolute interaction value between features j and p is defined as Eq. (11):

   $$\begin{array}{}\text{(11)}& \stackrel{\mathrm{‾}}{\left|{\mathit{\varphi }}_{jp}\right|}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\left|{\mathit{\varphi }}_{jp}\left(f,x_i\right)\right|,\end{array}$$

   where ϕ_jp(f,x_i) is the SHAP interaction value between features j and p for sample x_i, quantifying the portion of their joint effect on the prediction.

2. A variance-based global sensitivity analysis was performed using the Sobol method (Saltelli et al., 2010; Sobol, 2001). Sobol analysis decomposes the variance of the model output into contributions from individual input factors and their interactions over a prescribed input space and probability distribution. It thereby provides a global measure of sensitivity that is not determined by the empirical distribution of the training data.

   For a model output Y=f(x), where x₁, x₂, …, x_N denote N mutually independent input factors, the total output variance, Var(Y), can be decomposed into contributions from individual factors and their interactions (Eq. 12):

   $$\begin{array}{}\text{(12)}& \mathrm{Var}\left(Y\right)={\sum }_{j=\mathrm{1}}^N{V}_j+{\sum }_{j<k}{V}_{jk}+\mathrm{\dots },\end{array}$$

   where $V_j={\mathrm{Var}}_{x_j}\left[E_{x_{\sim j}}\left(Y|x_j\right)\right]$ is the variance attributable to the main effect of x_j, whereas V_jk represents the additional variance attributable to the interaction between x_j and x_k. Higher-order terms represent interactions involving three or more input factors.

   Two Sobol sensitivity indices were computed for each input factor. The first-order index, S_j, quantifies the proportion of output variance attributable to the main effect of x_j (Eq. 13):

   $$\begin{array}{}\text{(13)}& S_j={\displaystyle \frac{V_j}{\mathrm{Var}\left(Y\right)}},\end{array}$$

   The total-order index, $S_{T_j}$, quantifies the contribution of x_j to the output variance, including its main effect and all interaction effects involving x_j (Eq. 14):

   $$\begin{array}{}\text{(14)}& S_{T_j}={\displaystyle \frac{E_{x_{\sim j}}\left[{\mathrm{Var}}_{x_j}\left(Y|x_{\sim j}\right)\right]}{\mathrm{Var}\left(Y\right)}}=\mathrm{1}-{\displaystyle \frac{{\mathrm{Var}}_{x_{\sim j}}\left[E_{x_j}\left(Y|x_{\sim j}\right)\right]}{\mathrm{Var}\left(Y\right)}},\end{array}$$

   where x_∼j denotes all input factors except x_j. Under the assumed independent input distributions, both indices range from 0 to 1. The difference, $S_{T_j}-S_j$, quantifies the cumulative contribution of interaction effects involving x_j. Thus, values of $S_{T_j}$ substantially larger than S_j indicate that x_j participates strongly in non-additive interactions with other predictors, an effect that cannot be identified from main-effect rankings alone.

   Input samples for the Sobol sensitivity analysis were generated using the Saltelli sampling scheme based on Sobol' quasi-random sequences (Saltelli, 2002). For each continuous input factor, the sampling bounds were defined by the minimum and maximum values in the training dataset, and samples were generated independently over these bounds. A base sample size of 32 768 was used. Following the Saltelli sampling scheme with second-order effects enabled, this resulted in 524 288 surrogate-model evaluations for the seven input factors. This large evaluation budget was made computationally efficient by the trained XGBoost surrogate model, which evaluates each scenario several orders of magnitude faster than the corresponding box-model simulation. For ARI, which was encoded as an ordinal variable with integer values of 0, 1, and 2, the continuous Saltelli samples generated over [0, 2] were rounded to the nearest integer and clipped to the valid range to match the representation used for model training. Sobol indices, including first-order and total-order indices, were estimated using the SALib library (Herman and Usher, 2017). Bootstrap 95 % confidence intervals were estimated for all Sobol indices.

3. Partial dependence plots (PDPs) were used to visualize the model-predicted average relationship between selected factors and the target variable (Friedman, 2001). For a factor x_p and the complementary factor set x_c, the one-dimensional partial dependence function is defined in Eq. (15):

   $$\begin{array}{}\text{(15)}& {\widehat{f}}_p\left(x_p\right)=E_{x_c}\left[f\left(x_p,x_c\right)\right]\approx {\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^{n}f\left(x_p,x_c^{\left(i\right)}\right),\end{array}$$

   where the expectation is approximated by averaging the surrogate model predictions (${\widehat{f}}_p$) over all n training samples, while fixing x_p at each point of a prescribed grid.

   Two-dimensional PDPs were further constructed for selected factor pairs to visualize their joint dependence on the target variable. For a factor pair (x_p, x_q), the two-dimensional partial dependence function is defined in Eq. (16):

   $$\begin{array}{}\text{(16)}& \begin{array}{rl}{\widehat{f}}_{pq}\left(x_p,x_q\right)& =E_{x_c}\left[f\left(x_p,x_q,x_c\right)\right]\\ & \approx {\displaystyle \frac{\mathrm{1}}{n_{\mathrm{MC}}}}{\sum }_{i=\mathrm{1}}^{n_{\mathrm{MC}}}f\left(x_p,x_q,x_c^{\left(i\right)}\right),\end{array}\end{array}$$

   where n_MC denotes a Monte Carlo subsample of the training data used to reduce computational cost over the two-dimensional grid.

4. Individual conditional expectation (ICE) plots were used to complement the PDP analysis by examining the model response at the individual-sample level (Goldstein et al., 2015). Whereas a PDP describes the average model-predicted relationship between a factor and the target variable across samples, an ICE plot displays a separate conditional response curve for each sample. For each curve, the factor of interest is varied across a prescribed grid, while all remaining factors are held at their sample-specific values.

Figure 1 SHAP beeswarm summary plots for the XGBoost surrogate models of $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ under (a) MCM v3.3.1 and (b) MCM v3.3.1g. Each point represents one simulation scenario, positioned according to its SHAP value (molec. cm⁻³ s⁻¹) for the corresponding feature. Point color denotes the feature value from low to high. Features are ordered by mean absolute SHAP value, so that those at the top exert the strongest overall influence on the predicted sCI-driven SO₂ oxidation rate. (c) Comparison of SHAP-value distributions for individual features between the XGBoost surrogate models trained on MCM v3.3.1 and MCM v3.3.1g. Boxplots summarize the full distribution of SHAP values across all simulation scenarios. Differences between mechanism versions reflect how the kinetic update changes the magnitude of feature effects. (d) Relative feature contributions to the XGBoost surrogate models of $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ under MCM v3.3.1 and MCM v3.3.1g, quantified as the normalized mean absolute SHAP value of each feature. Contributions are expressed as fractions of the total sum of mean absolute SHAP values across all features for each model. This comparison highlights how the mechanism update redistributes the overall importance among ozonolysis precursors and competing sinks.

3 Results and discussion

3.1 Effects of key controlling factors on sCI-mediated SO₂ oxidation

sCIs are short-lived reactive intermediates generated from the ozonolysis of alkenes. Their contribution to gas-phase SO₂ oxidation is governed jointly by their atmospheric steady-state concentration and their bimolecular rate coefficient with SO₂. The steady-state concentration is itself determined by the initial ozonolysis rate coefficient, the total sCI yield, and the structure-dependent reactivities of individual sCI species toward other competing sinks – most notably H₂O, (H₂O)₂, and NO₂ – which collectively regulate the fraction of sCIs available to react with SO₂. We therefore introduce the sCI-mediated SO₂ oxidation rate, $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$, as the target variable, which simultaneously accounts for sCI production and the effective SO₂ oxidation capacity. Figure 1a and b show that, in both MCM v3.3.1 and MCM v3.3.1g, alkene and O₃ concentrations are positively correlated with $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$, whereas RH and NO_x exhibit marked negative correlations. This pattern is consistent with alkenes and O₃ driving sCI production, while H₂O/(H₂O)₂ and NO₂ act as competing sinks that reduce the fraction of sCIs available to react with SO₂. Notably, updating the mechanism substantially increases the magnitude of the absolute SHAP (|SHAP|) values across all features (Fig. 1c), with the average |SHAP| value ranging from 211.9 to 11 177.3 molec. cm⁻³ s⁻¹ in MCM v3.3.1, versus 694.0 to 30 429.5 molec. cm⁻³ s⁻¹ in MCM v3.3.1g. This indicates that the updated mechanism yields a substantially stronger sensitivity of $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ to perturbations in environmental and precursor variables, such that the same change in a given feature produces a larger response in $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$. Importantly, the update does not amplify the influence of all features to the same degree; rather, it also reshapes the relative importance attributed to individual alkenes (Fig. 1d). In MCM v3.3.1, the alkene relative importance follows the order: trans-but-2-ene (TBUT2ENE, 21.5 %) > cis-but-2-ene (CBUT2ENE, 17.6 %) > isoprene (C₅H₈, 9.5 %) > 2-methylpropene (MEPROPENE, 4.3 %) > ethene (C₂H₄, 0.7 %) > propene/but-1-ene (C₃C₄, 0.5 %). In MCM v3.3.1g, isoprene becomes the most important alkene feature (23.8 %), followed by trans-but-2-ene (11.0 %), cis-but-2-ene (9.7 %), and 2-methylpropene (6.0 %), while ethene (0.9 %) and propene/but-1-ene (0.7 %) remain minor contributors. The approximately sixfold increase in the average isoprene |SHAP| value, rising from 3773.1 to 24 500.7 molec. cm⁻³ s⁻¹ after the mechanism update, is markedly larger than the corresponding shift observed for any other alkene. Indeed, prior studies have shown that regions with intense biogenic emissions tend to exhibit elevated sCI-mediated SO₂ oxidation rates, an effect attributed primarily to sCIs derived from isoprene (Hata et al., 2023; Kukui et al., 2021). It should be noted, however, that the relative importance ranking obtained here reflects each alkene's attributed contribution when all alkenes are perturbed over an identical concentration range (0.1–20 ppb); by contrast, isoprene's prominence in previous field and modeling studies may partly reflect its typically higher ambient concentration relative to other alkenes, rather than solely an intrinsically higher sCI reactivity. Within the present perturbation ranges (O₃: 2–90 ppb; each alkene: 0.1–20 ppb; RH: 10 %–80 %), RH emerges as the most influential feature overall, exceeding the attributed importance of every alkene (Fig. 1d). This is consistent with prior kinetic evidence that H₂O/(H₂O)₂ dominates the bimolecular loss of sCIs, such that most sCIs have limited capacity to survive under high humidity (Cox et al., 2020; Lade et al., 2024).

Figure 2 Heat maps of pairwise SHAP interaction values among the nine input features for the XGBoost surrogate models of $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ under (a) MCM v3.3.1 and (b) MCM v3.3.1g. Colours denote the mean absolute SHAP interaction value for each feature pair across all simulation scenarios. Larger values indicate stronger non-additive interactions in the surrogate prediction, thereby highlighting the precursor–sink combinations that most strongly modulate sCI-driven SO₂ oxidation.

The changes in mean absolute SHAP interaction values before and after the update (Fig. 2) further clarify the mechanistic origin of the differential amplification magnitudes in alkene attribution. In MCM v3.3.1g, the elevated SHAP interaction values involving RH with both O₃ and the individual alkenes indicate that MCM v3.3.1 substantially underestimated the competitive scavenging of sCIs by H₂O/(H₂O)₂, thereby masking the strong regulatory role of humidity in sCI-mediated SO₂ oxidation. In parallel, the interaction effects between O₃ and the individual alkenes are also enhanced in MCM v3.3.1g, indicating that MCM v3.3.1 likewise underestimated the sCI + SO₂ oxidation capacity. Critically, the relative degree of underestimation in these two competing removal pathways differs across alkene systems. For isoprene, the RH–C₅H₈ interaction increased from 928.0 to 5094.0 molec. cm⁻³ s⁻¹ (a 5.5-fold increase), whereas the O₃–C₅H₈ interaction increased from 697.3 to 4435.3 molec. cm⁻³ s⁻¹ (a 6.4-fold increase). This indicates that MCM v3.3.1 underestimated the reactivity of its derived sCIs toward SO₂ to a greater extent than their reactivity toward H₂O/(H₂O)₂; a comparable pattern is observed for 2-methylpropene. For the remaining alkenes, by contrast, the underestimation of the sCI+H₂O/(H₂O)₂ pathway predominates.

Based on the MCM v3.3.1g results, isoprene contributes more than twice as much as cis-/trans-but-2-ene, 3.8 times as much as 2-methylpropene, and roughly thirty times more than ethene and propene/but-1-ene in terms of their attributed influence on $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ (Fig. 1c). At 298.15 K, the ozonolysis rate coefficients span more than two orders of magnitude (Table S1), following the order: trans-but-2-ene ($\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{16}}$ cm³ molec.⁻¹ s⁻¹) > cis-but-2-ene ($\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{16}}$ cm³ molec.⁻¹ s⁻¹) > isoprene ($\mathrm{1.28}\times {\mathrm{10}}^{-\mathrm{17}}$ cm³ molec.⁻¹ s⁻¹) ≈ 2-methylpropene ($\mathrm{1.15}\times {\mathrm{10}}^{-\mathrm{17}}$ cm³ molec.⁻¹ s⁻¹) ≈ propene ($\mathrm{1.05}\times {\mathrm{10}}^{-\mathrm{17}}$ cm³ molec.⁻¹ s⁻¹) ≈ but-1-ene ($\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{17}}$ cm³ molec.⁻¹ s⁻¹) > ethene ($\mathrm{1.56}\times {\mathrm{10}}^{-\mathrm{18}}$ cm³ molec.⁻¹ s⁻¹). Total sCI yields follow the sequence (Tables S1 and S3): isoprene (0.66; 0.55 CH₂OO +0.08 (CH=CH₂)(CH₃)COO +0.03 (C(CH₃)=CH₂)CHOO) > ethene (0.42; CH₂OO) > propene (0.28; 0.20 CH₂OO +0.08 CH₃CHOO) > but-1-ene (0.23; 0.15 CH₂OO +0.08 C₂H₅CHOO) ≈ trans-but-2-ene (0.22; CH₃CHOO) ≈ cis-but-2-ene (0.20; CH₃CHOO) ≈ 2-methylpropene (0.20; 0.13 CH₂OO +0.07 (CH₃)₂COO). Beyond overall yield, sCI speciation varies considerably: aside from cis- and trans-but-2-ene, CH₂OO is the dominant sCI across all remaining alkene systems. Critically, the bimolecular reactivities of sCIs are strongly structure-dependent. With respect to SO₂ oxidation, CH₃CHOO, C₂H₅CHOO, (CH₃)₂COO, and (C(CH₃)=CH₂)CHOO exhibit rate coefficients three to four times higher than those of CH₂OO and (CH=CH₂)(CH₃)COO. The reactivity order toward H₂O/(H₂O)₂ is markedly different, however, broadly following CH₃CHOO ≈ C₂H₅CHOO > CH₂OO > (CH=CH₂)(CH₃)COO ≈ (CH₃)₂COO ≫ (C(CH₃)=CH₂)CHOO. These contrasts in sCI fate, including ozonolysis kinetics, effective sCI yield (i.e., the fraction of sCIs that persist to undergo bimolecular reactions), and the relative importance of sCI+SO₂ vs. sCI+H₂O/(H₂O)₂ pathways, collectively shape the relationship between alkenes and $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$. Hence, the capacity of different alkene ozonolysis systems to oxidize SO₂ cannot be reduced to ozonolysis kinetics and sCI yield alone; the competitive removal of sCIs by H₂O/(H₂O)₂ should also be considered (Fig. 2b). A previous H₂SO₄ proxy study similarly noted that the apparent rate constant describing the alkene ozonolysis source spans nearly three orders of magnitude in the real atmosphere, reflecting the underlying complexity of sCI chemistry (Yang et al., 2021).

Figure 3 SHAP dependence plots for O₃ and the six alkene precursors – ethene (a), propene/but-1-ene (b), 2-methylpropene (c), cis-but-2-ene (d), trans-but-2-ene (e), and isoprene (f) – in the XGBoost surrogate models of $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$. The x-axis shows the feature value, the y-axis shows the SHAP value of the same feature, and point color represents RH. Each point corresponds to one simulation scenario. These plots illustrate the marginal effect of each alkene precursor on the predicted sCI-driven SO₂ oxidation rate and how this effect varies with humidity.

Figure3 further shows the effects of various alkenes on $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ at different RH under MCM v3.3.1g, illustrating the regulatory role of the sCI+H₂O/(H₂O)₂ pathway. Consistent with the average |SHAP| ranking, isoprene, trans-but-2-ene, cis-but-2-ene, and 2-methylpropene all exhibit substantially positive impacts of varying degrees as their concentrations increase, whereas ethene and propene/but-1-ene display near-zero or weakly positive dependence across the investigated concentration range. The substantial vertical dispersion at fixed alkene concentrations, observed for all alkenes, indicates that identical alkene abundances can yield markedly different impacts on $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$, signifying the non-linear responses arising from interaction effects. For most alkenes (excluding propene/but-1-ene), this vertical spread increases with concentration, indicating that higher alkene concentrations expand the range of possible $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ values, while the realized magnitude remains strongly contingent on ambient conditions. Propene/but-1-ene, by contrast, shows large vertical dispersion across the full concentration range without a significant concentration dependence, suggesting that variability in its SHAP values is driven almost entirely by interaction effects. Ethene similarly shows a significant degree of dispersion in its SHAP values at low concentrations. Taken together, the relationship between alkene concentration and $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ is nonlinear for all alkenes; however, the degree of nonlinearity differs among alkenes, reflecting differences in the relative strength of their interactions with other feature variables.

We use RH coloring in Fig. 3 to visualize the strength of humidity's interactive effects. Across all alkenes, higher RH consistently suppresses the contributions of alkenes to $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$. The degree of suppression, however, varies markedly with alkene identity, and can be rationalized directly from the kinetic properties of their produced sCIs. Under high-humidity conditions, the contribution of propene/but-1-ene to the variation of $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ is weak. This is attributable to the relatively low rate coefficients for the reactions of propene/but-1-ene with O₃, which result in significantly lower sCI production than that from ethene and isoprene. The primary sCIs generated from propene/but-1-ene ozonolysis are CH₂OO, CH₃CHOO, and C₂H₅CHOO, all of which react rapidly with H₂O/(H₂O)₂. Additionally, CH₂OO, which constitutes the largest proportion, has relatively low reactivity toward SO₂. Consequently, the limited sCIs produced from these two alkenes make a further restricted contribution to sCI-mediated SO₂ oxidation. Ethene, despite having the slowest ozonolysis rate coefficient, produces CH₂OO with a substantially higher yield (0.42) than propene or but-1-ene, partially compensating for its sluggish kinetics. However, since CH₂OO is efficiently scavenged by H₂O/(H₂O)₂, ethene's contribution remains small and sensitive to humidity, consistent with the mixed positive and negative SHAP values observed at low concentrations in Fig. 3a. In contrast, cis-but-2-ene and trans-but-2-ene benefit from both the highest ozonolysis rate coefficients among the investigated alkenes and the exclusive production of CH₃CHOO, which reacts rapidly with SO₂. These advantages translate into strong positive contributions to $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$; however, because CH₃CHOO is also highly reactive toward H₂O/(H₂O)₂, the realized contribution is strongly RH-dependent, as evidenced by the pronounced RH stratification in the dependence plots. 2-methylpropene generates both CH₂OO and (CH₃)₂COO. (CH₃)₂COO reacts rapidly with SO₂ but comparatively slowly toward H₂O/(H₂O)₂, resulting in weak competition from H₂O/(H₂O)₂. This explains why 2-methylpropene shows the smallest vertical dispersion growth with concentration among the more influential alkenes, and why its relationship with $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ is more nearly linear, since (CH₃)₂COO maintains its reactivity toward SO₂ even under elevated humidity. Isoprene exhibits the strongest positive contribution to $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ across the full humidity range, arising from the favorable combination of relatively fast ozonolysis, the highest total sCI yield (0.66), and a diverse sCI speciation that includes not only CH₂OO but also (CH=CH₂)(CH₃)COO and (C(CH₃)=CH₂)CHOO. (C(CH₃)=CH₂)CHOO reacts rapidly with SO₂ while undergoing only slow scavenging by H₂O/(H₂O)₂, sustaining efficient SO₂ oxidation even at high RH. The interaction between isoprene and RH is thus primarily attributable to the CH₂OO fraction, whereas (CH=CH₂)(CH₃)COO and (C(CH₃)=CH₂)CHOO enable isoprene to remain a dominant driver of $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ across a wider humidity range. Taken together, when the investigated alkenes are perturbed over the same concentration ranges, ethene and propene/but-1-ene exert negligible direct control on $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ relative to the remaining alkenes. Isoprene and the C₄ alkenes, cis-/trans-but-2-ene and 2-methylpropene, represent the dominant alkene precursors, but their effective impacts are jointly shaped by RH. This finding underscores the need for surrogate modeling to resolve these complex, nonlinear atmospheric reactions into a clear relationship between precursors and $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$. Such modeling enables the identification of key alkene precursors, quantification of the dependence of their contributions on ambient environmental conditions, and elucidation of the key reactions governing their roles in sCI-mediated SO₂ oxidation.

3.2 Controls on the sCIs fractional contribution and regime-dependent sensitivities of SO₂ oxidation

The surrogate model developed in Sect. 3.1 quantifies how key factors govern $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$, the absolute rate of SO₂ oxidation by sCIs. A large $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$, however, does not necessarily imply that sCIs are important for SO₂ oxidation, because their importance also depends on the strength of the parallel OH-driven oxidation pathway. OH and sCIs share certain precursors, but their formation pathways are not entirely identical: alkene ozonolysis simultaneously produces both, so alkenes and O₃ drive both pathways, yet OH is additionally sustained through routes that do not involve sCIs, including RO_x–NO_x radical cycling and the direct photolysis of O₃, HONO, and other species. OH is therefore governed by a broader set of factors, including NO_x, NO₂ %, and the wider VOC composition, several of which exert little or no direct influence on sCI chemistry. As a result, the two oxidation pathways do not respond in tandem to environmental perturbations. Variations in the relative importance of sCIs therefore shape the sensitivity of $L_{{\mathrm{SO}}_{\mathrm{2}}}$ to its controlling factors. To systematically diagnose the atmospheric conditions under which sCIs play a significant role in SO₂ oxidation, and to elucidate how this role modulates the precursor sensitivity of $L_{{\mathrm{SO}}_{\mathrm{2}}}$, we constructed dedicated surrogate models for μ_sCIs and $L_{{\mathrm{SO}}_{\mathrm{2}}}$. Their feature sets capture the key controls on both pathways (see Sect. 2.3), and the models were trained on a large ensemble of constrained AtChem simulations spanning atmospheric conditions representative of anthropogenically and biogenically influenced environments. Because photolysis exerts decisive control over OH production, daytime (08:00–18:00 LT) and nighttime (22:00–05:00 LT) conditions were analyzed separately.

Figure 4 (a) First-order (S₁) and total-order (S_T) Sobol sensitivity indices for the seven input features in the XGBoost surrogate model of μ_sCIs during daytime. S₁ quantifies the main effect of each feature acting alone, whereas S_T includes both the main effect and all interaction effects involving that feature. The difference between S_T and S₁ therefore reflects the extent to which a feature participates in non-additive interactions. (b–i) Two-dimensional partial dependence plots showing the joint effect of O₃norm×NO_xnorm, VOCsnorm×NO_xnorm, O₃norm×NO₂ %norm, VOCsnorm×NO₂ %norm, O₃norm×alkene %norm, VOCsnorm×alkene %norm, VOCsnorm×ARIcode, and RHnorm×alkene %norm on the surrogate-model prediction of μ_sCIs during daytime. The surface represents the average model response after marginalizing over the remaining features. Warmer colors indicate larger predicted μ_sCIs.

Sobol sensitivity analysis was applied to identify the dominant atmospheric controls on μ_sCIs and to quantify the extent to which each feature exerts its influence independently or through interaction with other features during the daytime period (Fig. 4a). Here, the first-order index S₁ measures the independent main-effect contribution of a single feature, whereas the total-order index S_T additionally encompasses all higher-order interactions involving that feature; a markedly larger S_T relative to S₁ therefore indicates that a feature operates predominantly through coupling with other features rather than through its own variation alone. The analysis identifies ARI, O₃, alkene %, and VOCs as the dominant controls on μ_sCIs, with independent contributions ranked as ARI (S₁=0.17) > alkene % (S₁=0.16) > O₃ (S₁=0.15) > VOCs (S₁=0.13). RH ranks immediately behind (S₁=0.09), while NO_x and NO₂ % register the weakest direct contributions (S₁=0.03 and 0.04, respectively). The low-ARI tier (ARI =0) is dominated by ethene and propene, representative of typical urban anthropogenic emissions (Zhao et al., 2020); the mid-ARI tier (ARI =1) features an elevated proportion of but-2-ene, characteristic of industrial source environments or marine ecosystems (Giorio et al., 2022); and the high-ARI tier (ARI =2) is driven primarily by isoprene, representative of biogenic (Rhew et al., 2017) or petrochemical-influenced conditions (Guo et al., 2022). The alkene mixture proportions at each ARI level are detailed in Fig. S1. Comparing S₁ and S_T across all features reveals a clear differentiation in how each feature exerts its influence. For RH, S_T≈S₁ (S_T=0.11), indicating negligible interaction effects; RH influences μ_sCIs almost entirely through its independent main effect. For ARI (S_T=0.27), alkene % (S_T=0.26), O₃(S_T=0.22), and VOCs (S_T=0.23), S_T substantially exceeds S₁, indicating that although their independent main effects dominate, non-negligible interaction effects are also present. Among these, ARI exhibits both the largest S₁ and S_T, reinforcing that alkene speciation not only exerts the strongest independent control on μ_sCIs but also acts as a key moderator of the effects of O₃, alkene %, and VOCs. The contrast is most pronounced for NO_x and NO₂ %: their S₁ values account for less than half of their respective S_T values (NO_x: S_T=0.09; NO₂ %: S_T=0.09), indicating that their influence on μ_sCIs is overwhelmingly mediated through interactions with other features rather than through direct action. This strong interaction dependence reflects the fact that the influence of NO_x and NO₂ % on μ_sCIs is conditional on the values of other features.

The directional effects of these features are further illustrated by the 2D PDPs (Fig. 4b–i). Hereafter, feature levels denote normalized values rather than the actual concentrations; corresponding time-varying concentrations are shown in Fig. S1. Among the features that directly govern $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$, ARI, O₃, VOCs, and alkene % exhibit mutually reinforcing positive effects on μ_sCIs, suggesting that conditions favorable for alkene ozonolysis simultaneously enhance the relative competitiveness of the sCI pathway. RH negatively affects μ_sCIs by suppressing the sCI+SO₂ pathway through enhanced sCI loss to H₂O/(H₂O)₂. Despite NO₂ being a known sCI scavenger, NO₂ % exhibits a synergistic co-enhancement with O₃ and VOCs on μ_sCIs. These synergistic effects are expected to arise from the influence of NO₂ % on OH, through its role in governing how NO_x is partitioned between radical termination and cycling (Sillman and He, 2002). A higher NO₂ fraction, at fixed total NO_x, is expected to reduce the NO available for HO₂/RO₂→OH conversion, thereby suppressing the SO₂+OH oxidation channel more strongly than it scavenges sCIs, particularly under high-O₃ and high-VOCs conditions. Correspondingly, the higher the NO₂ %, the more strongly μ_sCIs increase with rising O₃ and VOCs. NO_x, by contrast, exerts predominantly suppressive interactions with O₃ and VOCs: μ_sCIs peaks at the lowest NO_x levels and declines sharply as NO_x increases only slightly from its minimum, before leveling off or increasing slightly at higher concentrations, such that μ_sCIs at the highest NO_x level remains well below the peak observed at the lowest NO_x level. The sensitivity of μ_sCIs to O₃ and VOCs is modulated by NO_x levels. At the lowest NO_xlevel, increases in O₃ and VOCs produce the largest increments in μ_sCIs. Across the remainder of the NO_x range, raising O₃ or VOCs from their minimum to maximum levels increases μ_sCIs by only 0.03–0.04. The effect of NO_x on μ_sCIs arises primarily from its control over the steady-state OH concentration: at extremely low NO_x level, radical recycling efficiency is greatly diminished, leading to suppressed OH; at excessively high NO_x levels, OH is likewise suppressed, but through accelerated termination reactions between radicals and NO_x. Under both conditions, the attenuation of OH-mediated SO₂ oxidation shifts a greater share of the total oxidation burden to sCIs, resulting in elevated μ_sCIs.

Taken together, daytime μ_sCIs tends to be elevated under atmospheric conditions characterized by high O₃ concentrations, abundant VOCs, a large alkene fraction, and relatively low RH. With respect to NO_x, μ_sCIs is only appreciably enhanced under extremely low NO_x levels; across the remainder of the NO_x range, its influence on μ_sCIs is comparatively weak. This finding is consistent with previous studies reporting that sCI chemistry plays a more prominent role in SO₂ oxidation in biogenic-rich forested and rural environments (Kim et al., 2015) than in polluted urban settings (Dada et al., 2020). Among the identified controls, O₃, VOCs, alkene %, and ARI, which collectively characterize conditions favorable for sCI formation, not only directly promote $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ but also exhibit pronounced synergistic positive effects on μ_sCIs. This pattern implies that, as these variables increase, OH increases less strongly than sCI production and may even decrease modestly. Consequently, although O₃ and alkenes jointly initiate the chemical processes leading to the formation of both OH and sCIs, their effects on the OH-mediated and sCI-mediated SO₂ oxidation pathways differ in both magnitude and direction. Under most daytime conditions, the OH pathway dominates SO₂ oxidation, such that the sensitivity of $L_{{\mathrm{SO}}_{\mathrm{2}}}$ to precursor species is primarily shaped by how those precursors affect OH. However, as shown by the NO_x dependence of μ_sCIs, the response of OH to precursor perturbations is itself condition-dependent, exhibiting a non-monotonic relationship with NO_x. This gives rise to distinct OH response patterns under high- and low-μ_sCIs conditions and, by extension, distinct sensitivities of $L_{{\mathrm{SO}}_{\mathrm{2}}}$ to its precursors. To quantify the magnitude of these differences, the dataset was partitioned into high- and low-μ_sCIs subsets based on the median daytime μ_sCIs value (1.1 %), and the relationship between each feature and $L_{{\mathrm{SO}}_{\mathrm{2}}}$ was characterized separately for each subset.

Figure 5 Individual conditional expectation (ICE) plots for the regime-specific $L_{{\mathrm{SO}}_{\mathrm{2}}}$ models during daytime, showing the responses of predicted total gas-phase SO₂ oxidation to VOCsnorm, alkene %norm, and O₃norm in the low- and high-μ_sCIs regimes. The two rows correspond to the two μ_sCIs regimes, and the three columns correspond to the three selected predictors. Each thin curve represents one sampled scenario, colored by normalized NO_x. The curves are displayed in absolute prediction space (i.e. without centering), so both the slope and vertical spread reflect heterogeneity in the conditional model responses across different chemical backgrounds.

Figure 5 illustrates how total gas-phase SO₂ oxidation rate varies with O₃, VOCs, and alkene % across the high- and low-μ_sCIs regimes. The mean ICE curves show that $L_{{\mathrm{SO}}_{\mathrm{2}}}$ increases monotonically with all three features under both regimes; however, the response is systematically amplified under high-μ_sCIs conditions, with the mean curve exhibiting a steeper and more sustained increase across the full feature range. Beyond the mean response, the spread of individual ICE curves provides insight into response heterogeneity, namely, the extent to which the effect of a given feature on $L_{{\mathrm{SO}}_{\mathrm{2}}}$ varies across different atmospheric conditions. Notably, negative $L_{{\mathrm{SO}}_{\mathrm{2}}}$ responses to increases in VOCs and alkene % are observed under extremely low NO_x level. To quantify this directional heterogeneity, we computed the fraction of individual ICE curves exhibiting a net negative slope between the minimum and maximum feature values, termed the “neg-slope fraction”. For O₃, positive associations with $L_{{\mathrm{SO}}_{\mathrm{2}}}$ are nearly universal in both regimes, with a neg-slope fraction of 0 % under the high-μ_sCIs regime and 9 % under the low-μ_sCIs regime. For alkene %, the neg-slope fraction increases from 7 % under the high-μ_sCIs regime to 21 % under the low-μ_sCIs regime. A similar pattern is evident for VOCs, with the neg-slope fraction rising from 21 % under high-μ_sCIs to 35 % under low-μ_sCIs regimes. This regime contrast warrants further examination. At an extremely low NO_x level, high μ_sCIs is associated with abundant alkene concentrations. Although higher alkene concentrations may promote radical termination and thereby suppress OH and reduce $L_{{\mathrm{SO}}_{\mathrm{2}}}$, the high-μ_sCIs subset nevertheless exhibits a smaller neg-slope fraction than the low-μ_sCIs subset. This suggests that the alkene ozonolysis pathway offsets part of the OH suppression caused by enhanced radical termination, thereby dampening the negative $L_{{\mathrm{SO}}_{\mathrm{2}}}$ response. At moderate-to-high NO_x levels, positive $L_{{\mathrm{SO}}_{\mathrm{2}}}$ responses to increases in VOCs, alkene %, and O₃ dominate and are considerably stronger under the high-μ_sCIs regime than under the low-μ_sCIs regime. This amplified sensitivity likely reflects the fact that high-μ_sCIs regime is associated with abundant alkenes, which simultaneously promote OH propagation through RO_x cycling pathway and drive alkene ozonolysis.

Figure 6 Two representative daytime scenarios illustrating the budgets of OH and sCIs at normalized NO_x levels of (a) 0 and (b) 0.5. Numbers denote reaction rates (10⁶ molec. cm⁻³ s⁻¹). Pink, green, and purple lines indicate the termination, initiation, and propagation pathways of OH and sCIs, respectively.

In summary, $L_{{\mathrm{SO}}_{\mathrm{2}}}$ exhibits both positive and negative responses to VOCs and alkene % across the full scenario space, with the key distinction being that under the high-μ_sCIs regime, negative responses are less frequent and positive responses are stronger. To clarify the mechanistic basis of this contrast, we examined the OH and sCI budgets for two representative scenarios at normalized NO_x levels of 0 and 0.5 (Fig. 6a, b). This analysis was used to explain why negative $L_{{\mathrm{SO}}_{\mathrm{2}}}$ responses to increases in VOCs and alkene % are more frequent in the low-μ_sCIs regime. OH production is partitioned into initiation (∑OH_new), which encompasses all primary OH production via photolysis and alkene ozonolysis (Sheehy et al., 2010), and propagation (OH_propag), defined as the rate of HO₂→OH conversion. At an extremely low NO_x level (NO_x, norm = 0), the ozonolysis pathway accounts for 39.65 % of ∑OH_new (see Sect. S5 for detailed calculations), and the ratio OH_propag/∑OH_new is 3.67, indicating that radical propagation contributes approximately 3.7 times as much OH as primary initiation. At a moderate NO_x level (NO_x, norm =0.5), the ozonolysis share of ∑OH_new decreases to 29.28 %, whereas the OH_propag/∑OH_new ratio rises to 6.96, reflecting substantially enhanced radical cycling efficiency, as evidenced by the HO₂→OH conversion rate reaching 7.04×10⁸ molec. cm⁻³ s⁻¹. At an extremely low NO_x level, the RO_x recycling is severely suppressed, thus primary OH initiation contributes a correspondingly larger fraction to the total OH budget, and the ozonolysis pathway in turn accounts for a larger share of OH production. Under this condition, the contribution of alkene ozonolysis to OH production can partly offset the suppressive effect of enhanced radical termination on OH, which explains why the neg-slope fraction for $L_{{\mathrm{SO}}_{\mathrm{2}}}$ is smaller in the high-μ_sCIs subset than in the low-μ_sCIs subset. At a moderate NO_x level, OH is sustained predominantly through radical cycling rather than through photolysis or ozonolysis, and the increase in NO_x additionally enhances HONO photolysis, further diminishing the relative contribution of ozonolysis to OH generation. Under these conditions, perturbations in VOCs and alkene % influence $L_{{\mathrm{SO}}_{\mathrm{2}}}$ primarily through their modulation of OH propagation. Because high-μ_sCIs conditions at this NO_x level are associated with abundant alkenes, both the RO_x cycling and ozonolysis pathways are strongly activated, and $L_{{\mathrm{SO}}_{\mathrm{2}}}$ therefore exhibits a pronounced positive response to increases in precursor concentrations.

Figure 7 (a) First-order (S₁) and total-order (S_T) Sobol sensitivity indices for the seven input features in the XGBoost surrogate model of μ_sCIs during nighttime. S₁ quantifies the main effect of each feature acting alone, whereas S_T includes both the main effect and all interaction effects involving that feature. The difference between S_T and S₁ therefore reflects the extent to which a feature participates in non-additive interactions. (b–i) Two-dimensional partial dependence plots showing the joint effect of O₃norm×NO_xnorm, VOCsnorm×NO_xnorm, O₃norm×NO₂ %norm, VOCsnorm×NO₂ %norm, O₃norm×alkene %norm, VOCsnorm×alkene %norm, VOCsnorm×ARIcode, and RHnorm×alkene %norm on the surrogate-model prediction of μ_sCIs during nighttime. The surface represents the average model response after marginalizing over the remaining features. Warmer colors indicate larger predicted μ_sCIs.

The nighttime Sobol sensitivity indices (Fig. 7a) reveal controls on μ_sCIs that differ markedly from those identified during the day. ARI, RH, and NO_x emerge as the dominant controls, with independent contributions ranked as ARI (S₁=0.21) > RH (S₁=0.19) > NO_x (S₁=0.17). VOCs (S₁=0.13), alkene % (S₁=0.11), and O₃ (S₁=0.08) rank immediately behind, while NO₂ % exerts a negligible direct influence (S₁=0.002). Comparing S₁ and S_T across all features, S_T only slightly exceeds S₁ for most features, indicating that independent main effects dominate and interaction effects are limited, with the total-order ranking following ARI (S_T=0.25) > RH (S_T=0.22) > NO_x (S_T=0.20) > VOCs (S_T=0.17) > alkene % (S_T=0.14) > O₃ (S_T=0.10). NO₂ % retains the smallest S_T (S_T=0.02), reinforcing that its regulatory effect on μ_sCIs at night is negligible. The close agreement between S_T and S₁ across most features at night stands in sharp contrast to the daytime results, where S_T substantially exceeded S₁ for several key features. This collapse of interaction effects reflects a fundamental mechanistic shift: in the absence of photolysis, the primary OH initiation pathway is alkene ozonolysis, which is simultaneously the sole sCI formation pathway. Perturbations to precursors therefore affect both oxidation pathways in the same direction and with similar magnitude, suppressing the conditional interactions that characterize the daytime sensitivity structure. This finding is broadly consistent with previous studies demonstrating that Criegee intermediate chemistry constitutes a non-negligible OH source under dark or low-light conditions (Khan et al., 2018). The rise of RH from a minor factor during the day to the second-ranked control at night reflects the substantially lower OH concentrations characteristic of nocturnal conditions, which elevate the relative contribution of sCIs to gas-phase SO₂ oxidation (Kukui et al., 2021) and thereby amplify the suppressive effect of H₂O/(H₂O)₂ on μ_sCIs. The contrasting behavior of NO_x and NO₂ % between day and night also warrants attention. At night, photolytic regeneration of NO from NO₂ is effectively absent; instead, NO within NO_x is oxidized to NO₂ through reactions with O₃ and peroxy radicals, leading to progressive NO₂ accumulation. This enhances the termination reaction between NO₂ and OH, suppressing steady-state OH and indirectly elevating μ_sCIs. As a result, nighttime μ_sCIs is governed primarily by total NO_x abundance rather than by the NO/NO₂ partitioning, which explains why NO₂ % exerts a negligible effect at night.

The directional effects of these features are further elucidated by the 2D PDPs (Fig. 7b–i). ARI, O₃, VOCs, and alkene % exhibit mutually reinforcing positive effects on μ_sCIs, consistent with their roles as sCI precursors: conditions favorable for $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ simultaneously enhance the relative competitiveness of the sCI oxidation pathway. RH exerts a pronounced negative effect on μ_sCIs. With respect to NO_x, μ_sCIs reaches a locally elevated value (approximately 0.1) at the lowest NO_x level, then drops sharply as NO_x increases only slightly from its minimum, before rising significantly and monotonically across the remainder of the range to reach its global maximum (approximately 0.16) at the highest NO_x level, such that a strong positive relationship between NO_x and μ_sCIs prevails across most of the range. This near-monotonic behavior reflects the role of NO_x in modulating steady-state OH through radical termination, as described above.

Taken together, nighttime μ_sCIs tends to be elevated under atmospheric conditions characterized by high O₃ concentrations, abundant VOCs, a large alkene fraction, and relatively low RH. With respect to NO_x, the highest μ_sCIs values occur under elevated NO_x levels, although a smaller enhancement is also evident under extremely low NO_x conditions. This is consistent with a previous study reporting that elevated NO_x levels in power plant environments at night are highly conducive to a large sCI contribution fraction (Meidan et al., 2019). Furthermore, ARI, O₃, VOCs, and alkene %, which collectively characterize conditions favorable for sCI formation, not only strongly promote $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ but also exert pronounced synergistic effects on μ_sCIs. This pattern indicates that O₃ and alkenes do not affect the OH- and sCI-mediated SO₂ oxidation pathways to the same extent. Notably, unlike during the day, when OH is sustained primarily by radical propagation, nocturnal OH production is more strongly linked to primary initiation, alkene ozonolysis, because the absence of photolysis limits OH production from photolytic sources and radical recycling. Alkene ozonolysis is also the sole formation pathway for sCIs. Therefore, differences in the sensitivity of $L_{{\mathrm{SO}}_{\mathrm{2}}}$ across μ_sCIs regimes are expected to reflect differences in the strength of the alkene ozonolysis pathway. To examine these differences, the dataset was partitioned into high- and low-μ_sCIs subsets based on the median nighttime μ_sCIs value (8.0 %), and the relationship between each feature and $L_{{\mathrm{SO}}_{\mathrm{2}}}$ was characterized separately for each subset.

Figure 8 Individual conditional expectation (ICE) plots for the regime-specific $L_{{\mathrm{SO}}_{\mathrm{2}}}$ models during nighttime, showing the responses of predicted total gas-phase SO₂ oxidation to VOCsnorm, alkene %norm, and O₃norm in the low- and high-μ_sCIs regimes. The two rows correspond to the two μ_sCIs regimes, and the three columns correspond to the three selected predictors. Each thin curve represents one sampled scenario, colored by normalized NO_x. The curves are displayed in absolute prediction space (i.e. without centering), so both the slope and vertical spread reflect heterogeneity in the conditional model responses across different chemical backgrounds.

Figure 8 shows how $L_{{\mathrm{SO}}_{\mathrm{2}}}$ responds to precursor perturbations across the high-μ_sCIs and low-μ_sCIs regimes at night. The mean ICE curves reveal that $L_{{\mathrm{SO}}_{\mathrm{2}}}$ increases monotonically with all three features under both regimes; however, the response is systematically enhanced under the high-μ_sCIs regime, with the mean curve exhibiting a steeper and more sustained increase across the full feature range. Beyond the mean response, the spread of individual ICE curves provides insight into response heterogeneity across different atmospheric conditions. For O₃, near-universal positive associations with $L_{{\mathrm{SO}}_{\mathrm{2}}}$ are observed in both regimes, with neg-slope fractions of 1 % and 3 % under the high- and low-μ_sCIs regimes, respectively, indicating that O₃ consistently promotes $L_{{\mathrm{SO}}_{\mathrm{2}}}$ regardless of the μ_sCIs regime. For alkene %, the neg-slope fraction increases from 1 % under the high-μ_sCIs regime to 10 % under the low-μ_sCIs regime. The contrast is most pronounced for VOCs, for which the neg-slope fraction rises from 5 % under the high-μ_sCIs regime to 29 % under the low-μ_sCIs regime, with the excess negative responses in the low-μ_sCIs regime arising predominantly from high-NO_x scenarios.

Figure 9 Two representative nighttime scenarios illustrating the budgets of OH and sCIs corresponding to (a) the low-μ_sCIs regime (μ_sCIs = 0.008) and (b) the high-μ_sCIs regime (μ_sCIs = 0.14). Numbers denote reaction rates (10⁶ molec. cm⁻³ s⁻¹). Pink, green, and purple lines indicate the termination, initiation, and propagation pathways of OH and sCIs, respectively.

To elucidate the role of alkene ozonolysis in driving these response differences, we examined the OH and sCI budgets for two representative high-NO_x scenarios corresponding to the low-μ_sCIs (μ_sCIs=0.008) and high-μ_sCIs (μ_sCIs=0.14) regimes, respectively (Fig. 9a, b). Under the low-μ_sCIs scenario, alkene ozonolysis accounts for 30.18 % of ∑OH_new, which is 5.56×10⁵ molec. cm⁻³ s⁻¹. In contrast, under the high-μ_sCIs scenario, alkene ozonolysis contributes 94.24 % of ∑OH_new, which reaches 2.29×10⁷ molec. cm⁻³ s⁻¹. These results show that, when alkene ozonolysis is efficient at night, it acts not only as the sole formation pathway for sCIs but also as the dominant OH initiation pathway. When alkene ozonolysis is weak, the ratio OH_propag/∑OH_new is 16.23, indicating that OH is sustained predominantly through radical cycling. In this regime, OH initiation is supplied primarily by photolysis in the early morning hours, with a smaller contribution from ozonolysis. Reducing total VOCs does not affect the photolytic OH initiation source, but decreases the OH–to–RO₂ flux and thus weakens OH loss through VOC oxidation. Provided that radical recycling efficiency changes little, this reduction in OH reactivity can slightly increase steady-state OH and $L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{OH}}$. Because the alkene-derived sCI flux is negligible in this regime, $\mathrm{\Delta }{L}_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ is close to zero. The OH-mediated response can therefore dominate, producing a net increase in $L_{{\mathrm{SO}}_{\mathrm{2}}}$ as VOCs decrease (a negative slope). Conversely, when alkene ozonolysis is strong, it simultaneously serves as the dominant OH initiation source and the principal sCI formation pathway. The ratio OH_propag/∑OH_new decreases to 2.68, far below that under weak ozonolysis, indicating a much greater relative contribution of initiation to OH. Under these conditions, VOC perturbations influence the OH- and sCI-mediated SO₂ oxidation pathways in the same direction, thereby reducing the occurrence of negative $L_{{\mathrm{SO}}_{\mathrm{2}}}$ responses. Between the two scenarios, the rates of ozonolysis-driven sCI and OH formation increase by factors of 127.5 and 123.3, respectively. Accordingly, the rate at which sCIs form H₂SO₄ increases by a factor of 133, whereas the rate at which OH forms H₂SO₄ increases by only a factor of six. This difference arises because OH reacts with alkenes and is incorporated into the RO_x cycle, such that part of the additional OH production supports radical cycling rather than directly increasing steady-state OH. The contrasting neg-slope fractions for VOCs and alkene % can thus be explained by their distinct chemical roles. Reducing alkene % simultaneously weakens a primary OH initiation source and the sCI formation pathway, such that $L_{{\mathrm{SO}}_{\mathrm{2}}}$ declines almost universally. By contrast, reducing total VOCs primarily lowers OH reactivity rather than suppressing OH initiation under alkene-poor conditions corresponding to the low-μ_sCIs regime. The resulting increase in OH-mediated oxidation can outweigh the negligible change in the sCI pathway, producing a negative $L_{{\mathrm{SO}}_{\mathrm{2}}}$ response (neg-slope =29 %).

Figure 10 Constrained AtChem–MCM v3.3.1 simulations for a summer episode in Wuhai City (1 June–15 July 2021), showing time series of H₂SO₄–AtChem (AtChem-simulated sulfuric acid concentration), μ_sCIs–AtChem (AtChem-simulated fractional contribution of sCIs to total SO₂ oxidation), SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (concurrently measured ambient particulate sulfate mass concentration), and μ_sCIs-surrogate model (daytime- and nighttime-averaged μ_sCIs predicted by the surrogate model).

3.3 Observational–model corroboration for regime-dependent sensitivities of SO₂ oxidation

To validate the predictive capability of the surrogate model developed in Sect. 3.2 for μ_sCIs, and to evaluate how the sensitivity of SO₂ oxidation to its precursors depends on the relative contributions of sCIs in ambient air, we selected a summer episode in Wuhai City as a real-atmosphere validation case. Figure 10 presents the constrained simulation results derived from AtChem–MCM v3.3.1, including time series of sulfuric acid concentrations (H₂SO₄–AtChem) and the daytime- and nighttime-averaged fractional contributions of sCIs to total SO₂ oxidation (μ_sCIs–AtChem) over the period 1 June to 15 July 2021. The corresponding pollutant mixing ratios and meteorological parameters are shown in Fig. S3. SO${}_{\mathrm{4}}^{\mathrm{2}-}$ denotes the ambient particulate sulfate mass concentrations measured concurrently by online ion chromatography. The term “μ_sCIs–surrogate model” represents the daytime- and nighttime-averaged μ_sCIs predictions generated by the surrogate model described in Sect. 3.2. The strong agreement (R²=0.75) between the surrogate predictions (μ_sCIs–surrogate model) and the explicit AtChem simulations (μ_sCIs–AtChem) demonstrates the robust predictive performance of the constructed machine learning surrogate for μ_sCIs.

Figure 11 (a) SHAP beeswarm summary plot for the observation-based XGBoost model of ambient SO${}_{\mathrm{4}}^{\mathrm{2}-}$. Each point represents one hourly observation, positioned according to the SHAP value of the corresponding predictor. Positive SHAP values indicate that the predictor increases the model-predicted SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentration, whereas negative values indicate a suppressing effect. Point colour denotes the predictor magnitude from low to high. Features are ordered by mean absolute SHAP value, so the plot summarizes both the relative importance and the directionality of the physically interpretable features. (b) Comparison of SHAP beeswarm plots for the observation-based XGBoost model under daytime low- and high-μ_sCIs regimes. Each point represents one hourly observation in the corresponding subset. The horizontal position gives the SHAP value, and the color indicates the feature value. Differences between the two panels illustrate how the relative importance and directional effects of observed predictors on ambient SO${}_{\mathrm{4}}^{\mathrm{2}-}$ vary between daytime chemical regimes with low and high inferred sCI contributions. (c) Comparison of SHAP beeswarm plots for the observation-based XGBoost model under nighttime low- and high-μ_sCIs regimes. Each point represents one hourly observation in the corresponding subset. The horizontal position gives the SHAP value, and the color indicates the feature value. Differences between the two panels illustrate how the relative importance and directional effects of observed predictors on ambient SO${}_{\mathrm{4}}^{\mathrm{2}-}$ vary between nighttime chemical regimes with low and high inferred sCI contributions.

To evaluate whether ambient SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations serve as a reliable proxy for gas-phase H₂SO₄ production during the target episode, an XGBoost model was trained on a comprehensive long-term observational dataset spanning 9 October 2019 to 30 June 2022 (summary statistics of all variables are provided in Table S4), with hourly SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentration as the target variable. The feature set included PM_2.5, NO${}_{\mathrm{3}}^{-}$, SO₂, O₃, VOCs, alkene %, NO_x, NO₂ %, RH, wind speed (WS), and Fe. Training over this extended period ensures that the model captures the full range of variability in sulfate formation, including seasonal cycles, episodic pollution events, and varying meteorological regimes, thereby establishing a representative and well-generalized mapping between precursor conditions and SO${}_{\mathrm{4}}^{\mathrm{2}-}$. SHAP values were then extracted exclusively for the target episode (1 June to 15 July 2021) to attribute feature contributions during this specific period. As shown in Fig. 11a, the three highest-ranked features are NO${}_{\mathrm{3}}^{-}$, SO₂, and PM_2.5, all of which exhibit clear positive correlations with SO${}_{\mathrm{4}}^{\mathrm{2}-}$. The elevated NO${}_{\mathrm{3}}^{-}$ co-occurs with SO${}_{\mathrm{4}}^{\mathrm{2}-}$, together with the concurrent increase in PM_2.5, consistent with the hallmark of secondary inorganic aerosol formation (Gao et al., 2021). The strong positive SHAP contribution of SO₂ directly corroborates its role as the primary precursor for SO${}_{\mathrm{4}}^{\mathrm{2}-}$ production. Notably, WS also exerts a discernible positive contribution, suggesting a degree of regional transport influence during this episode, although its impact remains secondary to in situ chemical production. By contrast, RH shows no significant positive correlation with SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and ranks relatively low in feature importance. Fe, a key catalyst for aqueous-phase SO₂ oxidation via Fenton-type chemistry (Ye et al., 2023), shows no clear positive association with SO${}_{\mathrm{4}}^{\mathrm{2}-}$, with elevated ambient Fe concentrations not corresponding to enhanced sulfate formation. Taken together, these results indicate that aqueous-phase oxidation was not the dominant SO${}_{\mathrm{4}}^{\mathrm{2}-}$ formation pathway during the target episode. Meanwhile, O₃, VOCs, and alkene %, which serve as key proxies of gas-phase oxidant pathways, exhibit positive correlations with SO${}_{\mathrm{4}}^{\mathrm{2}-}$, providing evidence that gas-phase oxidation constituted the primary sulfate formation pathway during this period. On this basis, ambient SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations can be regarded as a reliable observational proxy for the variability in total gas-phase H₂SO₄ production during the target episode.

To further examine whether the sensitivity of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ to key precursors is consistent with the regime-dependent behavior of $L_{{\mathrm{SO}}_{\mathrm{2}}}$ identified in Sect. 3.2, the target episode was partitioned into high- and low-μ_sCIs subsets for daytime and nighttime periods separately, using the same classification thresholds applied in Sect. 3.2 (daytime: 1.1 %; nighttime: 8.0 %). SHAP values for O₃, VOCs, and alkene % were extracted for each subset and compared. During the daytime (Fig. 11b), the high-μ_sCIs subset exhibits systematically amplified SHAP responses relative to the low-μ_sCIs subset across all three features. O₃ shows the most pronounced regime contrast: in the high-μ_sCIs subset, high feature values (red) are associated with distinctly larger positive SHAP contributions and a considerably wider spread, whereas the corresponding responses in the low-μ_sCIs subset are substantially attenuated. VOCs and alkene % follow the same pattern, with high feature values driving larger positive SHAP contributions under high-μ_sCIs conditions. This daytime regime contrast is consistent with the finding in Sect. 3.2 that the sensitivity of SO₂ oxidation to O₃, VOCs, and alkene % is amplified under high-μ_sCIs conditions, providing observational corroboration of that mechanistic result. At night (Fig. 11c), VOCs and alkene % exhibit stronger positive SHAP contributions to SO${}_{\mathrm{4}}^{\mathrm{2}-}$ in the high-μ_sCIs subset, consistent with their roles in promoting both sCIs and OH formation through nocturnal alkene ozonolysis. By contrast, O₃ exhibits negative SHAP contributions in both subsets, which can be attributed to the distinct behavior of O₃ at night: nocturnal O₃ is primarily of advective origin rather than in situ photochemical production, and elevated O₃ concentrations are often associated with suppressed local NO levels due to rapid O₃–NO titration. Taken together, the regime-dependent SO${}_{\mathrm{4}}^{\mathrm{2}-}$ sensitivity patterns observed during the target episode (Fig. 11b, c) provide observational validation for the mechanistic findings derived from the AtChem-based surrogate model simulations in Sect. 3.2, confirming that these findings are manifested in ambient atmospheric measurements. This cross-validation between model-diagnosed sensitivity and observation-based attribution reinforces the conclusion that under high-μ_sCIs regimes, alkene-targeted emission control strategies can effectively reduce sCIs and suppress the total SO₂ oxidation rate, while rarely triggering a rebound in OH-driven oxidation pathways.

4 Summary and conclusion

Sulfuric acid plays a central role in atmospheric aerosol formation, yet effective mitigation requires quantifying the sensitivity of H₂SO₄ formation to its precursors, which varies across atmospheric conditions as the relative contributions of different gas-phase SO₂ oxidation pathways shift. By treating sCIs as reactive intermediate species linking precursor emissions to H₂SO₄ formation, this study deploys an interpretable multi-target XGBoost modeling framework trained on ensemble box-model simulations incorporating updated CI chemistry. This framework determines the controls on sCI-mediated SO₂ oxidation ($L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$) and on the fractional sCI contribution (μ_sCIs), while establishing how the regime of μ_sCIs itself reshapes the sensitivity of total SO₂ oxidation ($L_{{\mathrm{SO}}_{\mathrm{2}}}$) to its precursors.

$L_{{\mathrm{SO}}_{\mathrm{2}},\mathrm{sCIs}}$ is dominated by RH (29.5 %) and isoprene (23.8 %), revealing that low humidity and high biogenic emissions primarily promote the sCI-mediated oxidation pathways. For μ_sCIs, the dominant controlling factors exhibit a pronounced diurnal shift: daytime control resides in alkene reactivity (ARI), O₃, alkene fraction, and total VOCs, whereas nighttime control transitions to ARI, RH, and NO_x. High-μ_sCIs regimes uniformly amplify positive $L_{{\mathrm{SO}}_{\mathrm{2}}}$ sensitivities to O₃, VOCs, and alkene %, while decreasing the frequency of occasional negative $L_{{\mathrm{SO}}_{\mathrm{2}}}$ responses to VOCs or alkene %. Observation-based modeling of a Wuhai summer episode with SHAP-based attribution of ambient sulfate provided support for this regime dependence.

Previous studies have extensively reported the contribution of the SO₂ + sCI reaction to total sCI loss (Cox et al., 2020), and the contribution of sCIs to overall SO₂ oxidation across diverse atmospheric conditions (Khan et al., 2018). Building on the understanding that sCI chemistry constitutes one pathway for SO₂ oxidation and that sCIs are secondarily formed atmospheric intermediates, this study no longer treats them as terminal oxidants but rather takes seriously their role as an intermediate species and establishes quantitative relationships linking controlling factors to both the sCI-mediated SO₂ oxidation rate and its fractional contribution to total SO₂ oxidation. The resulting picture goes beyond estimating the sCI contribution for an individual scenario, extending instead to a predictive understanding of the role sCIs may play across a wider range of atmospheric conditions. By further analyzing how the sensitivity of total SO₂ oxidation to factor perturbations varies among internal oxidation regimes characterized by different μ_sCIs values, this work provides a scientific basis for developing strategies to mitigate H₂SO₄ formation. Beyond these mechanistic-sensitivity insights, the interpretable surrogate modeling framework developed here offers a transferable approach for distilling quantitative relationships between reactants and products from complex reaction networks. By linking targeted box-model scenario design to machine-learning-based sensitivity diagnostics, this methodology enables efficient identification of dominant kinetic controls – providing a basis for prioritizing kinetic processes for explicit representation when developing lumped chemical mechanisms for chemical transport models.

It should be noted that applying these findings to real-world atmospheric phenomena requires recognizing specific applicability conditions. While key feature variables were explicitly modeled, capturing full atmospheric variability also requires future integration of transport, deposition, multiphase processes, temperature dependence, and other factors. Additionally, our 100 % SO₃ yield assumption for the sCI + SO₂ reaction establishes an upper bound for absolute sCI-mediated H₂SO₄ production and fractional contributions. Because this yield functions as a linear scaling factor, absolute magnitudes may be overestimated, but the underlying functional relationships between controlling factors and sCI chemistry remain unchanged.