Driving factors of oxalic acid and enhanced role of gas-phase oxidation under cleaner conditions: insights from 2007–2018 field observations in the Pearl River Delta

Secondary organic aerosol (SOA) is a dominant constituent of fine particulate matter, exerting significant impacts on both climate and human health. Oxalic acid (C₂), a key end-product formed from the oxidation of volatile organic compounds, can provide insights into the formation mechanism of SOA. Thus, long-term measurements of C₂ and related compounds help understand the changes in SOA formation with decreasing pollutant levels. In this study, C₂ and its homologs, along with five primary anthropogenic source markers and three SOA markers, were measured in the Pearl River Delta (PRD) during 2007–2018. The concentrations of C₂ did not exhibit significant downward trends, despite substantial reductions in anthropogenic emissions, such as biomass burning (−11 % yr⁻¹), vehicle emissions (−17 % yr⁻¹), and cooking emissions (−7 % yr⁻¹). Correlation analysis revealed that aerosol liquid water content (ALWC) and O_x (O₃ + NO₂) were the main drivers of C₂ variations. Moreover, the relative contribution of biogenic SOA increased under cleaner conditions. A machine learning model was applied to quantify the impacts of changes in anthropogenic precursor emissions, biogenic precursor emissions, aqueous-phase oxidation processes, and gas-phase oxidation processes on C₂ variability. As pollution levels declined, the impacts of gas-phase oxidation increased from 37 % to 55 %, whereas that of aqueous-phase oxidation declined from 42 % to 30 %. This shift indicated a transition from aqueous-phase to gas-phase pathways in C₂ and SOA formation. Our findings highlight the increasing importance of gas-phase oxidation under low-pollution conditions and underscore the need for effective ozone control strategies to further reduce SOA in the future.

1 Introduction

Low molecular weight dicarboxylic acids (DCA) represent one of the most abundant classes of organic compounds in aerosols and are ubiquitously distributed in the lower troposphere (Kawamura and Bikkina, 2016). In urban areas, DCA and related compounds (e.g., glyoxal (Gly), methylglyoxal (mGly), glyoxylic acid (ωC₂), and pyruvic acid (Pyr)) account for 0.2 %–6 % of total carbon in aerosols (Kawamura and Ikushima, 1993; Kawamura and Watanabe, 2004), while this proportion could rise to 16 % in remote marine (Kawamura and Sakaguchi, 1999; Wang et al., 2006). Due to their high solubility and hygroscopicity (Bilde et al., 2015), DCA can promote aerosol hygroscopic growth and facilitate the activation of cloud condensation nuclei (CCN). Thus, they play a crucial role in the Earth's climate and radiative forcing by scattering solar radiation and forming clouds (Kawamura et al., 2013a; Meng et al., 2018; 2023).

Previous studies have suggested that DCA can be directly emitted from biogenic sources (Rinaldi et al., 2011), fossil fuel combustion (Kawamura and Kaplan, 1987), biomass burning (Kawamura et al., 2013b; Cao et al., 2017), waste incineration (Jung et al., 2010), and cooking (He et al., 2004). However, increasing evidence has demonstrated that most DCA are generated through photochemical oxidation of volatile organic compounds (VOCs), such as olefins (Bikkina et al., 2021), cyclic alkanes (Kawamura and Usukura, 1993), and unsaturated fatty acids (Kawamura et al., 1996), followed by partitioning into wet aerosols. Therefore, DCA have been widely used as SOA indicators to investigate the secondary reactions and aging processes of VOCs in the atmosphere (Kawamura and Bikkina, 2016; Meng et al., 2023). As a key end-product of numerous formation pathways in aerosols, oxalic acid (C₂) is the most abundant compound among DCA homologs (58 %–89 %) (Ho et al., 2007; Xu et al., 2023). Extremely high concentrations of C₂ (∼ 1000 ng m⁻³) have been observed in polluted areas (Poore, 2000; Wang et al., 2012; Xu et al., 2023), suggesting fast photochemical oxidation of organic precursors. Thus, a growing number of studies have focused on identifying the sources and elucidating the formation mechanism of C₂ over the past decades.

Previous studies have shown that the concentrations of DCA and C₂ in urban areas are generally higher than those in coastal, mountain, and marine areas (Table S1 in the Supplement). The elevated concentrations have been attributed to higher emissions of organic precursors in urban environments with intensive human activities. However, an unexpected increase in C₂ was reported during the COVID-19 period (Meng et al., 2023), a period characterized by a dramatic reduction in anthropogenic emissions. They identified higher oxidant levels and stronger solar radiation as drivers of this C₂ enhancement. These findings indicate that effective mitigation of C₂ and SOA requires the simultaneous control of both anthropogenic emissions and atmospheric oxidants. As anthropogenic emissions decline, the relative contributions of biogenic source to C₂ are expected to increase. For example, Xu et al. (2022) found that biogenic sources accounted for only 22 % of C₂ in continental outflows characterized by high anthropogenic emissions. In contrast, this proportion increased to 61 % in oceanic air masses with minimal anthropogenic influence. In addition, relatively low anthropogenic emissions and enhanced biogenic activity in summer led to a higher non-fossil contribution to C₂ (46 %) compared to that in winter (33 %) (Xu et al., 2023).

Numerous studies have reported strong correlations between C₂ and sulfate across diverse regions, including East Asia (Yu et al., 2005; Wang et al., 2017; Zhang et al., 2022a; Meng et al., 2023), the USA (Hilario et al., 2021), the Pacific and Atlantic Oceans (Hilario et al., 2021), the Bay of Bengal (Bikkina et al., 2017) and Mongolia (Jung et al., 2010). This widespread correlation suggests that C₂ and sulfate may share common formation pathways across different atmospheric environments. Assuming that sulfate primarily originates from aqueous-phase processing, C₂ has been considered a secondary product formed via aqueous reactions (Yu et al., 2005; Hilario et al., 2021). However, recent research found that gas-phase oxidation also played an important role in C₂ formation. For instance, Ding et al. (2021) observed a strong correlation between C₂ and ozone (O₃) but not with sulfate during a haze event in the Yangtze River Delta, indicating C₂ was mainly produced via gas-phase oxidation. During COVID-19, lower aerosol liquid water content (ALWC) and elevated O₃ shifted the dominant formation pathway of C₂ from aqueous-phase oxidation of ωC₂ and Pyr to gas-phase photochemical decomposition of longer-chain DCA (malonic (C₃) and succinic (C₄)). As a result, the contribution of gas-phase oxidation increased from 12.3 % to 50.5 % (Meng et al., 2023). In addition, Xu et al. (2023) found that the gas-phase photochemical aging dominated the formation of C₂ at five emission hotspots in China during summer, when pollutant levels were relatively low and solar radiation and oxidant concentrations were high. Recently, a laboratory experiment also proposed new reaction schemes for high levels of C₂ via gas-phase chemical reactions with O₃ and hydroxyl radicals, which could explain the high levels of C₂ in the free troposphere (Bikkina et al., 2021). Despite these findings, the formation mechanism of C₂ under different pollution conditions remains incompletely understood and requires further investigation.

Over the past decades, the implementation of stringent air pollution control policies has led to a significant improvement in air quality across the PRD region (Bian et al., 2019; Geng et al., 2019; Zhang et al., 2019; Yan et al., 2020). Although VOC emissions increased during this period (Bian et al., 2019; Guo et al., 2024), a reduction in SOA was observed (Yan et al., 2020), which was likely related to the changes in ALWC and atmospheric oxidation capacity. Long-term measurements of DCA can provide critical insights into the SOA formation processes from high-pollution to low-pollution conditions. However, until now, such measurements remain scarce in the PRD. In this study, we conducted long-term measurements of aliphatic DCA (C₂, C₄–C₁₀), primary anthropogenic source markers, as well as SOA markers at a background station in the PRD during 2007–2018. Based on annual trends and correlation analysis, the driving factors of C₂ were identified. Then, a machine learning model was applied to quantify the impacts of anthropogenic precursor emissions, biogenic precursor emissions, aqueous-phase oxidation processes, and gas-phase oxidation processes on C₂ variability. Our results highlight an enhanced role of gas-phase oxidation in C₂ formation under cleaner conditions.

2 Materials and methods

2.1 Field sampling

Particulate matter with aerodynamic diameter less than 2.5 µm (PM_2.5) samples were collected on the rooftop of a five-storey building (approximately 30 m above the ground) in the campus of a middle school at Wanqingsha (WQS; 22.42° N, 113.32° E). The sampling site is in a small town in Guangzhou, which is located in the center of the PRD region (Fig. S1 in the Supplement). The influence of local anthropogenic emissions on this site was limited due to the low traffic and less industrial emissions from surrounding areas. During the summertime, southwesterly winds bring clean and moist air masses from the ocean. In contrast, northerly winds transport polluted and dry air masses from urban clusters during fall and winter. This makes WQS an ideal background station for investigating regional air pollution.

PM_2.5 samples were collected using prebaked (450°, 4 h) quartz filters (20.3 cm × 25.4 cm, QMA, Whatman, UK). Each sample lasted for 24 h using a high-volume air sampler (HVPM2.5, Tisch Environmental Inc., USA) at an airflow rate of 1.1 m³ min⁻¹. Field blank samples were also collected by placing the blank filter on the sampler for 10 min without turning on the sampler. In this study, a total of 462 PM_2.5 samples were collected mostly during the wintertime (October, November, and December of each year) from 2007 to 2018. The detailed information about sampling can be found in Table S2 in the Supplement. After the collection, each filter was wrapped in an aluminum foil, zipped in Teflon bags, and stored in a freezer (−20 °C) prior to analysis. The samples were analyzed within several months after collecting. Resch et al. (2023) reported that ambient aerosol samples stored at −20 °C for more than one month largely retained their molecular profiles. In addition, this study demonstrates that low molecular weight carboxylic acids are more stable than high molecular weight carboxylic acids during storage. Thus, the aliphatic DCA measured in this study are expected to remain largely stable during storage.

Gaseous pollutant data (e.g., NO₂ and O₃) and meteorological parameters (e.g., temperature, relative humidity) were obtained from an air quality monitoring station operated at WQS. The station was under maintenance during 2012–2013, therefore data for this period are unavailable. The surface net solar radiation (SR) and boundary layer height (BLH) data used in this study were obtained from the ERA5 reanalysis dataset provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) via the Copernicus Climate Data Store (CDS, https://cds.climate.copernicus.eu/datasets/, last access: 5 December 2025). The concentrations of PM_2.5 and its main components, as well as ALWC and pH, can be found in our previous study (He et al., 2025).

2.2 Chemical analysis

For analysis of oxalate (C₂${\mathrm{O}}_{\mathrm{4}}^{\mathrm{2}-}$), a punch (5.06 cm²) of the filter was extracted twice with 10 mL ultrapure Milli–Q water (18.2 MΩcm/25 °C) each for 15 min using an ultrasonic ice–water bath. The total water extracts (20 mL) were filtered through a 0.22 µm pore size filter and then stored in a pre-cleaned HDPE bottle. An ion-chromatography system (Metrohm, 883 Basic IC plus) was used to identify oxalate.

For organic tracers, a quarter of the filter was cut into small pieces and then mixed with 400 µL isotope-labeled internal standards, which included hexadecane-D₃₄, tetracosane-D₅₀, triacontane-D₆₂, naphthalene-D₈, acenaphthene-D₁₀, phenanthrene-D₁₀, chrysene-D₁₂, perylene-D₁₂, C₂₇aaa(20R)-cholestane-D₄, lauric acid-D₂₃, hexadecanoic acid-D₃₁, docosanoic acid-D₄₃, levoglucosan-¹³C₆, and phthalic acid 3,4,5,6-D₄. Samples were sonicated in an ice-water bath using mixed solvents, twice with a 1:1 ($v/v$) dichloromethane (DCM)/hexane solution, followed by twice with a 1:1 ($v/v$) DCM/methanol solution. The extraction solutions of each sample were combined, filtered, and concentrated to ∼ 2 mL, and then split into two parts for silylation and methylation, respectively. All derivatized solutions were analyzed by an Agilent gas chromatography-mass spectrometry (GC/MS, 7890/5975C) in the scan mode with an HP-5MS capillary column (30 m × 0.25 mm × 0.25 µm). The detailed information about sample treatment and GC-MS analysis was described by previous studies (Ding et al., 2012, 2014). In this study, we measured 9 dicarboxylic acids (succinic acid, glutaric acid, adipic acid, pimelic acid, suberic acid, azelaic acid, sebacic acid, phthalic acid, terephthalic acid), 5 hopanes (17α(H)-22,29,30-trisnorhopane, 17α(H),21β(H)-30-norhopane, 17α(H),21β(H)-30-hopane, 17α(H),21β(H)-22R-homohopane, and 17α(H),21β(H)-22S-homohopane), levoglucosan, octadecanoic acid, picene, 2,3-dihydroxy-4-oxopentanioic acid (DHOPA), and malic acid. Table S3 in the Supplement shows detailed information about target compounds and their internal standards. Due to the lack of a commercial standard, DHOPA was quantified using an isomer surrogate, citramalic acid. The analytical procedure of DHOPA is described elsewhere (Ding et al., 2017).

2.3 Quality assurance/quality control

We analyzed the field blanks and found that oxalate, succinic acid, octadecanoic acid, phthalic acid, and terephthalic acid were detected in trace amounts. Thus, the data reported in this study were corrected by corresponding field blanks. We also analyzed the spiked samples (n=7) to test the recoveries of analytical procedures. The results indicated that the recoveries were 95 %–103 % for oxalate, 100 %–121 % for levoglucosan, 66 %–118 % for dicarboxylic acid and octadecanoic acid, 89 %–113 % for hopanes, 76 %–97 % for picene, 88 %–94 % for citramalic acid. The method detection limits (MDLs) were 0.04 µg m⁻³ for oxalate, 6.7–40.5 ng m⁻³ for dicarboxylic acids, 0.07 ng m⁻³ for levoglucosan, 0.010–0.017 ng m⁻³ for hopanes, 0.18 ng m⁻³ for octadecanoic acid, 0.11 ng m⁻³ for picene, and 0.11 ng m⁻³ for citramalic acid.

Examining long-term trends of molecular markers is inherently challenging, as even minor changes in sampling or analytical protocols may introduce discontinuities. Consequently, it was necessary to assess the stability of both the sampling procedures and the analytical system over the entire study period. Table S4 in the Supplement presents the response factors (RFs) of all measured species derived from the annual calibration curves. The RF values remained highly consistent across years, indicating that neither instrument sensitivity nor analytical performance experienced significant drift during the measurement period.

Before data analysis, all data were manually inspected and the outliers (i.e., ${\mathrm{X}}_{\mathrm{75}\mathrm{\%}}+\mathrm{3}({\mathrm{X}}_{\mathrm{75}\mathrm{\%}}-{\mathrm{X}}_{\mathrm{25}\mathrm{\%}})$) were removed to rule out the influence of extreme concentrations on the overall trends. The change rates were calculated using the slopes derived from Theil–Sen regression and evaluated for statistical significance via the non-parametric Mann–Kendall test, providing a robust and reliable assessment of temporal variations.

2.4 Evaluate the impacts of various factors to oxalic acid variation by the machine learning model

Extreme Gradient Boosting (XGBoost), an advanced ensemble machine learning method based on gradient boosting decision trees, is known for its high computational efficiency, robust predictive performance (Chen et al., 2016) and thus has been applied in air pollutant research recently (Hou et al., 2022; Peng et al., 2023; Liu et al., 2025). In this study, XGBoost was employed to assess the impacts of various factors to oxalic acid variation. The implementation and Python package of XGBoost algorithm are publicly available online (https://github.com/dmlc/xgboost, last access: 5 December 2025). A total of 14 variables were used as input features to train the model, including levoglucosan, hopanes, octadecanoic acid, picene, terephthalic acid (tPh), O_x, photolysis frequencies of O₃ (JO¹D) and NO₂ (JNO₂), ALWC, pH, sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), temperature (Temp), solar radiation (SR), and relative humidity (RH). To avoid redundant and confounding explanations, the secondary organic molecular markers, such as DHOPA, phthalic acid (Ph), and malic acid, were excluded in the model training. They are influenced by VOC emissions and secondary oxidation processes, which are already represented by the factors mentioned above. Our results showed that there were great agreements between the observations and simulations for C₂ and other DCA (Fig. S2 in the Supplement), which indicated the model predictions were reliable.

To further explain the results of XGBoost, Python package of SHapley Additive ExPlanation (SHAP), a game theory-based approach proposed by Lundberg and Lee (2017), was used in this study (https://github.com/slundberg/shap, last access: 5 December 2025). The SHAP approach can quantify the impact of each feature to the model output while also accounting for the complex interdependencies among features. Thus, this approach provides deeper insights into the key drivers of oxalic acid formation. The fundamental principle of SHAP is to fairly allocate the prediction output among all input features by evaluating their marginal contributions across all possible feature subsets. For a given prediction f(x_i) generated by the model, the SHAP decomposes the prediction into the sum of feature attributions (SHAP values) and a baseline expectation (Eq. 1):

$$\begin{array}{}\text{(1)}& f\left(x_i\right)={\mathrm{\varnothing }}_{\mathrm{0}}(f,x)+{\sum }_{j=\mathrm{1}}^K{\mathrm{\varnothing }}_j(f,x_i)\end{array}$$

where ∅_j(f,x_i) is the SHAP value representing the impact of variable j on the prediction of the model for input x_i. The base value, ${\mathrm{\varnothing }}_{\mathrm{0}}(f,x)=E\left[f\right(x\left)\right]$, is the expected value of the model output over the data set. More details about this approach can be found in previous literature (Lundberg and Lee, 2017; Hou et al., 2022).

SHAP value represents the impact of a unit change in an individual variable on the predicted C₂ concentration, while holding all other variables constant. The sign of the SHAP value reflects the direction of the feature's effect on the target variable, while its magnitude quantifies the strength of that effect. For example, a positive (negative) SHAP value indicates that the change in a given feature increases (decreases) the predicted value relative to the model's baseline. The mean absolute SHAP (|SHAP|) values can remove the influence of directional effects that may cancel each other out when averaged, thereby providing a more robust measure of each feature's net contribution to model output. Using |SHAP| values facilitates a comprehensive ranking of feature importance based solely on contribution intensity, regardless of whether the effect is positive or negative. We further calculated the impact factor (IF) to evaluate the impact of individual factor (i) to oxalic acid variation as Eq. (2), where N is the number of samples:

$$\begin{array}{}\text{(2)}& \text{IF}\left(\mathit{\%}\right)={\displaystyle \frac{\mathrm{|}\text{SHAPvalue}\left(i\right)\mathrm{|}}{{\sum }_{j=\mathrm{1}}^N\mathrm{|}\text{SHAPvalue}\left(j\right)\mathrm{|}}}\end{array}$$

3 Results and discussion

3.1 Annual trends of the molecular tracers and aliphatic DCA in the PRD

Annual variations in the organic molecular markers, aliphatic DCA and meteorological parameters are summarized in Table S5 in the Supplement and presented in Fig. 1. Five molecular markers representing different primary sources were analyzed in this study, including levoglucosan, hopanes, octadecanoic acid, picene, and terephthalic acid. These compounds are commonly used as tracers for biomass burning (Giannoni et al., 2012), vehicle emissions (Riddle et al., 2007), cooking (Zhao et al., 2007), coal combustion (Yu et al., 2020), and waste incineration (Simoneit et al., 2005), respectively. Between 2007 and 2018, the concentrations of levoglucosan, hopanes and octadecanoic acid exhibited significant declining trends (−11 % yr⁻¹, −17 % yr⁻¹, and −7 % yr⁻¹, respectively; p < 0.01). In contrast, no clear trends were observed for picene and terephthalic acid, suggesting that control measures targeting biomass burning, vehicle emissions, and cooking were effective, while those for coal combustion and waste incineration may have been less effective or less implemented. Phthalic acid has been identified as a SOA tracer derived from naphthalene (Kleindienst et al., 2012), while DHOPA is a tracer for SOA formed from aromatic hydrocarbons (Ding et al., 2017). Due to the substantial anthropogenic sources of naphthalene and aromatic hydrocarbons, phthalic acid and DHOPA can be used as anthropogenic SOA (ASOA) markers. Malic acid is a typical secondary product formed through photooxidation of both anthropogenic and biogenic VOCs (AVOCs and BVOCs). However, a recent study estimated that malic acid produced through the oxidation of BVOCs (α-pinene and isoprene) was higher than that formed by AVOCs (toluene and naphthalene) (Sato et al., 2021). In addition, malic acid was also found to be strongly correlated with monoterpene tracers (R² = 0.87–0.95) in field measurements (Hu and Yu, 2013; Cheng et al., 2021). Given high BVOC emissions (Wang et al., 2021a) in the PRD region, large amount of malic acid was likely produced from biogenic precursors in this study, especially monoterpenes. Thus, we used malic acid to reflect the variations of SOA (BSOA). Although phthalic acid, DHOPA, and malic acid decreased from 51.9 ± 14.9, 1.85 ± 1.35, and 24.2 ± 19.4 ng m⁻³ to 16.7 ± 5.7, 1.05 ± 0.88, 5.9 ± 4.9 ng m⁻³, respectively, their declining trends were not statistically significant (p > 0.05). This indicated that the influence of reductions in emissions of anthropogenic organic precursors on SOA was limited.

Figure 1 Annual variations in different molecular markers (ng m⁻³) in the PRD during 2007 to 2018. The shaded area represents the 95 % prediction band.

Figure 2 (a) Annual variations in aliphatic DCA. (b) Annual variations in oxalic acid. Due to the absence of aliphatic DCA and oxalic acid measurements in 2009, their concentrations for that year are not presented.

Aliphatic DCA, a group of typical SOA molecular markers (Kawamura and Bikkina, 2016), were analyzed to further investigate the changes in SOA formation under decreasing pollutant levels. Although their concentrations in 2007 (864 ± 283 ng m⁻³) were much higher than those in 2018 (307 ± 122 ng m⁻³) (Fig. 2a), this downward trend was not statistically significant (p > 0.05), which was similar to the trends of other SOA markers. Unexpectedly, a rebound in the concentrations of molecular markers and aliphatic DCA was observed in 2013, which was also observed in PM_2.5 and its major component (He et al., 2025). This rebound was unlikely driven by meteorological variability, as inter-annual differences in key meteorological parameters were relatively small (Table S5). Similarly, another study covering multiple stations across the PRD also reported a post-2013 increase in PM_2.5 (Yan et al., 2020), supporting that the rebound was not attributable to newly emerging pollution sources near our sampling site. Unfortunately, until now, the underlying cause of this rebound remained unclear.

C₂ was the most abundant compound among aliphatic DCA, accounting for 80 %–91 %, followed by C₄ (4 %–13 %), and C₉ (1 %–4 %). Therefore, the overall trend of aliphatic DCA was primarily driven by C₂ (Fig. 2b), and subsequent discussions will focus on C₂. Its concentration declined from 692 ± 243 (2007) to 274 ± 114 (2018), but did not exhibit a clear trend (p > 0.05). Carbon monoxide (CO) can be used as a normalization tracer to assess the influence of atmospheric dilution. As shown in Fig. S3 in the Supplement, the temporal trends of DCA and C₂ normalized by CO are consistent with their original trends, indicating that atmospheric dilution had a limited influence on their observed variations. To further explore the changes of SOA formation under different pollution conditions, our samples were divided into five categories according to interim targets recommended by the Worle Health Organization (WHO) in 2021 (World Health Organization, 2021): IT0 ( PM_2.5 > 75 µg m⁻³), IT1 (75 µg m⁻³ > PM_2.5 > 50 µg m⁻³), IT2 (50 µg m⁻³ > PM_2.5 > 37.5 µg m⁻³), IT3 (37.5 µg m⁻³ > PM_2.5 > 25 µg m⁻³), and IT4 (PM_2.5 < 25 µg m⁻³). We found that the molecular markers and C₂ decreased significantly (p < 0.01) from IT0 to IT4 (Table S6 in the Supplement). However, the ratio of C₂ to PM_2.5 (${\mathrm{C}}_{\mathrm{2}}/{\mathrm{PM}}_{\mathrm{2.5}}$) increased from 6.8 × 10⁻³ to 10.3 × 10⁻³ (p < 0.01, Fig. S4 in the Supplement), suggesting that the relative importance of SOA increased as pollution levels decreased.

3.2 Sources and driving factors of oxalic acid

The ratio of C₂ to C₄ (${\mathrm{C}}_{\mathrm{2}}/{\mathrm{C}}_{\mathrm{4}}$) has been used as an indicator for the aging degree of DCA and C₂ (Bikkina et al., 2015; Bikkina et al., 2017; Zhao et al., 2020; Meng et al., 2023), because the photochemical degradation of C₄ can lead to C₃, and C₃ can be photochemically oxidized into C₂ via intermediates (Kawamura and Bikkina, 2016). In this study, ${\mathrm{C}}_{\mathrm{2}}/{\mathrm{C}}_{\mathrm{4}}$ (21.9 ± 18.0) was significantly higher than those reported for primary sources, such as fossil fuel emission (4.1) (Kawamura and Kaplan, 1987) and biomass burning (1.6) (Cao et al., 2017). This indicated that C₂ was highly aged and mainly formed through secondary processes rather than being directly emitted from primary sources in this study. The ratio was comparable to those observed at a rural site of the PRD (15.4–21.0) (Xu et al., 2022), but substantially higher than those reported in urban areas of North China, such as Jinan (8.4 ± 3.4) (Meng et al., 2023), Tianjin (4.5–6.1) (Pavuluri et al., 2021), Xi'an (4.6–8.0) (Wang et al., 2012), and Beijing (4.7–8.6) (Zhao et al., 2018). This was likely related to elevated temperature (above 20 °C), solar radiation (95.3–161.3 W m⁻²), and relative humidity (∼ 60 %) in the PRD (Table S5), which led to a higher degree of aerosol aging.

Figure 3 Correlations between C₂ and primary anthropogenic source markers, SOA markers, indicators for gas- and aqueous-phase oxidation, as well as meteorological parameters. Blank cells indicate no significant correlations. One asterisk, two asterisks denote p value < 0.05, 0.01, respectively. Due to the unavailability of O_x data in 2012 and 2013, correlation analysis was not conducted for these 2 years.

The results of correlation analysis are presented in Fig. 3. The correlation coefficients between C₂ and SOA tracers (phthalic acid, DHOPA, and malic acid) were relatively higher (r = 0.58, 0.41, and 0.57, respectively; p < 0.01), further supporting that C₂ was primarily formed via secondary oxidation processes. Meng et al. (2023) reported an unexpected enhancement of C₂ during the COVID-19 pandemic, when anthropogenic emissions were substantially reduced. This reflected limited influence of reductions in anthropogenic organic precursors on formation of C₂. Similarly, we found that anthropogenic emissions experienced substantial reductions during our campaign period (discussed in Sect. 3.1), while C₂ did not show a corresponding significant decreasing trend. Although strong correlations between C₂ and primary anthropogenic source markers were observed in certain individual years, the correlations remained weak across the entire dataset. These findings implied that the changes in anthropogenic emissions were not the driving factors for C₂ formation in this study. Because our field measurements were conducted in the same season each year (from October to December), the inter-annual differences in meteorological conditions were small. This resulted in consistently weak correlations observed between C₂ and meteorological parameters such as temperature, SR, and RH. Therefore, we concluded that the changes in meteorology were too small to be the driving factors for C₂ formation.

Notably, the correlations between C₂ and malic acid strengthened progressively with the reductions in anthropogenic emissions (from 2007 to 2018). This trend became more apparent when the data were categorized by pollution levels, with the correlation coefficients increased from 0.33 (IT0) to 0.72 (IT4) (Table S7 in the Supplement). As discussed previously, malic acid can be produced by photooxidation of both anthropogenic and biogenic precursors. However, no corresponding increasing trends were observed in the correlations between C₂ and ASOA tracers (Ph and DHOPA), supporting that anthropogenic precursors were not the dominant source of malic acid in this study. Thus, these results indicated that the relative contributions of biogenic sources to SOA become more important under cleaner conditions.

ALWC not only regulates the gas-particle partitioning of semi-volatile VOCs and their reaction rates by acting as a medium (Nenes et al., 2021), but also serves as a nucleophile that participates in reactive uptake of SOA intermedium (Zhang et al., 2022b). Aerosol pH plays a crucial role in governing acid-catalyzed reactions during aqueous-phase processing (Cooke et al., 2024). In addition, sulfate is also an important secondary product formed through aqueous oxidation (Liu et al., 2021). O_x, a proxy of atmospheric oxidants, facilitates secondary photochemical oxidation of VOCs. J(O¹D) and J(NO₂) represent photolysis frequencies of O₃ and NO₂ (Ehhalt and Rohrer, 2000). Accordingly, ALWC, pH, and sulfate were employed as indicators for aqueous-phase oxidation , while O_x, J(O¹D), and J(NO₂) were used as indicators for gas-phase oxidation. In this study, C₂ showed strong correlations with ALWC (r = 0.50, p < 0.01), sulfate (r = 0.62, p < 0.01), and O_x (r = 0.64, p < 0.01) across the whole datasets, suggesting secondary oxidation processes were the dominant drivers of C₂ variability between 2007 and 2018.

3.3 Gas-phase and aqueous-phase pathways of oxalic acid formation

Atmospheric oxidation of VOCs, such as aromatics, isoprene, and monoterpenes, with ozone and/or hydroxyl radicals produces semi-volatile gaseous organic precursors (e.g., mGly and Gly). These α-dicarbonyls partition into ALWC-enriched aerosols and the subsequent aqueous-phase reactions produce Pyr and ωC₂, which are further oxidized to form C₂ (Ervens et al., 2008; Bikkina et al., 2017). In addition, several studies have also demonstrated that C₂ can be formed through the photochemical degradation of longer-chain DCA (e.g., C₃ and C₄) in the gas phase (Kawamura and Bikkina, 2016; Zhao et al., 2020; Meng et al., 2023). Due to the reductions in sulfate and nitrate (the highly hygroscopic compounds in PM_2.5), ALWC experienced a fast downward trend between 2007 and 2012 (Table S5, Fig. S5 in the Supplement). This hindered the partitioning of semi-volatile gaseous organic precursors from the gas phase into the aqueous phase, thereby slowing the aqueous-phase production of C₂. As a result, an obvious decline in C₂ concentrations was observed during 2007–2012 (Fig. 2). Because ALWC exhibited minor fluctuations between 2013 and 2018, the influence of aqueous-phase pathway on C₂ became weaker.

Figure 4 The correlation between C₂ and ALWC, as well as O_x. The values in brackets indicate the 95 % confidence intervals (CIs) of the correlation coefficients. One, two asterisks denote p values less than 0.05, 0.01, respectively. With decreasing pollution levels, the correlation between C₂ and ALWC weakens, whereas that between C₂ and O_x strengthens.

As presented in Fig. 3, the correlations between C₂ and ALWC weakened, whereas that between C₂ and O_x strengthened from 2007 to 2018. It is noteworthy that no significant correlations were observed between C₂ and O_x in 2007 and 2010, when ALWC concentrations were high. In contrast, a strong correlation (r = 0.62, p < 0.01) emerged in 2008, coinciding with a sharp decrease in ALWC (Fig. S5). This suggested that the gas-phase formation pathway of C₂ might be enhanced when ALWC was low. The pattern became more pronounced with decreasing pollution levels (Fig. 4). From IT0 to IT4, ALWC decreased from 20.9 ± 11.0 to 7.2 ± 3.0 µg m⁻³ (Table S6). The Pearson's r values between C₂ and ALWC dropped from 0.48 to 0.15, whereas that between C₂ and O_x increased from 0.28 to 0.68. The opposite trends implied the roles of gas-phase and aqueous-phase oxidation in C₂ formation might change. However, the precision and stability of Pearson's r values are strongly influenced by sample size when the variations appear small or when sample sizes differ among groups. Therefore, the differences in Pearson's r values do not necessarily imply statistically significant changes, especially when they are very close (IT1–IT4). To assess the statistical significance of these differences, we compared correlation coefficients between groups using the method described in Sect. S1 in the Supplement. As shown in Tables S8 and S9 in the Supplement, significant differences in the C₂-O_x correlation were observed only between IT0 and the other pollution levels. For the C₂-ALWC correlation, a significant difference was found only between IT0 and IT4. Given that IT1–IT4 represents a continuous evolution of atmospheric conditions, rather than discrete and independent regimes, large differences in correlation coefficients among these categories are not expected. Although the correlation between C₂ and sulfate was strong, it did not show the similar trends as that between C₂ and ALWC. In contrast, the correlations between C₂ and primary anthropogenic markers remained generally weak across all pollution categories (Table S7), indicating that changes in anthropogenic emissions exert only limited influence on C₂ variations. Therefore, the significant and opposite changes in correlations of C₂ with O_x and ALWC between high pollution level (IT0) and low pollution level (IT4) suggested a shift in the dominant C₂ formation pathway from aqueous-phase oxidation to gas-phase photochemical oxidation under lower pollution conditions.

The Pearson correlation coefficients reflect the strength of linear associations but do not directly represent the quantitative contribution of each factor. Thus, the lower Pearson's r values between C₂ and ALWC than those between C₂ and O_x under IT1-IT4 do not necessarily imply a smaller impact of aqueous-phase pathways compared with gas-phase pathways. Here, correlation analysis is used primarily to identify potential drivers and their changing patterns, while a machine learning model is further applied to quantify their impacts to C₂ variations.

Figure 5 (a) Bar plot of the mean |SHAP| values representing the overall importance of each feature in predicting C₂ concentrations. (b) Beeswarm plot of individual SHAP values for each feature across all samples. Red (blue) represents high (low) value in each feature. Positive (negative) SHAP values indicate that the feature contributes to an increase (decrease) in the C₂ prediction.

3.4 Evaluate the impacts of different factors on oxalic acid formation

We identified ALWC and O_x as the dominant drivers for C₂ variation based on annual variations and linear correlation analysis. Although pH has been suggested to play an important role in heterogeneous uptake and acid-catalyzed reactions (Xu et al., 2016; Cooke et al., 2024), it did not exhibit a strong correlation with C₂ (r = −0.24, p < 0.01). This was likely due to their non-linear relationship. In addition, correlation analysis is unable to quantify the influences of changes in each factor. Hence, we used a machine learning model to further evaluate the impacts of changes in individual factor on C₂ variation during our study period. The details about the model used in this study can be seen in Sect. 2.4 and elsewhere (Chen et al., 2016; Lundberg and Lee, 2017).

The rationality for selecting the variables used to train the model need to be clarified to ensure the reliability of the results. Levoglucosan, hopanes, octadecanoic acid, picene, and tPh serve as source-specific molecular markers for biomass burning, vehicle emission, cooking, coal combustion, and waste incineration, respectively. These species are used to represent changes in AVOC emissions. As two of the most important BVOCs globally, isoprene emission is highly dependent on temperature and solar radiation, while monoterpenes emission is sensitive to temperature (Guenther et al., 1993). Their emissions rate can be estimated using Eqs. (3)–(5) and Eq. (6), respectively:

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{i}}=I_{\mathrm{s}}\cdot C_{\mathrm{L}}\cdot C_{\mathrm{T}}\text{(4)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{L}}={\displaystyle \frac{\mathit{\alpha }{c}_{\mathrm{L}\mathrm{1}}L}{\sqrt{\mathrm{1}+{\mathit{\alpha }}^{\mathrm{2}}{L}^{\mathrm{2}}}}}\text{(5)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{T}}={\displaystyle \frac{\text{exp}\frac{c_{\mathrm{T}\mathrm{1}}(T-T_{\mathrm{S}})}{R{T}_{\mathrm{S}}T}}{\mathrm{1}+\text{exp}\frac{c_{\mathrm{T}\mathrm{2}}(T-T_{\mathrm{M}})}{R{T}_{\mathrm{S}}T}}}\end{array}$$

where E_i is isoprene emission rate at a temperature T (K) and photosynthetically active radiation (PAR) flux L ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), I_s is isoprene emission rate at a standard temperature T_S and a standard PAR flux (1000 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). α = 0.0027 and c_L1 = 1.066 are empirical coefficients determined by measurements. L can be calculated as multiplying solar radiation (W m⁻²) by photon flux efficacy (1.86 µmol J⁻¹). R is a constant −8.314 $\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{mol}}^{-\mathrm{1}}$, and c_T1 = 95 000 J mol⁻¹, c_T2 = 230 000 J mol⁻¹, and T_M = 314 K are empirical coefficients estimated by measurements.

$$\begin{array}{}\text{(6)}& E_{\mathrm{m}}=M_{\mathrm{S}}\cdot \exp \left(\mathit{\beta }\right(T-T_{\mathrm{S}}\left)\right)\end{array}$$

where E_m is monoterpenes emission rate at temperature T (K), M_s is monoterpenes emission rate at a standard temperature T_S, β (K⁻¹) is an empirical coefficient ranging from 0.057 to 0.144 K⁻¹. In addition, inadequate moisture can have significantly decreased stomatal conductance and photosynthesis (Guenther et al., 2006). Therefore, temperature, solar radiation, and RH are important factors determining BVOC emissions. As discussed previously, O_x, J(O¹D), and J(NO₂) can be employed as indicators of gas-phase oxidation, while ALWC, pH, and sulfate are used as indicators of aqueous-phase processes. To avoid redundant and confounding explanations, the secondary organic molecular markers, such as DHOPA, Ph, and malic acid, were excluded from the model training. These species are influenced by both VOC emissions and secondary oxidation processes, which are already represented by the factors mentioned above.

The feature importance is presented in Fig. 5a. O_x, sulfate, and J(O¹D), which represent secondary oxidation processes, exhibited the three highest |SHAP| values, indicating their dominant impacts on C₂ variation. Although pH and ALWC exhibited relatively high feature importance among all variables, their |SHAP| values were lower than sulfate. This is because pH and ALWC in this study was calculated by a thermodynamic equilibrium model, ISORROPIA II (Nenes et al., 1998), in which sulfate plays a crucial role and partly reflects variations in both pH and ALWC. In contrast, the feature importance of anthropogenic emission markers and meteorological parameters were relatively lower, suggesting that the impacts of AVOC and BVOC emissions were smaller compared to that of secondary processes. As shown in Fig. 5b, O_x, sulfate, and J(O¹D) exhibited obviously positive correlations with their SHAP values, indicating that higher values of these variables contributed to increases in C₂ concentrations. However, pH showed a negative correlation with its SHAP values, suggesting that lower pH levels were associated with higher C₂ concentrations. Notably, the influence of extremely low pH on C₂ formation appeared to be more pronounced.

Figure 6 (a) Impacts of changes in each variable on C₂ variation during the whole study period. (b) Impact factor of individual variable under different pollution conditions.

To further quantify the impacts of changes in all factors on C₂, IF (discussed in Sect. 2.4) was calculated and presented in Fig. 6. O_x accounted for the highest IF (35 %), followed by sulfate (24 %) and J(O¹D) (9 %). All factors were classified into four groups according to their representativeness mentioned before: (1) AVOC emissions (levoglucosan, hopanes, octadecanoic acid, picene, and tPh); (2) BVOC emissions (Temp, SR, and RH); (3) gas-phase oxidation pathways (O_x, J(O¹D), and J(NO₂)); (4) aqueous-phase oxidation pathways (ALWC, pH, and sulfate). Due to the minor fluctuations of meteorological conditions in each year, the impacts of changes in BVOC emissions on C₂ were small (7 %). Although AVOC emissions showed an obvious decreasing trend over the study period, the impacts of these changes (14 %) were significantly lower than that of gas-phase oxidation processes (45 %) and aqueous-phase oxidation processes (34 %). The results were consistent with correlation analysis, underscoring the dominant role of secondary oxidation processes in C₂ formation.

The IF values for each variable are presented in Table S10 in the Supplement. From IT0 to IT4, IF values for gas-phase oxidation processes increased from 37 % to 55 %, whereas those for aqueous-phase oxidation processes decreased from 42 % to 30 % (Fig. 6b). Meanwhile, IF values for AVOC (10 %–15 %) and BVOC emissions (5 %–8 %) remained at a low and stable level. These findings indicated that the gas-phase oxidation pathway became increasingly important as pollution levels decreased. A possible explanation is that under cleaner conditions, lower ALWC levels would suppress the partitioning of semi-volatile C₂ precursors (e.g., Gly and mGly) from the gas-phase into the particle-phase (Hu et al., 2022). In addition, less ALWC participates in heterogeneous reactions of SOA intermedium as a nucleophile. Thus, their aqueous-phase pathway was hindered, and more C₂ was formed via photochemical degradation of longer-chain DCA (Kawamura and Bikkina, 2016; Meng et al., 2023). This indicated the growing importance of gas-phase oxidation processes in the formation of C₂ and SOA under cleaner conditions. Although O_x did not exhibit a clear trend at our measurement station, Cao et al. (2024) reported a rapid increase in O₃ concentration across the PRD region over the past decade. This may promote SOA formation through enhanced gas-phase oxidation pathways. Therefore, coordinated control of VOCs and NO_x should be emphasized (Wang et al., 2021b) in the future to reduce ozone pollution and further mitigate SOA formation.

4 Conclusions

In this study, long-term measurements were conducted at a regional background site in the PRD region during 2007–2018. Primary and secondary molecular markers were measured, including five primary anthropogenic source molecular markers (levoglucosan, hopanes, octadecanoic acid, picene, and terephthalic acid), three SOA markers (phthalic acid, DHOPA, and malic acid), and a group of general SOA indicators (aliphatic DCA). C₂, the most abundant compound among aliphatic DCA, was investigated to identify the key drivers for SOA variability. Previous studies have attributed elevated C₂ to higher emissions of its organic precursors (Cao et al., 2017; Ding et al., 2021; Yu et al., 2021) and emphasized the dominance of aqueous-phase pathways in C₂ formation (Yu et al., 2005; Hilario et al., 2021). However, our results showed that reductions in anthropogenic organic precursors had a limited influence on C₂ formation. In addition, we found a growing importance of gas-phase oxidation to C₂ and SOA formation under cleaner conditions.

Between 2007 and 2018, C₂ and other SOA molecular markers did not show clear trends, despite substantial reductions in anthropogenic emissions, such as biomass burning (−11 % yr⁻¹), vehicle emissions (−17 % yr⁻¹), and cooking emissions (−7 % yr⁻¹). Correlation analysis revealed that ALWC and O_x were the dominant drivers for C₂ variations. In addition, the correlation between malic acid and C₂ became stronger from high to low pollution conditions (from IT0 to IT4), underscoring an increasing relative contribution of BSOA in SOA. We applied a machine learning model to further quantify the impacts of changes in AVOC emissions, BVOC emissions, aqueous-phase oxidation processes, and gas-phase oxidation processes. The model results showed that, although AVOC emissions showed an obvious decreasing trend over the study period, the impacts of these changes (14 %) were significantly lower than that of gas-phase oxidation processes (45 %) and aqueous-phase oxidation processes (34 %). With decreasing pollutant levels, the impacts of gas-phase oxidation processes increased from 37 % to 55 %, whereas that of aqueous-phase oxidation processes declined from 42 % to 30 %. This probably resulted from the low ALWC hindered partitioning of semi-volatile precursors (e.g., Gly and mGly) from gas phase to particle phase, thereby suppressing their aqueous-phase reactions and enhancing the gas-phase photochemical degradation of longer-chain DCA (e.g., C₃ and C₄). Meanwhile, the impacts of changes in AVOC (10 %–15 %) and BVOC emissions (5 %–8 %) were minor. Our results provide critical insights that improve the accuracy of atmospheric modeling of C₂, advance understanding of SOA formation mechanisms under different pollution levels, and support the formulation of more targeted pollution mitigation strategies in the future. Meanwhile, we acknowledge there are several limitations in this study. First, our measurements were mainly conducted in wintertime, which may not represent summertime conditions when photochemical activity is higher. Second, there are other factors influencing the formation of aqueous-phase products other than pH and ALWC, similar to gas-phase products. However, due to unavailability of related data in this study, such as transition metals and hydroxyl radical in aqueous phase, we were unable to quantify their impacts on variations of C₂, which may introduce uncertainties.

The continuous increase of O₃ over the past decade in the PRD poses a growing challenge to further improvements in particulate pollution. As gas-phase oxidation processes become more important in SOA formation under a cleaner environment, coordinated control of both NO_x and VOCs is essential for mitigating O₃ pollution and achieving further reductions in SOA in the PRD region.